[GRASS-SVN] r65146 - in grass-addons/grass7/raster/r.green: . r.green.hydro r.green.hydro/libhydro r.green.hydro/r.green.hydro.closest r.green.hydro/r.green.hydro.delplants r.green.hydro/r.green.hydro.economic r.green.hydro/r.green.hydro.legal r.green.hydro/r.green.hydro.optimal r.green.hydro/r.green.hydro.recommended r.green.hydro/r.green.hydro.structure r.green.hydro/r.green.hydro.theoretical
svn_grass at osgeo.org
svn_grass at osgeo.org
Sun Apr 26 23:36:52 PDT 2015
Author: Giulia
Date: 2015-04-26 23:36:52 -0700 (Sun, 26 Apr 2015)
New Revision: 65146
Added:
grass-addons/grass7/raster/r.green/r.green.hydro/
grass-addons/grass7/raster/r.green/r.green.hydro/Makefile
grass-addons/grass7/raster/r.green/r.green.hydro/libhydro/
grass-addons/grass7/raster/r.green/r.green.hydro/libhydro/Makefile
grass-addons/grass7/raster/r.green/r.green.hydro/libhydro/__init__.py
grass-addons/grass7/raster/r.green/r.green.hydro/libhydro/basin.py
grass-addons/grass7/raster/r.green/r.green.hydro/libhydro/optimal.py
grass-addons/grass7/raster/r.green/r.green.hydro/libhydro/plant.py
grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.closest/
grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.closest/Makefile
grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.closest/r.green.hydro.closest.html
grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.closest/r.green.hydro.closest.py
grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.delplants/
grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.delplants/Makefile
grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.delplants/r.green.hydro.delplants.html
grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.delplants/r.green.hydro.delplants.py
grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.economic/
grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.economic/Makefile
grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.economic/r.green.hydro.economic.html
grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.economic/r.green.hydro.economic.py
grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.html
grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.legal/
grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.legal/Makefile
grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.legal/r.green.hydro.legal.html
grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.legal/r.green.hydro.legal.py
grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.optimal/
grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.optimal/Makefile
grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.optimal/availablestreams.jpg
grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.optimal/potentialplants.jpg
grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.optimal/r.green.hydro.optimal.html
grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.optimal/r.green.hydro.optimal.py
grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.recommended/
grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.recommended/Makefile
grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.recommended/r.green.hydro.recommended.html
grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.recommended/r.green.hydro.recommended.py
grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.structure/
grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.structure/Makefile
grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.structure/plantstructure.png
grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.structure/potentialplants.jpg
grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.structure/r.green.hydro.structure.html
grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.structure/r.green.hydro.structure.py
grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.theoretical/
grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.theoretical/Makefile
grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.theoretical/r.green.hydro.theoretical.html
grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.theoretical/r.green.hydro.theoretical.py
Modified:
grass-addons/grass7/raster/r.green/Makefile
Log:
Add r.green.hydro modules
Modified: grass-addons/grass7/raster/r.green/Makefile
===================================================================
--- grass-addons/grass7/raster/r.green/Makefile 2015-04-27 03:59:19 UTC (rev 65145)
+++ grass-addons/grass7/raster/r.green/Makefile 2015-04-27 06:36:52 UTC (rev 65146)
@@ -2,7 +2,9 @@
PGM=r.green
-SUBDIRS = r.green.biomassfor
+SUBDIRS = r.green.biomassfor \
+ r.green.hydro
+
include $(MODULE_TOPDIR)/include/Make/Dir.make
Added: grass-addons/grass7/raster/r.green/r.green.hydro/Makefile
===================================================================
--- grass-addons/grass7/raster/r.green/r.green.hydro/Makefile (rev 0)
+++ grass-addons/grass7/raster/r.green/r.green.hydro/Makefile 2015-04-27 06:36:52 UTC (rev 65146)
@@ -0,0 +1,19 @@
+MODULE_TOPDIR = ../../..
+
+PGM=r.green.hydro
+
+SUBDIRS = r.green.hydro.theoretical \
+ r.green.hydro.legal \
+ r.green.hydro.recommended \
+ r.green.hydro.economic \
+ r.green.hydro.closest \
+ r.green.hydro.optimal \
+ r.green.hydro.structure \
+ r.green.hydro.delplants \
+ libhydro
+
+include $(MODULE_TOPDIR)/include/Make/Dir.make
+
+default: parsubdirs htmldir
+
+install: installsubdirs
Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/Makefile
___________________________________________________________________
Added: svn:mime-type
+ text/x-makefile
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.green/r.green.hydro/libhydro/Makefile
===================================================================
--- grass-addons/grass7/raster/r.green/r.green.hydro/libhydro/Makefile (rev 0)
+++ grass-addons/grass7/raster/r.green/r.green.hydro/libhydro/Makefile 2015-04-27 06:36:52 UTC (rev 65146)
@@ -0,0 +1,25 @@
+MODULE_TOPDIR = ../../../..
+
+include $(MODULE_TOPDIR)/include/Make/Other.make
+include $(MODULE_TOPDIR)/include/Make/Python.make
+
+MODULES = basin plant optimal
+
+
+LIBGREENDIR = $(ETC)/r.green
+ETCDIR = $(LIBGREENDIR)/libhydro
+
+
+PYFILES := $(patsubst %,$(ETCDIR)/%.py,$(MODULES) __init__)
+PYCFILES := $(patsubst %,$(ETCDIR)/%.pyc,$(MODULES) __init__)
+
+default: $(PYFILES) $(PYCFILES) $(ETCDIR)/__init__.py $(ETCDIR)/__init__.pyc
+
+$(ETCDIR):
+ $(MKDIR) $@
+
+$(ETCDIR)/%: % | $(ETCDIR)
+ $(INSTALL_DATA) $< $@
+
+install:
+ cp -r $(ETCDIR) $(INST_DIR)
Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/libhydro/Makefile
___________________________________________________________________
Added: svn:mime-type
+ text/x-makefile
Added: svn:eol-style
+ native
Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/libhydro/__init__.py
___________________________________________________________________
Added: svn:mime-type
+ text/x-python
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.green/r.green.hydro/libhydro/basin.py
===================================================================
--- grass-addons/grass7/raster/r.green/r.green.hydro/libhydro/basin.py (rev 0)
+++ grass-addons/grass7/raster/r.green/r.green.hydro/libhydro/basin.py 2015-04-27 06:36:52 UTC (rev 65146)
@@ -0,0 +1,721 @@
+#!/usr/bin/env python
+# -- coding: utf-8 --
+#
+############################################################################
+#
+# AUTHOR(S): Giulia Garegnani
+# PURPOSE: Definition of the object basin and functions
+# COPYRIGHT: (C) 2014 by the GRASS Development Team
+#
+# This program is free software under the GNU General Public
+# License (>=v2). Read the file COPYING that comes with GRASS
+# for details.
+#
+#############################################################################
+#
+
+# import system libraries
+import os
+
+from scipy import integrate
+import numpy as np
+import itertools
+#import pdb
+
+
+from grass.pygrass.vector import VectorTopo
+from grass.script import core as gcore
+from grass.pygrass.messages import get_msgr
+from grass.script import mapcalc
+from grass.pygrass.utils import set_path
+
+set_path('r.green', 'libhydro', '..')
+set_path('r.green', 'libgreen', os.path.join('..', '..'))
+# finally import the module in the library
+from libgreen.utils import dissolve_lines
+from libgreen.utils import raster2numpy
+from libgreen.utils import raster2compressM
+from libgreen.utils import get_coo
+
+
+def discharge_sum(list_basin, list_ID):
+ """
+ sum discharge of a list of basins
+ :param list_basin: list of
+ the object Basin
+ :type list_basin: list
+ :param list_ID: list of ID to sum the discharge
+ :type list_ID: list
+ """
+ discharge_tot = 0
+ for ID in list_ID:
+ discharge_tot = discharge_tot + list_basin[ID].discharge_own
+ return discharge_tot
+
+
+def E_hydro(delta, Q):
+ """
+ compute the energy in kW
+ :param Q
+ discharge m3/s
+ :param delta
+ gross head m
+ """
+ if isinstance(delta, float) or isinstance(delta, int):
+ if (isinstance(Q, float) or isinstance(Q, int) or
+ isinstance(Q, np.float32)):
+ return max((delta * 9.81 * Q), 0)
+ else:
+ return 0.0
+ else:
+ return 0.0
+
+
+def check_new(b, lista, new):
+ """
+ add the ID of the upper basins with a recursive
+ function in order to add the new ID upper to the added
+ ID
+ """
+# ID_basin = b.ID_all[indice]
+# for new in lista[ID_basin].up:
+# b.ID_all.add(new)
+# indice += 1
+# if (indice < len(b.ID_all)):
+# return check_new(b, lista, indice)
+# else:
+# return
+ if not new:
+ return
+ else:
+ b.ID_all = b.ID_all.union(new)
+ for ID in new:
+ check_new(b, lista, lista[ID].up)
+
+
+class Basin(object):
+
+ def __init__(self, ID, h_mean=None, h_closure=0, discharge_own=None,
+ discharge_tot=None, E_own=None, ID_all=None,
+ E_up=None, area=None, length=None, up=None):
+ """Instantiate a Basin object.: Define the class Basin in order
+ to create a list
+ of object basins with own discharge, mean elevation, closure elevation,
+ and the ID of the all the up basins
+ usefull in order to evaluate the total discharge
+ without running r.watershed twice. The aim os to have an object Basin
+ and a dictionary of basins with key equal to ID in order to
+ manage several operations with the basins.
+
+ :param ID: specify the ID basin
+ :type ID:
+ :param h_mean: mean elevation of the basin
+ :type h_mean: int or float
+ :param h_closure: elevation at the closure point of the basin
+ :type h_closure: int or float
+ :param discharge_own: discharge of the single basin withot the
+ share due to the upper basins
+ :type discharge_own: int or float
+ :param discharge_tot: total discharge coming from upper basins, can be
+ computed thanks to ID_all and the method check_ID_up, plus the own
+ discharge. It is the total outflow
+ :type discharge_tot: int or float
+ :param E_own: value of potential power in the basin
+ :type E_own: int or float
+ :param up: list of ID of the closest basins (maximum 3)
+ :type up: list
+ :param ID_all: list of ID of all the upper basins that generate
+ the total inflow
+ :type ID_all: list
+ :param E_up: dict with key the ID of the upper basin and value the
+ potential power computed as function of the total
+ discharge and the difference between h_closure of
+ the upper basin and h_closure of the current basin
+ :type E_up: dict
+ >>> B1=Basin(1,h_mean=20, h_closure=10,discharge_own=5)
+ >>> B1.E_own
+ 0.49050000000000005
+ >>> B2=Basin(2,h_mean=15, h_closure=10,discharge_own=7)
+ >>> B3=Basin(3,h_mean=7, h_closure=5,discharge_own=10)
+ >>> B3.add_E_up(B1.ID)
+ >>> B3.E_up[1]=E_hydro(B1.h_closure-B3.h_closure, B1.discharge_own)
+ >>> E2 = E_hydro(B2.h_closure-B3.h_closure, B2.discharge_own)
+ >>> B3.add_E_up(B2.ID,E2)
+ 0.34335000000000004
+ >>> B3.E_up
+ {1: 0.24525000000000002, 2: 0.34335000000000004}
+ >>> B3.E_up.keys()[0]
+ 1
+ >>> B3.up = set([B1.ID, B2.ID])
+ >>> B3.up
+ set([1, 2])
+ >>> B3.up.remove(2)
+ >>> B3.up
+ set([1])
+ >>> B3.up.add(4)
+ >>> B3.up
+ set([1, 4])
+ >>> B0 = Basin(ID=0, up=set([3, 1]), discharge_own=2)
+ >>> B1 = Basin(ID=1, up=set([2, 6]), discharge_own=9)
+ >>> B2 = Basin(ID=2, discharge_own=3)
+ >>> B3 = Basin(ID=3, up=set([4, 5]), discharge_own=4)
+ >>> B4 = Basin(ID=4, discharge_own=3)
+ >>> B5 = Basin(ID=5, discharge_own=2)
+ >>> B6 = Basin(ID=6, up=set([7, 8]), discharge_own=6)
+ >>> B7 = Basin(ID=7, discharge_own=1)
+ >>> B8 = Basin(ID=8, discharge_own=6)
+ >>> B_tot = {0: B0, 1: B1,2: B2,3: B3,4: B4,5: B5,6: B6}
+ >>> B_tot.update({7: B7, 8: B8})
+ >>> B0.check_ID_up(B_tot)
+ >>> set([3, 1, 4, 5, 2, 6, 7, 8]).difference(B_tot[0].ID_all)
+ set([])
+ >>> B5.check_ID_up(B_tot)
+ >>> B3.check_ID_up(B_tot)
+ >>> B3.ID_all
+ set([4, 5])
+ >>> discharge_sum(B_tot,B_tot[0].ID_all)
+ 34
+ >>> temp=discharge_sum(B_tot,B_tot[0].ID_all)+B_tot[0].discharge_own
+ >>> B_tot[0].discharge_tot=temp
+ >>> B_tot[0].discharge_tot
+ 36
+
+ @author: GGaregnani
+ """
+ self.ID = ID
+ self.h_mean = h_mean
+ self.h_closure = h_closure
+ self.discharge_own = discharge_own
+ self.discharge_tot = discharge_tot
+ self.area = area
+ self.length = length
+
+ if ((E_own is None) and (h_mean is not None) and
+ (h_closure is not None)):
+
+ delta = (self.h_mean-self.h_closure)
+ self.E_own = E_hydro(delta, self.discharge_own)
+ else:
+ self.E_own = E_own
+
+ if up is None:
+ self.up = set()
+ elif isinstance(up, set):
+ self.up = up
+ else:
+ raise TypeError("ID_all must be a set")
+
+ if ID_all is None:
+ self.ID_all = up
+ elif isinstance(ID_all, set):
+ self.ID_all = ID_all
+ else:
+ raise TypeError("ID_all must be a set")
+
+ if E_up is None:
+ self.E_up = {}
+ elif isinstance(E_up, dict):
+ self.E_up = E_up
+ else:
+ raise TypeError("E_up must be a dict")
+
+# # provo ad assegnare un magic method in modo che ID_up
+# # sia sempre un set
+# def _get_up(self):
+# return self.up
+#
+# def _set_up(self, up):
+# if isinstance(up, set):
+# self.up = up
+# else:
+# raise TypeError("up must be a set")
+#
+# up = property(fget=_get_up, fset=_set_up)
+
+ def _get_E_up(self):
+ return self._E_up
+
+ def _set_E_up(self, E_up):
+ if isinstance(E_up, dict):
+ self._E_up = E_up
+ else:
+ raise TypeError("E_up must be a dict")
+
+ E_up = property(fget=_get_E_up, fset=_set_E_up)
+
+ def add_E_up(self, ID, E=None):
+ self.E_up[ID] = E
+ return self.E_up[ID]
+
+ def check_ID_up(self, lista):
+ """
+ :param lista: dictionary of all the basins belong to the study area
+ with key the ID of the basin
+ :type h_mean: dictionary of object Basin, key = ID basin, value =
+ object Basin
+ """
+ self.ID_all = self.up
+ #self.ID_all.union(set([99]))
+ # indice = 0
+ if not self.ID_all:
+ return
+ else:
+ for ID in self.up:
+ new = lista[ID].up
+ check_new(self, lista, new)
+
+ def E_spec(self):
+ """
+ compute the specific energy for length unit of the basin
+ """
+ #pdb.set_trace()
+ E_spec = self.E_own/self.area # kW/km2
+ if self.up:
+ for i in self.E_up:
+ E_spec = E_spec + self.E_up[i]/(self.length/1000.0)
+ return E_spec
+
+
+class Station(object):
+
+ def __init__(self, ID, ID_bas=None, days=None, discharge=None,
+ area=None, coord=None):
+ """Instantiate a Basin object.: Define the class Station in order
+ to create a list
+ of object station with duration curve and point coordinates.
+
+ :param ID: specify the ID of the station
+ :type ID: string or int
+ :param ID_bas: ID of the basins
+ :type h_mean: string or int
+ :param days: days of the duration curve [days]
+ :type h_closure: int
+ :param discharge: discharge realted to days [m3/s]
+ :type discharge_own: float
+ :param coord: coordinate of the station point
+ :type discharge_tot: tuple
+ >>> B0 = Station(ID=1, area=100, days=[1, 10, 365],
+ ... discharge=[15, 4, 1])
+ >>> B0.mean()
+ 2.6657534246575341
+
+ @author: GGaregnani
+ """
+ self.ID = ID
+ self.ID_bas = ID_bas
+ self.days = days
+ self.discharge = discharge
+ self.coord = coord
+
+ def mean(self):
+ """Return the mean value of discharge in the duration curve
+ """
+ Q = integrate.trapz(self.discharge, self.days)/365
+ return Q
+
+
+def write_results2newvec(stream, E, basins_tot, inputs):
+ """
+ Create the stream vector and
+ write the basins object in a vector with the same cat value
+ of the ID basin
+ """
+ gcore.run_command("r.thin", input=stream, output='stream_thin')
+ gcore.run_command("r.to.vect", input='stream_thin',
+ flags='v',
+ output='vec_clean', type="line")
+ gcore.run_command("v.edit", map='vec_clean', tool='delete', cats='0')
+ #pdb.set_trace()
+ gcore.run_command('v.build', map='vec_clean')
+ dissolve_lines('vec_clean', E)
+ # TODO: dissolve the areas with the same cat
+ # adding columns
+ gcore.run_command("v.db.addcolumn", map=E,
+ columns="E_spec double precision,"
+ "Qown double precision,"
+ "Qtot double precision, Hmean double precision,"
+ "H0 double precision, Eown_kW double precision,"
+ "IDup1 int, Eup1_kW double precision,"
+ "IDup2 int, Eup2_kW double precision,"
+ "IDup3 int, Eup3_kW double precision")
+ gcore.run_command("db.dropcolumn", flags="f",
+ table=E, column="label")
+ # Open database connection
+ vec = VectorTopo(E)
+ vec.open("rw")
+ link = vec.dblinks[0]
+ conn = link.connection()
+ # prepare a cursor object using cursor() method
+ cursor = conn.cursor()
+ # I modify with specific power (kW/km) see
+ # 4._Julio_Alterach_-_Evaluation_of_the_residual_potential_
+ # hydropower_production_in_Italy
+ # compute the lenght of the river in a basin
+ #import ipdb; ipdb.set_trace()
+ for ID in inputs:
+ length = 0
+ for l in vec.cat(ID, 'lines'):
+ length += l.length()
+ basins_tot[ID].length = length
+ db = [basins_tot[ID].E_spec(),
+ basins_tot[ID].discharge_own,
+ basins_tot[ID].discharge_tot,
+ basins_tot[ID].h_mean,
+ basins_tot[ID].h_closure,
+ basins_tot[ID].E_own]
+ if len(basins_tot[ID].E_up) == 0:
+ db = db + [0, 0.0, 0, 0.0, 0, 0.0]
+ elif len(basins_tot[ID].E_up) == 1:
+ db = (db + [basins_tot[ID].E_up.keys()[0],
+ basins_tot[ID].E_up.values()[0],
+ 0, 0.0, 0, 0.0])
+ elif len(basins_tot[ID].E_up) == 2:
+ db = (db + [basins_tot[ID].E_up.keys()[0],
+ basins_tot[ID].E_up.values()[0],
+ basins_tot[ID].E_up.keys()[1],
+ basins_tot[ID].E_up.values()[1],
+ 0, 0])
+ elif len(basins_tot[ID].E_up) == 3:
+ #pdb.set_trace()
+ db = (db + [basins_tot[ID].E_up.keys()[0],
+ basins_tot[ID].E_up.values()[0],
+ basins_tot[ID].E_up.keys()[1],
+ basins_tot[ID].E_up.values()[1],
+ basins_tot[ID].E_up.keys()[2],
+ basins_tot[ID].E_up.values()[2]])
+ else:
+ db = db + [0, 0.0, 0, 0.0, 0, 0.0]
+ #print db
+ # TODO: change values, give only key and vals without key
+ vec.table.update(basins_tot[ID].ID, db, cursor)
+
+ # disconnect from server
+ conn.commit()
+ conn.close()
+ vec.close()
+
+
+def add_results2shp(basins, rivers, basins_tot, E, inputs):
+ """
+ Add the attribute of the object basins to
+ a vector with the category equal to the ID of the basin
+ """
+ gcore.run_command("r.to.vect", input=basins,
+ flags='v',
+ output='basins', type="area")
+ gcore.run_command('v.overlay', overwrite=True,
+ ainput=rivers, binput=basins, operator='and',
+ output=E)
+ gcore.run_command("v.db.addcolumn", map=E,
+ columns="E_spec double precision")
+ for ID in inputs:
+ #pdb.set_trace()
+ E_spec = str(basins_tot[ID].E_spec())
+ cond = 'b_cat=%i' % (ID)
+ gcore.run_command('v.db.update', map=E,
+ layer=1, column='E_spec',
+ value=E_spec, where=cond)
+
+
+def init_basins(basins):
+ """
+ Create a dictionary of Basin objects with the keys equal
+ to the ID of each basins
+ """
+ # I use r.stats because of the case with MASK
+ info = gcore.parse_command('r.stats', flags='n', input=basins)
+ inputs = map(int, info.keys())
+ basins_tot = {}
+
+ for count in inputs:
+ basins_tot[count] = Basin(ID=count) # discharge_own=Q_basin)
+ return basins_tot, inputs
+
+
+def check_compute_drain_stream(drain, stream, dtm, threshold=100000):
+ """
+ Compute the stream and drain map with r.watersheld
+ """
+ msgr = get_msgr()
+ if (not(drain) or not(stream)):
+ drain = 'drainage'
+ stream = 'stream'
+ gcore.run_command('r.watershed',
+ elevation=dtm,
+ threshold=threshold,
+ drainage=drain,
+ stream=stream)
+ else:
+ info1 = gcore.parse_command('r.info', flags='e', map=drain)
+ info2 = gcore.parse_command('r.info', flags='e', map=stream)
+ #pdb.set_trace()
+ in1 = '\"%s\"' % (dtm)
+ in2 = '\"%s\"' % (dtm)
+ if ((in1 != info1['source1'] or (info1['description']
+ != '\"generated by r.watershed\"'))):
+ warn = ("%s map not generated "
+ "by r.watershed starting from %s") % (drain, dtm)
+ msgr.warning(warn)
+ if ((in2 != info2['source1'] or (info2['description']
+ != '\"generated by r.watershed\"'))):
+ warn = ("%s map not generated "
+ "by r.watershed starting from %s") % (stream, dtm)
+ msgr.warning(warn)
+ return drain, stream
+
+
+def check_compute_basin_stream(basins, stream, dtm, threshold):
+ """
+ Compute the stream and basin map with r.watersheld
+ """
+ msgr = get_msgr()
+ if (not(basins) or not(stream)):
+ basins = 'basins'
+ stream = 'stream'
+ gcore.run_command('r.watershed',
+ elevation=dtm,
+ threshold=threshold,
+ basin=basins,
+ stream=stream)
+ else:
+ info1 = gcore.parse_command('r.info', flags='e', map=basins)
+ info2 = gcore.parse_command('r.info', flags='e', map=stream)
+ #pdb.set_trace()
+ in1 = '\"%s\"' % (dtm)
+ if ((in1 != info1['source1'] or (info1['description']
+ != '\"generated by r.watershed\"'))):
+ warn = ("%s map not generated "
+ "by r.watershed starting from %s") % (basins, dtm)
+ msgr.warning(warn)
+ if ((in1 != info2['source1'] or (info2['description']
+ != '\"generated by r.watershed\"'))):
+ warn = ("%s map not generated "
+ "by r.watershed starting from %s") % (stream, dtm)
+ msgr.warning(warn)
+ return basins, stream
+
+
+def build_network(stream, dtm, basins_tot):
+ """
+ Build the network of streams with the ID of all the basins
+ in order to know the dependencies among basins
+ """
+ river = raster2numpy(stream)
+ river_comp = raster2compressM(stream).tocoo()
+ gcore.run_command('r.neighbors', input=stream,
+ output="neighbors", method="minimu", size='5',
+ quantile='0.5')
+
+ formula = 'closure = neighbors-%s' % (stream)
+ mapcalc(formula)
+ # del the map down, it should be not necessary
+ gcore.run_command('r.stats.zonal',
+ base=stream,
+ cover='closure',
+ output='down',
+ method='min')
+ #pdb.set_trace()
+ dtm_n = raster2numpy(dtm)
+ clos = raster2numpy('closure')
+ ID_down = raster2numpy('down')
+ #pdb.set_trace()
+ for i, j, v in itertools.izip(river_comp.row,
+ river_comp.col, river_comp.data):
+ up = clos[i][j]
+ if up < 0:
+ ID = river[i, j]
+ basins_tot[(ID+int(ID_down[i, j]))].up.add(ID)
+ basins_tot[ID].h_closure = dtm_n[i, j]
+
+
+def area_of_basins(basins, count, dtm):
+ """
+ By knowing the basin map compute the area of a given basin ID=count
+ """
+ #TODO: check the similar function in r.green.discharge
+ command = 'area_stat=if(%s == %i, %s, null())' % ((basins, count, dtm))
+ mapcalc(command, overwrite=True)
+ return area_of_watershed('area_stat')
+
+
+def area_of_watershed(watershed):
+ info = gcore.parse_command('r.report', map=watershed, flags='hn',
+ units='k')
+ temp = 0
+ for somestring in info.keys():
+ if ('TOTAL') in somestring:
+ temp = float(somestring.split('|')[-2])
+ return temp
+
+
+def fill_energyown(bas, h_mean, discharge_n, stream_n):
+ """
+ Fill basins dictionary with discharge, h_mean and compute the power
+ """
+ msgr = get_msgr()
+ ttt = discharge_n[stream_n == bas.ID]
+ #import ipdb; ipdb.set_trace()
+ bas.discharge_own = ttt.sort()
+ #FIXME: take the second bgger value to avoid to take the value of
+ # another catchment, it is not so elegant
+ if len(ttt) > 1:
+ bas.discharge_own = float(ttt[-2])
+ elif ((len(ttt) < 1) and (len(ttt) > 0)):
+ bas.discharge_own = float(ttt[0])
+ warn = ("Only one value of discharge for the river ID %i") % bas.ID
+ msgr.warning(warn)
+ else:
+ warn = ("Only one value of discharge for the river ID %i") % bas.ID
+ msgr.warning(warn)
+ bas.h_mean = h_mean
+ # basins_tot[count].h_closure = float(info_c[count])
+ delta = bas.h_mean - bas.h_closure
+ #TODO: modify power with potential
+ bas.E_own = (E_hydro(delta, bas.discharge_own))
+ #pdb.set_trace()
+
+
+def fill_discharge_tot(basins_tot, b):
+ """
+ Fill the total discharge for the basin b by knowing
+ the relation among basins
+ thank to the ID_all attribute in the object Basin
+ """
+ basins_tot[b].discharge_tot = basins_tot[b].discharge_own
+ #import pdb; pdb.set_trace()
+ if basins_tot[b].up:
+ basins_tot[b].check_ID_up(basins_tot)
+ temp = discharge_sum(basins_tot, basins_tot[b].ID_all)
+ basins_tot[b].discharge_tot = basins_tot[b].discharge_own + temp
+
+
+def fill_Eup(basins_tot, b):
+ """
+ Fill the energy coming from discharge in the upper basins
+ """
+ for ID in basins_tot[b].up:
+ deltaH = basins_tot[ID].h_closure - basins_tot[b].h_closure
+ power = E_hydro(deltaH, basins_tot[ID].discharge_tot)
+ # if area > 1km2 I use the up-up basin
+ #TODO: change the up basins for the energy, something
+ # similar to an optimization problem
+ if basins_tot[ID].area > 1:
+ basins_tot[b].add_E_up(ID, power)
+ else:
+ if basins_tot[ID].up:
+ for id_up in basins_tot[ID].up:
+ deltaH = (basins_tot[id_up].h_closure -
+ basins_tot[b].h_closure)
+ power = E_hydro(deltaH, basins_tot[id_up].discharge_tot)
+ basins_tot[b].add_E_up(id_up, power)
+ else:
+ basins_tot[b].add_E_up(ID, power)
+
+
+def fill_basins(inputs, basins_tot, basins, dtm, discharge, stream):
+ """
+ Fill the dictionary with the basins
+ """
+ info_h = gcore.parse_command('r.category', map='dtm_mean', separator='=')
+ #pdb.set_trace()
+ for count in inputs:
+ if info_h[str(count)] != '':
+ area = area_of_basins(basins, count, dtm)
+ basins_tot[count].area = float(area)
+ h_mean = float(info_h[str(count)])
+ fill_energyown(basins_tot[count], h_mean, discharge, stream)
+
+ for b in inputs:
+ fill_discharge_tot(basins_tot, b)
+
+ for b in inputs:
+ fill_Eup(basins_tot, b)
+
+
+def compute_river_discharge(drain, stream, string, **kwargs):
+ """
+ Given a stream network and drainage map
+ it computes a rester with the
+ the given resolution
+ with the area in km2 of the upper basin
+ for each pixel (bas_area)
+ and a statistic on the bas_area for another series of raster
+ if q_spec string=sum
+ if piedmont case kwargs-> a=a, dtm=dtm and string=mean
+ """
+ msgr = get_msgr()
+ info = gcore.parse_command('r.info', flags='g', map=stream)
+ raster_out = {}
+ for name, value in kwargs.items():
+ raster_out[name] = raster2numpy(stream)
+
+ bas_area = raster2numpy(stream)
+ river_comp = raster2compressM(stream).tocoo()
+ count = 0
+ for i, j in itertools.izip(river_comp.row, river_comp.col):
+ count = count+1
+ msgr.message("\n %i \n" % count)
+ p_x, p_y = get_coo(stream, i, j)
+ #pdb.set_trace()
+ coo = '%f, %f' % (p_x, p_y)
+ gcore.run_command('r.water.outlet', overwrite=True, input=drain,
+ output='area_temp',
+ coordinates=coo)
+ for name, value in kwargs.items():
+ formula = 'temp = if(not(area_temp),0, %s)' % (value)
+ #pdb.set_trace()
+ mapcalc(formula, overwrite=True)
+ #pdb.set_trace()
+ info = gcore.parse_command('r.univar', map='temp', flags='g')
+ #pdb.set_trace()
+ raster_out[name][i][j] = float(info[string])
+ bas_area[i][j] = area_of_watershed('area_temp')
+ #pdb.set_trace()
+ return raster_out, bas_area
+
+
+def dtm_corr(dtm, river, lake=None):
+ temp_vec = []
+ temp_rast = []
+ msgr = get_msgr()
+ info = gcore.parse_command('g.region', flags='pg')
+ """ change dtm accordingo to the river and lake vectors"""
+ if lake:
+ inputs = '%s,%s' % (lake, river)
+ gcore.run_command('v.patch',
+ input=inputs,
+ output='hydro_network')
+ temp_vec.append('hydro_network')
+ river = 'hydro_network'
+
+ msgr.warning("The DTM will be temporarily modified")
+ distance = [float(info['nsres']), float(info['nsres'])*1.5,
+ float(info['nsres'])*3]
+ for i, val in enumerate(distance):
+ output = 'buff_%i' % i
+ gcore.run_command('v.buffer',
+ input=river,
+ output=output,
+ distance=val)
+ temp_vec.append(output)
+ gcore.run_command('v.to.rast',
+ input=output,
+ output=output,
+ use='val',
+ value=val,
+ overwrite=True)
+ temp_rast.append(output)
+ command = ('%s_c = if(isnull(%s),0,%s)') % (output, output,
+ output)
+ mapcalc(command, overwrite=True)
+ temp_rast.append('%s_c' % output)
+
+ command = (('DTM_corr = if(%s,%s-buff_0_c-buff_1_c-buff_2_c)')
+ % (dtm, dtm))
+ mapcalc(command, overwrite=True)
+ return 'DTM_corr', temp_vec, temp_rast
+
+
+if __name__ == "__main__":
+ import doctest
+ doctest.testmod()
Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/libhydro/basin.py
___________________________________________________________________
Added: svn:mime-type
+ text/x-python
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.green/r.green.hydro/libhydro/optimal.py
===================================================================
--- grass-addons/grass7/raster/r.green/r.green.hydro/libhydro/optimal.py (rev 0)
+++ grass-addons/grass7/raster/r.green/r.green.hydro/libhydro/optimal.py 2015-04-27 06:36:52 UTC (rev 65146)
@@ -0,0 +1,349 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Jan 12 14:29:46 2015
+
+ at author: ggaregnani
+"""
+# import system libraries
+from __future__ import print_function
+import os
+import sys
+from scipy import interpolate
+from scipy.optimize import fsolve
+import numpy as np
+
+# import scientific libraries
+from scipy import optimize
+
+#from grass.script import mapcalc
+from grass.script import core as gcore
+from grass.pygrass.messages import get_msgr
+from grass.pygrass.raster import RasterRow
+
+#from grass.pygrass.raster.buffer import Buffer
+from grass.pygrass.gis.region import Region
+from grass.pygrass.vector import VectorTopo
+
+#from grass.pygrass.raster.buffer import Buffer
+from grass.pygrass.utils import set_path
+
+set_path('r.green', 'libhydro', '..')
+set_path('r.green', 'libgreen', os.path.join('..', '..'))
+
+from libhydro.plant import Plant
+from libhydro.plant import Intake
+from libhydro.plant import Restitution
+from libhydro.plant import COLS
+from libhydro.plant import COLS_points
+
+
+def f(x, *params):
+ #msgr = get_msgr()
+ prog, h, q = params
+ #msgr.message("\%f, %f\n" % (x[0], x[1]))
+ #msgr.message("\init int\n")
+ fun_h = interpolate.interp1d(prog, h)
+ fun_q = interpolate.interp1d(prog, q)
+ #msgr.message("\end int\n")
+ dh = fun_h(x[0]) - fun_h(x[0]+x[1])
+ # pdb.set_trace()
+ return -dh*fun_q(x[0])
+
+
+def find_s(s, *params):
+ #msgr = get_msgr()
+ x0, prog, h, q, p_max = params
+ fun_h = interpolate.interp1d(prog, h, bounds_error=False, fill_value=0)
+ fun_q = interpolate.interp1d(prog, q, bounds_error=False, fill_value=0)
+ if fun_q(x0) == 0:
+ #import pdb; pdb.set_trace()
+ return 99999
+ else:
+ return p_max-(fun_q(x0)*9.81)*(fun_h(x0) - fun_h(x0+s))
+
+
+def min_len_plant(prog, h, q, p_max, range_plant,
+ (start, end)):
+ q = np.array(q)
+ q[np.isnan(q)] = 0
+# if q.sum() == 0:
+# return (0, 0)
+ len_plant = []
+ fun_q = interpolate.interp1d(prog, q, bounds_error=False, fill_value=0)
+ if end - start > 10:
+ num = (end-start-10)/10
+ for x0 in np.linspace(start, end-1, num=num):
+ #import pdb; pdb.set_trace()
+ if fun_q(x0) == 0:
+ #import pdb; pdb.set_trace()
+ len_plant.append(0)
+ else:
+ #import pdb; pdb.set_trace()
+ sol, infodict, ier, mesg = fsolve(find_s, [0], args=(x0,
+ prog, h, q, p_max),
+ full_output=True)
+ if ier == 1:
+ #import pdb; pdb.set_trace()
+ len_plant.append(sol[0])
+ else:
+ len_plant.append(0)
+
+ len_plant = np.array(len_plant)
+ order = np.argsort(len_plant)
+ x_ink = np.linspace(start, end-10, num=num)[order]
+ x_rest = x_ink + len_plant[order]
+ s = len_plant[order]
+ aaa = x_rest < end
+ bbb = s > 0
+ ccc = aaa * bbb
+ if ccc.any():
+ if len_plant.sum() == 0:
+ return (start, end-start)
+ else:
+ return (x_ink[ccc][0],
+ s[ccc][0])
+ else:
+ return (start, end-start)
+ else:
+ return (start, end-start)
+
+
+#TODO: to move in another lib
+def check_multilines(vector):
+ vector, mset = vector.split('@') if '@' in vector else (vector, '')
+ msgr = get_msgr()
+ vec = VectorTopo(vector, mapset=mset, mode='r')
+ vec.open("r")
+ info = gcore.parse_command('v.category',
+ input=vector, option='print')
+ for i in info.keys():
+ vec.cat(int(i), 'lines', 1)
+# if i == '28':
+# import ipdb; ipdb.set_trace()
+ if len(vec.cat(int(i), 'lines', 1)) > 1:
+ # import ipdb; ipdb.set_trace()
+ warn = ("Multilines for the same category %s" % i)
+ msgr.warning(warn)
+ vec.close()
+
+
+def build_array(line, raster_q, raster_dtm):
+ """
+ build arrays with discharge, elevation and
+ coordinate + order the line according to dtm and discharge
+ """
+ prog = []
+ h = []
+ q = []
+ coordinate = []
+
+ for point in line:
+ coordinate.append(point)
+ q.append(raster_q.get_value(point))
+ h.append(raster_dtm.get_value(point))
+ #FIXME: the second and the second to last because
+ # we should avoid to take discharge in another branch
+ #import ipdb; ipdb.set_trace()
+ if len(h) > 3:
+ h_diff = np.array(h[0: -2])-np.array(h[1: -1])
+ if h_diff.sum() < 0:
+ q = q[::-1]
+ h = h[::-1]
+ coordinate = coordinate[::-1]
+ line.reverse()
+
+ prog.append(0)
+ # note that i is i-1
+ for i, point in enumerate(coordinate[1::]):
+ dist = point.distance(coordinate[i])
+ #pdb.set_trace()
+ prog.append(prog[i]+dist)
+ return line, prog, h, q
+
+
+def find_optimal(args, range_plant, ranges):
+
+ eps = 0.5 # for the grid
+ rranges = ((ranges[0], ranges[1] - range_plant[1] - eps),
+ range_plant)
+ x_min = optimize.brute(f, rranges, args=args, finish=None)
+ return x_min
+
+
+def recursive_plant(args, range_plant, distance,
+ start, end, rank, cat, line, plants, count, p_max):
+ msgr = get_msgr()
+ count = count + 1
+ #import ipdb; ipdb.set_trace()
+ msgr.message("\n%i\n" % cat)
+ #import pdb; pdb.set_trace()
+ res = check_plant(args, range_plant, distance, start, end, rank,
+ cat, line, count, p_max)
+ if res:
+ plants.append(res[0])
+ plants = check_segments(args, range_plant, distance, res[1], start,
+ end, rank, cat, line, plants, count, p_max)
+ return plants
+
+
+def check_plant(args, range_plant, distance, start, end, rank, cat,
+ line, count, p_max):
+ #import pdb; pdb.set_trace()
+ prog, h, q = args
+ fun_h = interpolate.interp1d(prog, h, bounds_error=False, fill_value=0)
+ fun_q = interpolate.interp1d(prog, q, bounds_error=False, fill_value=0)
+ len_plant = range_plant[1]
+ len_min = range_plant[0]
+ if count < 6:
+ if (end-start) > len_plant + 2*distance:
+ if not(p_max):
+ x_min = find_optimal(args, range_plant,
+ (start+distance, end-distance))
+ else:
+ x_min = min_len_plant(prog, h, q, p_max, range_plant,
+ (start+distance, end-distance))
+ seg = line.segment(x_min[0], x_min[0]+x_min[1])
+ if seg:
+ ink = Intake(id_point=str(1), id_plants=rank, id_stream=cat,
+ point=x_min[0], discharge=fun_q(x_min[0]),
+ elevation=fun_h(x_min[0]))
+ res = Restitution(id_point=str(2), id_plants=rank,
+ id_stream=cat,
+ point=x_min[0]+x_min[1],
+ discharge=fun_q(x_min[0]),
+ elevation=fun_h(x_min[0]+x_min[1]))
+ p = Plant(id_plant=rank, id_stream=cat,
+ restitution=res, intakes=(ink,),
+ line=seg)
+ return p, x_min
+ elif ((end-start-2*distance) > len_min
+ and ((end-start) < len_plant + 2*distance)):
+ len_plant = end-start-2*distance
+ range_plant = (len_min, len_plant)
+ if not(p_max):
+ x_min = find_optimal(args, range_plant,
+ (start+distance, end-distance))
+ else:
+ x_min = min_len_plant(prog, h, q, p_max, range_plant,
+ (start+distance, end-distance))
+ seg = line.segment(x_min[0], x_min[0]+x_min[1])
+ if seg:
+ ink = Intake(id_point=str(1), id_plants=rank, id_stream=cat,
+ point=x_min[0], discharge=fun_q(x_min[0]),
+ elevation=fun_h(x_min[0]))
+ #import pdb; pdb.set_trace()
+ res = Restitution(id_point=str(2), id_plants=rank,
+ id_stream=cat,
+ point=x_min[0]+x_min[1],
+ discharge=fun_q(x_min[0]),
+ elevation=fun_h(x_min[0]+x_min[1]))
+ p = Plant(id_plant=rank, id_stream=cat,
+ restitution=res, intakes=(ink,),
+ line=seg)
+ return p, x_min
+ #FIXME: check the minimum segment for hydroplants,
+ # perhaps it is better a minimum gross head
+ #return
+ # pdb.set_trace()
+# seg1 = line.segment(start, x_min[0])
+# check_plant(seg1, args, len_plant, new_vec, start, x_min[0],
+# rank+1)
+# seg2 = line.segment(x_min[0]+x_min[1], end)
+# check_plant(seg2, args, len_plant, new_vec, start, x_min[0],
+# rank+1)
+
+
+def check_segments(args, range_plant, distance, x_min, start_old,
+ end_old, rank, cat, line, plants, count, p_max):
+ #seg1 = line.segment(start_old, x_min[0])
+ rank = rank+'.1'
+ plants = recursive_plant(args, range_plant, distance, start_old, x_min[0],
+ rank, cat, line, plants, count, p_max)
+ #seg2 = line.segment(x_min[0]+x_min[1], end_old)
+ rank = rank+'.2'
+ plants = recursive_plant(args, range_plant, distance, x_min[0]+x_min[1],
+ end_old, rank, cat, line, plants, count, p_max)
+ return plants
+
+
+# else:
+# return
+#def find_residual(line, prog, h, q, x_min):
+# rest = []
+# rest[1] = line.segment(prog)
+
+
+def find_segments(river, discharge, dtm, range_plant, distance, p_max):
+ check_multilines(river)
+ #pdb.set_trace()
+ river, mset = river.split('@') if '@' in river else (river, '')
+ vec = VectorTopo(river, mapset=mset, mode='r')
+ vec.open("r")
+ raster_q = RasterRow(discharge)
+ raster_dtm = RasterRow(dtm)
+ raster_q.open('r')
+ raster_dtm.open('r')
+ reg = Region()
+ plants = []
+ for line in vec:
+ count = 0
+ # args is prog, h, q
+ line, prog, h, q = build_array(line, raster_q, raster_dtm)
+ #pdb.set_trace()
+ if len(line) > 2:
+# import ipdb; ipdb.set_trace()
+# else:
+ # import pdb; pdb.set_trace()
+ plants = recursive_plant((prog, h, q), range_plant, distance,
+ prog[0], prog[-1],
+ str(line.cat), line.cat, line, plants,
+ count, p_max)
+ #pdb.set_trace()
+ vec.close()
+ raster_q.close()
+ raster_dtm.close()
+ return plants
+
+
+def write_plants(plants, output, efficiency):
+ # create vetor segment
+ new_vec = VectorTopo(output)
+ #TODO: check if the vector already exists
+ new_vec.layer = 1
+ new_vec.open('w', tab_cols=COLS)
+ reg = Region()
+ for pla in plants:
+ new_vec.write(pla.line, (pla.id, pla.id_stream,
+ pla.potential_power(efficiency=efficiency)))
+
+ new_vec.table.conn.commit()
+ new_vec.comment = (' '.join(sys.argv))
+ #pdb.set_trace()
+ new_vec.close()
+
+
+def write_points(plants, output, efficiency):
+ # create vetor segment
+ new_vec = VectorTopo(output)
+ #TODO: check if the vector already exists
+ new_vec.layer = 1
+ new_vec.open('w', tab_cols=COLS_points)
+ reg = Region()
+ # import ipdb; ipdb.set_trace()
+ for pla in plants:
+ new_vec.write(pla.line[-1], (pla.restitution.id,
+ pla.id, 'restitution', pla.id_stream,
+ float(pla.restitution.elevation),
+ float(pla.restitution.discharge),
+ pla.potential_power()))
+ for ink in pla.intakes:
+ new_vec.write(pla.line[0], (ink.id,
+ pla.id, 'intake', pla.id_stream,
+ float(ink.elevation), float(ink.discharge),
+ pla.potential_power(efficiency=efficiency)))
+
+ new_vec.table.conn.commit()
+ new_vec.comment = (' '.join(sys.argv))
+ new_vec.write_header()
+ #pdb.set_trace()
+ new_vec.close()
Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/libhydro/optimal.py
___________________________________________________________________
Added: svn:mime-type
+ text/x-python
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.green/r.green.hydro/libhydro/plant.py
===================================================================
--- grass-addons/grass7/raster/r.green/r.green.hydro/libhydro/plant.py (rev 0)
+++ grass-addons/grass7/raster/r.green/r.green.hydro/libhydro/plant.py 2015-04-27 06:36:52 UTC (rev 65146)
@@ -0,0 +1,469 @@
+# -*- coding: utf-8 -*-
+import random
+from collections import namedtuple
+
+import numpy as np
+
+from grass.pygrass.modules.shortcuts import raster as r
+from grass.pygrass.gis.region import Region
+from grass.pygrass.vector import VectorTopo
+from grass.pygrass.vector.geometry import Line
+from grass.pygrass.vector.table import Link
+
+
+COLS = [(u'cat', 'INTEGER PRIMARY KEY'),
+ (u'plant_id', 'VARCHAR(10)'),
+ (u'stream_id', 'INTEGER'),
+ (u'potential_power', 'DOUBLE')]
+
+
+COLS_points = [(u'cat', 'INTEGER PRIMARY KEY'),
+ (u'kind', 'VARCHAR(10)'),
+ (u'plant_id', 'VARCHAR(10)'),
+ (u'kind_label', 'VARCHAR(10)'),
+ (u'stream_id', 'INTEGER'),
+ (u'elevation', 'DOUBLE'),
+ (u'discharge', 'DOUBLE'),
+ (u'potential_power', 'DOUBLE')]
+
+HydroStruct = namedtuple('HydroStruct',
+ ['intake', 'conduct', 'penstock', 'side'])
+
+
+def isinverted(line, elev, region):
+ el0, el1 = elev.get_value(line[0]), elev.get_value(line[-1])
+ return False if el0 > el1 else True
+
+
+def not_overlaped(line):
+ """The countur lines are always a ring even if we see tham as line
+ they are just overlaped, we want to avoid this situation
+ therefore we return only the part of the line that is not overlaped
+ """
+ if len(line) >= 2 and line[1] == line[-2]:
+ return Line(line[:len(line)//2+1])
+ return line
+
+
+def sort_by_elev(line, elev, region):
+ """Return a Line object sorted using the elevation, the first point will
+ be the higher and the last the lower.
+ """
+ if isinverted(line, elev, region):
+ line.reverse()
+ return line
+
+
+def sort_by_west2east(line):
+ """Return a Line object sorted using east coordinate, the first point will
+ be the wester and the last the easter.
+ """
+ if line[0].x > line[-1].x:
+ line.reverse()
+ return line
+
+
+def splitline(line, point):
+ """Split a line using a point. return two lines that start with the point.
+
+ point
+ *---------
+ \ \
+ \ -------->
+ \
+ --------|
+
+ """
+ dist = line.distance(point)
+ l0 = line.segment(0, dist.sldist)
+ l0.reverse()
+ l1 = line.segment(dist.sldist, line.length())
+ return l0, l1
+
+
+
+def read_plants(hydro, elev=None, restitution='restitution', intake='intake',
+ cid_plant='id_plant', cid_point='id_point', ckind_label='kind_label',
+ celevation='elevation', cdischarge='discharge'):
+ plants = {}
+ skipped = []
+ for pnt in hydro:
+ if pnt is None:
+ import ipdb; ipdb.set_trace()
+ if elev is None:
+ select = ','.join([cid_plant, cid_point, ckind_label, celevation, cdischarge])
+ id_plant, id_point, kind_label, el, disch = pnt.attrs[select]
+ else:
+ select = ','.join([cid_plant, cid_point, ckind_label, cdischarge])
+ id_plant, id_point, kind_label, disch = pnt.attrs[select]
+ el = elev.get_value(pnt)
+ plant = plants.get(id_plant, Plant(id_plant))
+ if kind_label == restitution:
+ plant.restitution = Restitution(id_point, pnt, el)
+ elif kind_label == intake:
+ plant.intakes.add(Intake(id_point=id_point, point=pnt,
+ elevation=el, id_plants=id_plant,
+ discharge=disch if disch else 1.))
+ else:
+ skipped.append((id_plant, id_point, kind_label))
+ plants[id_plant] = plant
+ return plants, skipped
+
+
+def write_plants(plants, output, stream, elev, overwrite=False):
+ """Write a vector map with the plant"""
+ with VectorTopo(output, mode='w', tab_cols=COLS,
+ overwrite=overwrite) as out:
+ for p in plants:
+ potential_power= plants[p].potential_power()
+ plant_id = plants[p].id
+ lines, ids = plants[p].plant(stream, elev)
+ for line, r_id in zip(lines, ids):
+ out.write(line, (plant_id, r_id, potential_power, ))
+ out.table.conn.commit()
+
+
+def write_structures(plants, output, elev, stream=None, overwrite=False):
+ """Write a vector map with the plant structures"""
+ def write_hydrostruct(out, hydro, plant):
+ #import pdb; pdb.set_trace()
+ pot = plant.potential_power(intakes=[hydro.intake, ])
+ (plant_id, itk_id, side,
+ disch, gross_head) = (plant.id, hydro.intake.id, hydro.side,
+ float(hydro.intake.discharge),
+ float(hydro.intake.elevation -
+ plant.restitution.elevation))
+ out.write(hydro.conduct,
+ (plant_id, itk_id, disch, 0., 0., 'conduct', side))
+ out.write(hydro.penstock,
+ (plant_id, itk_id, disch, gross_head, pot, 'penstock', side))
+ out.table.conn.commit()
+
+ tab_cols = [(u'cat', 'INTEGER PRIMARY KEY'),
+ (u'plant_id', 'VARCHAR(10)'),
+ (u'intake_id', 'INTEGER'),
+ (u'discharge', 'DOUBLE'),
+ (u'gross_head', 'DOUBLE'),
+ (u'power', 'DOUBLE'),
+ (u'kind', 'VARCHAR(10)'),
+ (u'side', 'VARCHAR(10)'), ]
+ with VectorTopo(output, mode='w', overwrite=overwrite) as out:
+ link = Link(layer=1, name=output, table=output,
+ driver='sqlite')
+ out.open('w')
+ out.dblinks.add(link)
+ out.table = out.dblinks[0].table()
+ out.table.create(tab_cols)
+ print('Number of plants: %d' % len(plants))
+ for p in plants:
+ plant = plants[p]
+ print(plant.id)
+ #import ipdb; ipdb.set_trace
+ for options in plant.structures(elev, stream=stream):
+ for hydro in options:
+ print('writing: ', hydro.intake)
+ write_hydrostruct(out, hydro, plant)
+
+
+class AbstractPoint(object):
+ def __init__(self, id_point, point, elevation, id_plants=None,
+ id_stream=None, discharge=1.):
+ self.id = id_point
+ self.point = point
+ self.elevation = elevation
+ self.id_plants = id_plants
+ self.id_stream = id_stream
+ self.discharge = discharge
+
+ def __repr__(self):
+ srepr = "%s(id_point=%r, point=%r, elevation=%r)"
+ return srepr % (self.__class__.__name__, self.id,
+ self.point, self.elevation)
+
+
+class Restitution(AbstractPoint):
+ pass
+
+
+class Intake(AbstractPoint):
+ pass
+
+
+class Plant(object):
+ def __init__(self, id_plant, id_stream=None, restitution=None,
+ intakes=None, line=None):
+ self.id = id_plant
+ self.id_stream = id_stream
+ self.restitution = restitution
+ self.intakes = set() if intakes is None else set(intakes)
+ self.line = line
+
+ def __repr__(self):
+ srepr = "%s(id_plant=%r, restitution=%r, intakes=%r)"
+ return srepr % (self.__class__.__name__, self.id,
+ self.restitution, self.intakes)
+
+ def potential_power(self, efficiency=1., intakes=None):
+ """Return the potential of the plant: input discharge [m3/s],
+ elevetion [m] and output [kW].
+
+ Parameters
+ ----------
+
+ efficiency: float
+ It is a float value between 0 and 1 to take into account
+ the efficiency of the turbine, default value is 1.
+ intakes: iterable of Intake instances
+ It is an iterable object to specify a subset of the intakes
+ that will be used to compute the potential.
+
+ Returns
+ -------
+
+ potential: float
+ It is the potential of the plant without the flow.
+ """
+ intakes = self.intakes if intakes is None else intakes
+ # create an array with all the elevation of the intakes points
+ # and the discharge
+ elev = np.array([ink.elevation for ink in intakes])
+ discharge = np.array([ink.discharge for ink in intakes])
+ return ((elev - self.restitution.elevation)
+ * discharge * 9.810 * efficiency).sum()
+
+ def plant(self, stream, elev, maxdist=1.0, region=None):
+ """Return a list with the segments involved by the plant.
+
+ Parameters
+ ----------
+
+ strem: vector map
+ A vector map instance that is already opened containing lines with
+ the streams.
+ maxdist: float
+ It is the maximum distance that will be used to find the closest
+ point, if no line is found for the point a TypeError is raised.
+ elev: raster map
+ A rester map instance that it is already opened with the digital
+ elevation values of the study area, if not given the function
+ will perform a network serach to find the path that link the
+ intake and restoration points.
+
+ Returns
+ -------
+
+ lines: a list of lines
+ A dictionary of vector features ids and tuple with the segments. ::
+
+ {id1: ([Line(...), ], [Line(...), ]), # plant intake
+ id2: (Line(...), []), # whole part of the plant
+ ...
+ idN: (Line(...), [Line(...), ]), } # plant restitution
+
+ An example could be an intake (x) and restitution (0) that are
+ on the same line. ::
+
+ <---------------------id--------------------->
+ |============x=================o===============|
+ |<---iseg--->|
+ |<------------rseg------------>|
+ |-----a------|--------b--------|--------c------|
+
+ the list returned is: [b, ]
+
+ For an intake (x) and restitution (o) that are on different lines
+
+ <------id1-----> <------id2----> <---id3---->
+ |======x=========||===============||========o===|
+ |<-i0->|
+ |<--r0-->|
+ |---a--|-----b---|--------c--------|----d---|-e-|
+
+ the list returned is: [b, c, d]
+
+ """
+ def error(pnt):
+ raise TypeError("Line not found for %r" % pnt)
+
+ def split_line(pnt):
+ """Return a tuple with three elements:
+ 0. the line
+ 1. the distance from the beginning of the line
+ 2. the total length of the line.
+ """
+ line = stream.find['by_point'].geo(pnt, maxdist=maxdist)
+ if line is None:
+ error(pnt)
+ (newpnt, _, _, seg) = line.distance(pnt)
+ return line, seg, line.length()
+
+ def find_path(line, res):
+ def lower(line):
+ """Return the lower node of the line"""
+ nfirst, nlast = line.nodes()
+ first = elev.get_value(nfirst, region=region)
+ last = elev.get_value(nlast, region=region)
+ return nlast if first > last else nfirst
+
+ def next_line(node, prev, resid):
+ """Return a generator with the lines until the restitution."""
+ for line in node.lines():
+ if line.id != prev.id:
+ nd = lower(line)
+ if nd.id != node.id:
+ if line.id != resid:
+ yield line
+ for l in next_line(nd, line, resid):
+ yield l
+ else:
+ yield line
+
+ lines = []
+ node = lower(line)
+ lines = [l for l in next_line(node, line, res.id)]
+ return lines
+
+ lines = []
+ ids = []
+ reg = Region() if region is None else region
+ #stream, mset = stream.split('q') if '@' in stream else (stream, '')
+ #with VectorTopo(stream, mapset=mset, mode='r') as stm:
+ res, rseg, rlen = split_line(self.restitution.point)
+ for intake in self.intakes:
+ itk, iseg, ilen = split_line(intake.point)
+ if itk.id == res.id:
+ start, end = (rseg, iseg) if rseg <= iseg else (iseg, rseg)
+ line = res.segment(start, end)
+ lines.append(sort_by_elev(line, elev, reg))
+ ids.append(itk.id)
+ else:
+ if isinverted(itk, elev, reg):
+ iseg = ilen - iseg
+ itk.reverse()
+ line = itk.segment(iseg, ilen)
+ lines.append(line)
+ ids.append(itk.id)
+ for line in find_path(itk, res):
+ l_id = line.id
+ if l_id == res.id:
+ if isinverted(line, elev, reg):
+ rseg = rlen - rseg
+ line.reverse()
+ line = line.segment(0, rseg)
+ lines.append(sort_by_elev(line, elev, reg))
+ ids.append(l_id)
+ return lines, ids
+
+ def structures(self, elev, stream=None):
+ """Return a tuple with lines structres options of a hypotetical plant.
+
+ ::
+
+ river
+ \ \
+ \i\-------_______
+ |\ \ \
+ | \ \ \
+ ) \ \ ) cond0
+ / \ \ /
+ / \ \ /
+ ( cond1 \ \ /
+ \ \ \ |
+ \ \ \ |
+ \ \ \ |
+ o--------\r\---o
+ pstk1 \ \ pstk0
+ \ \
+
+ Parameters
+ ----------
+
+ elev: raster
+ Raster instance already opened with the elevation.
+ intake_pnt: point
+ It is the point of the intake.
+ restitution_pnt: point
+ It is the point of the restitution.
+
+ Returns
+ -------
+
+ a list of tuple: [(HydroStruct(intake, conduct=cond0, penstock=pstk0),
+ HydroStruct(intake, conduct=cond1, penstock=pstk1))]
+ Return a list of tuples, containing two HydroStruct the first with
+ the shortest penstock and the second with the other option.
+ """
+ def get_struct(contur, respoint):
+ """Return the lines of the conduct and the penstock.
+
+ Parameters
+ ----------
+
+ contur: line segment
+ It is a line segment of the contur line splited with splitline
+ function, where the first point is the intake.
+ respoint: point
+ It is the point of the plant restitution.
+
+ Returns
+ -------
+
+ tuple: (conduct, penstock)
+ Return two lines, the first with the conduct and the second with
+ the penstock. Note: both conduct and penstock lines are coherent
+ with water flow direction.
+ """
+ dist = contur.distance(respoint)
+ conduct = contur.segment(0, dist.sldist)
+ penstock = Line([dist.point, respoint])
+ return conduct, penstock
+
+ def get_all_structs(contur, itk, res):
+ l0, l1 = splitline(contur, itk.point)
+ # get structs
+ c0, p0 = get_struct(l0, res.point)
+ c1, p1 = get_struct(l1, res.point)
+ s0, s1 = 'option0', 'option1'
+ # TODO: uncomment this to have left and right distinction...
+ # but sofar is not working properly, therefore is commented.
+ #if stream is not None:
+ # sitk = stream.find['by_point'].geo(itk.point, maxdist=100000)
+ # s0, s1 = (('right', 'left') if isinverted(sitk, elev, reg)
+ # else ('left', 'right'))
+ return (HydroStruct(itk, c0, p0, s0), HydroStruct(itk, c1, p1, s1))
+
+ result = []
+ reg = Region()
+ for itk in self.intakes:
+ tmpvect = 'tmpvect%d' % random.randint(1000, 9999)
+ # generate the contur line that pass to the point
+ # TODO: maybe we can compute the contour with all the levels taked
+ # from the all te intake, just to compute this once
+ r.contour(input='%s@%s' % (elev.name, elev.mapset),
+ output=tmpvect, step=0, levels=itk.elevation,
+ overwrite=True)
+ with VectorTopo(tmpvect, mode='r') as cnt:
+ # find the closest contur line
+ contur_res = cnt.find['by_point'].geo(self.restitution.point,
+ maxdist=100000.0)
+ contur = not_overlaped(contur_res)
+ contur = sort_by_west2east(contur)
+ # TODO: probably find the contur line for the intake and
+ # the restitution it is not necessary, and we could also remove
+ # the check bellow: contur_itk.id != contur_res.id
+ contur_itk = cnt.find['by_point'].geo(itk.point,
+ maxdist=100000.0)
+ if contur_itk is None or contur_res is None:
+ msg = ('Not able to find the contur line closest to the '
+ 'intake point %r, of the plant %r'
+ 'from the contur line map: %s')
+ raise TypeError(msg % (itk, self, cnt.name))
+ if contur_itk.id != contur_res.id:
+ print('=' * 30)
+ print(itk)
+ msg = "Contur lines are different! %d != %d, in %s"
+ print(msg % (contur_itk.id, contur_res.id, cnt.name))
+ result.append(get_all_structs(contur, itk, self.restitution))
+ # remove temporary vector map
+ cnt.remove()
+ return result
Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/libhydro/plant.py
___________________________________________________________________
Added: svn:mime-type
+ text/x-python
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.closest/Makefile
===================================================================
--- grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.closest/Makefile (rev 0)
+++ grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.closest/Makefile 2015-04-27 06:36:52 UTC (rev 65146)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../../../..
+
+PGM = r.green.hydro.closest
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.closest/Makefile
___________________________________________________________________
Added: svn:mime-type
+ text/x-makefile
Added: svn:eol-style
+ native
Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.closest/r.green.hydro.closest.html
___________________________________________________________________
Added: svn:mime-type
+ text/html
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.closest/r.green.hydro.closest.py
===================================================================
--- grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.closest/r.green.hydro.closest.py (rev 0)
+++ grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.closest/r.green.hydro.closest.py 2015-04-27 06:36:52 UTC (rev 65146)
@@ -0,0 +1,96 @@
+#!/usr/bin/env python
+# -- coding: utf-8 --
+#
+############################################################################
+#
+# MODULE: r.green.hydro.potential
+# AUTHOR(S): Pietro Zambelli
+# PURPOSE: Calculate the hydropower energy potential for each basin
+# COPYRIGHT: (C) 2014 by the GRASS Development Team
+#
+# This program is free software under the GNU General Public
+# License (>=v2). Read the file COPYING that comes with GRASS
+# for details.
+#
+#############################################################################
+#
+
+#%module
+#% description: Move points to the closest vector map
+#% keywords: vector
+#%end
+#%option G_OPT_V_INPUT
+#% key: points
+#% required: yes
+#%end
+#%option G_OPT_V_INPUT
+#% key: lines
+#% required: yes
+#%end
+#%option G_OPT_V_OUTPUT
+#% key: output
+#% required: yes
+#%end
+#%option
+#% key: max_dist
+#% type: double
+#% description: Maximum distance from points to river
+#% required: yes
+#% answer: 10
+#%end
+from __future__ import print_function
+
+# import system libraries
+import os
+import sys
+
+# import grass libraries
+from grass.script import core as gcore
+from grass.pygrass.messages import get_msgr
+from grass.pygrass.vector import VectorTopo
+
+
+if "GISBASE" not in os.environ:
+ print("You must be in GRASS GIS to run this program.")
+ sys.exit(1)
+
+
+def get_new_points(points, lines, output, maxdist=50):
+ skipped = []
+ ovwr = gcore.overwrite()
+ msgr = get_msgr()
+ points, pmset = points.split('@') if '@' in points else (points, '')
+ lines, lmset = lines.split('@') if '@' in lines else (lines, '')
+ with VectorTopo(points, mapset=pmset, mode='r') as pts:
+ cols = pts.table.columns.items() if pts.table else None
+ with VectorTopo(lines, mapset=lmset, mode='r') as lns:
+ with VectorTopo(output, mode='w', tab_cols=cols,
+ overwrite=ovwr) as out:
+ for pnt in pts:
+ line = lns.find['by_point'].geo(pnt, maxdist=maxdist)
+ if line is None:
+ msg = ("Not found any line in the radius of %.2f "
+ "for the point with cat: %d. The point "
+ "will be skipped!")
+ msgr.warning(msg % (maxdist, pnt.cat))
+ skipped.append(pnt.cat)
+ continue
+ # find the new point
+ newpnt, dist, _, _ = line.distance(pnt)
+ # get the attributes
+ attrs = (None if pnt.attrs is None
+ else pnt.attrs.values()[1:])
+ # write the new point in the new vector map
+ out.write(newpnt, attrs)
+ # save the changes on the output table
+ out.table.conn.commit()
+
+
+def main(opts, flgs):
+ get_new_points(opts['points'], opts['lines'], opts['output'],
+ float(opts['max_dist']))
+
+
+if __name__ == "__main__":
+ options, flags = gcore.parser()
+ sys.exit(main(options, flags))
Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.closest/r.green.hydro.closest.py
___________________________________________________________________
Added: svn:mime-type
+ text/x-python
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.delplants/Makefile
===================================================================
--- grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.delplants/Makefile (rev 0)
+++ grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.delplants/Makefile 2015-04-27 06:36:52 UTC (rev 65146)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../../../..
+
+PGM = r.green.hydro.delplants
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.delplants/Makefile
___________________________________________________________________
Added: svn:mime-type
+ text/x-makefile
Added: svn:eol-style
+ native
Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.delplants/r.green.hydro.delplants.html
___________________________________________________________________
Added: svn:mime-type
+ text/html
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.delplants/r.green.hydro.delplants.py
===================================================================
--- grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.delplants/r.green.hydro.delplants.py (rev 0)
+++ grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.delplants/r.green.hydro.delplants.py 2015-04-27 06:36:52 UTC (rev 65146)
@@ -0,0 +1,183 @@
+#!/usr/bin/env python
+# -- coding: utf-8 --
+#
+############################################################################
+#
+# MODULE: r.green.hydro.del_plants
+# AUTHOR(S): Pietro Zambelli & Giulia Garegnani
+# PURPOSE: Delete segmets where there is an existing plant
+# COPYRIGHT: (C) 2014 by the GRASS Development Team
+#
+# This program is free software under the GNU General Public
+# License (>=v2). Read the file COPYING that comes with GRASS
+# for details.
+#
+#############################################################################
+#
+#%Module
+#% description: Delete segmets where there is an existing plant
+#% overwrite: yes
+#%End
+#%option G_OPT_V_INPUT
+#% key: hydro
+#% description: Name of the vector map point of hydro plants.
+#% required: yes
+#%end
+#%option G_OPT_V_FIELD
+#% key: hydro_layer
+#% description: Name of the vector map layer of the hydro plant, with the following attributes: kind (water intake/turbine), discharge [l/s], working hours [hours], altitute[m], id point, id plant
+#% required: no
+#% answer: 1
+#%end
+#%option G_OPT_V_INPUT
+#% key: river
+#% description: Name of the vector map with the streams.
+#% required: yes
+#%end
+#%option G_OPT_V_OUTPUT
+#% description: Name of the vector map with the stream segments without plants.
+#% required: yes
+#%end
+#%option G_OPT_V_OUTPUT
+#% key: plants
+#% description: Name of the vector map with the stream segments already with plants.
+#% required: no
+#%end
+#%option
+#% key: hydro_kind_intake
+#% type: string
+#% description: Value contained in the column: hydro_kind that indicate the plant is an intake.
+#% required: no
+#% answer: intake
+#%end
+#%option
+#% key: hydro_kind_turbine
+#% type: string
+#% description: Value contained in the column: hydro_kind that indicate the plant is an intake.
+#% required: no
+#% answer: turbine
+#%end
+#%option G_OPT_R_ELEV
+#% required: yes
+#%end
+#%option G_OPT_V_MAP
+#% key: other
+#% description: Name of the vector map point of other plants such as irrigation, acqueducts, etc.
+#% required: no
+#%end
+#%option G_OPT_V_INPUT
+#% key: other_layer
+#% description: Name of the vector map layer of other plants, with the following attributes: kind (water intake/turbine), discharge [m3/year], id point, id plant
+#% required: no
+#% answer: 1
+#%end
+#%option
+#% key: other_kind_intake
+#% type: string
+#% description: Value contained in the column: other_kind that indicate the plant is an intake.
+#% required: no
+#% answer: intake
+#%end
+#%option
+#% key: other_kind_turbine
+#% type: string
+#% description: Value contained in the column: other_kind that indicate the plant is an intake.
+#% required: no
+#% answer: turbine
+#%end
+#%option
+#% key: efficiency
+#% type: double
+#% description: Plant efficiency.
+#% required: no
+#% answer: 0.8
+#%end
+#%flag
+#% key: d
+#% description: Debug with intermediate maps
+#%end
+from __future__ import print_function
+
+# import stadard libraries
+import os
+import atexit
+
+# import GRASS libraries
+from grass.exceptions import ParameterError
+from grass.script.core import parser, overwrite
+from grass.pygrass.vector import VectorTopo
+from grass.pygrass.raster import RasterRow
+from grass.pygrass.modules.shortcuts import vector as v
+from grass.pygrass.messages import get_msgr
+from grass.pygrass.utils import set_path
+
+set_path('r.green', 'libhydro', '..')
+set_path('r.green', 'libgreen', os.path.join('..', '..'))
+
+# finally import the module in the library
+from libgreen.utils import cleanup
+from libgreen.checkparameter import (check_required_columns,
+ exception2error)
+from libhydro.plant import read_plants, write_plants
+
+
+def main(opts, flgs):
+ TMPVECT = []
+ DEBUG = True if flgs['d'] else False
+ atexit.register(cleanup, vect=TMPVECT, debug=DEBUG)
+ # check input maps
+ rhydro = ['kind', 'discharge', 'id_point', 'id_plant']
+ rother = ['kind', 'discharge', 'id_point', 'id_plant']
+ ovwr = overwrite()
+
+ try:
+ hydro = check_required_columns(opts['hydro'], int(opts['hydro_layer']),
+ rhydro, 'hydro')
+ if opts['other']:
+ other = check_required_columns(opts['other'], opts['other_layer'],
+ rother, 'other')
+ else:
+ other = None
+ #minflow = check_float_or_raster(opts['minflow'])
+ except ParameterError as exc:
+ exception2error(exc)
+ # TODO: do we really nead the efficiency of the plant here?
+ #efficiency = float(opts['efficiency'])
+
+ # start working
+ hydro.open('r')
+ el, mset = (opts['elevation'].split('@') if '@' in opts['elevation']
+ else (opts['elevation'], ''))
+ elev = RasterRow(name=el, mapset=mset)
+ elev.open('r')
+ plants, skipped = read_plants(hydro, elev=elev,
+ restitution=opts['hydro_kind_turbine'],
+ intake=opts['hydro_kind_intake'])
+ hydro.close()
+ rvname, rvmset = (opts['river'].split('@') if '@' in opts['river']
+ else (opts['river'], ''))
+
+ vplants = opts['plants'] if opts['plants'] else 'tmpplants'
+ #FIXME: I try with tmpplants in my mapset and it doesn'work
+ if opts['plants'] == '':
+ TMPVECT.append(vplants)
+ with VectorTopo(rvname, rvmset, mode='r') as river:
+ write_plants(plants, vplants, river, elev, overwrite=ovwr)
+
+ if skipped:
+ for skip in skipped:
+ print("Plant: %r, Point: %r, kind: %r" % skip)
+ elev.close()
+
+ # compute a buffer around the plants
+ buff = vplants + 'buff'
+ v.buffer(input=vplants, type='line', output=buff, distance=0.1,
+ overwrite=ovwr)
+ TMPVECT.append(buff)
+ # return all the river segments that are not already with plants
+ v.overlay(flags='t', ainput=opts['river'], atype='line', binput=buff,
+ operator='not', output=opts['output'], overwrite=ovwr)
+
+
+if __name__ == "__main__":
+ main(*parser())
Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.delplants/r.green.hydro.delplants.py
___________________________________________________________________
Added: svn:mime-type
+ text/x-python
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.economic/Makefile
===================================================================
--- grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.economic/Makefile (rev 0)
+++ grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.economic/Makefile 2015-04-27 06:36:52 UTC (rev 65146)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../../../..
+
+PGM = r.green.hydro.economic
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.economic/Makefile
___________________________________________________________________
Added: svn:mime-type
+ text/x-makefile
Added: svn:eol-style
+ native
Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.economic/r.green.hydro.economic.html
___________________________________________________________________
Added: svn:mime-type
+ text/html
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.economic/r.green.hydro.economic.py
===================================================================
--- grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.economic/r.green.hydro.economic.py (rev 0)
+++ grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.economic/r.green.hydro.economic.py 2015-04-27 06:36:52 UTC (rev 65146)
@@ -0,0 +1,916 @@
+#!/usr/bin/env python
+# -- coding: utf-8 --
+#
+############################################################################
+#
+# MODULE: r.green.hydro.economic
+# AUTHOR(S): Sandro Sacchelli (BASH-concept), Pietro Zambelli (Python)
+# PURPOSE: Assess the hydro plant costs
+# COPYRIGHT: (C) 2014 by the GRASS Development Team
+#
+# This program is free software under the GNU General Public
+# License (>=v2). Read the file COPYING that comes with GRASS
+# for details.
+#
+#############################################################################
+#
+#%Module
+#% description: Compute the economic costs and values
+#% overwrite: yes
+#%End
+#%option G_OPT_V_MAP
+#% key: plant
+#% description: Name of the input vector map with the plants
+#% required: yes
+#%end
+#%option G_OPT_V_INPUT
+#% key: plant_layer
+#% description: Name of the vector map layer of other plants, with the following attributes: kind (water intake/turbine), discharge [m3/year], id point, id plant
+#% required: no
+#% answer: 1
+#%end
+#############################################################################
+# DEFINE COLUMNS OF V INPUT
+#%option
+#% key: plant_id_column
+#% type: string
+#% description: Column name with the plant id.
+#% required: no
+#% answer: plant_id
+#% guisection: Input columns
+#%end
+#%option
+#% key: plant_power_column
+#% type: string
+#% description: Column name with power value.
+#% required: no
+#% answer: power
+#% guisection: Input columns
+#%end
+#%option
+#% key: plant_head_column
+#% type: string
+#% description: Column name with head value.
+#% required: no
+#% answer: head
+#% guisection: Input columns
+#%end
+#%option
+#% key: plant_side_column
+#% type: string
+#% description: Column name with a string that define the side of the plant
+#% required: no
+#% answer: side
+#% guisection: Input columns
+#%end
+#%option
+#% key: plant_kind_column
+#% type: string
+#% description: Column name with the strings that define if it is conduct or penstock.
+#% required: no
+#% answer: kind
+#% guisection: Input columns
+#%end
+#%option
+#% key: plant_kind_intake
+#% type: string
+#% description: Value contained in the column 'kind' which corresponds to the derivation (conduct).
+#% required: no
+#% answer: conduct
+#% guisection: Input columns
+#%end
+#%option
+#% key: plant_kind_turbine
+#% type: string
+#% description: Value contained in the column 'kind' which corresponds to the penstock.
+#% required: no
+#% answer: penstock
+#% guisection: Input columns
+#%end
+
+
+#############################################################################
+# DEFINE COMPENSATION COSTS
+# provide raster
+
+#%option
+#% key: interest_rate
+#% type: double
+#% description: Interest rate value
+#% required: no
+#% answer: 0.03
+#% guisection: Compensation
+#%end
+#%option
+#% key: gamma_comp
+#% type: double
+#% description: Coefficient
+#% required: no
+#% answer: 1.25
+#% guisection: Compensation
+#%end
+#%option
+#% key: life
+#% type: double
+#% description: Life of the hydropower plant [year]
+#% required: no
+#% answer: 30
+#% guisection: Compensation
+#%end
+
+
+#%option G_OPT_R_INPUT
+#% key: landvalue
+#% description: Name of the raster map with the land value [currency/ha].
+#% required: no
+#% guisection: Compensation
+#%end
+#%option G_OPT_R_INPUT
+#% key: tributes
+#% description: Name of the raster map with the tributes [currency/ha].
+#% required: no
+#% guisection: Compensation
+#%end
+#%option G_OPT_R_INPUT
+#% key: stumpage
+#% description: Name of the raster map with the stumpage value [currency/ha].
+#% required: no
+#% guisection: Compensation
+#%end
+#%option G_OPT_R_INPUT
+#% key: rotation
+#% description: Name of the raster map with the rotation period per landuse type [year].
+#% required: no
+#% guisection: Compensation
+#%end
+#%option G_OPT_R_INPUT
+#% key: age
+#% description: Name of the raster map with the average age [year].
+#% required: no
+#% guisection: Compensation
+#%end
+
+# or rule to transform landuse categories in raster maps
+#%option G_OPT_R_INPUT
+#% key: landuse
+#% description: Name of the raster map with the landuse categories.
+#% required: no
+#% guisection: Compensation
+#%end
+
+# RULES
+#%option
+#% key: rules_landvalue
+#% type: string
+#% description: Rule file for the reclassification of the landuse to associate for each landuse a value [currency/ha].
+#% required: no
+#% guisection: Compensation
+#%end
+#%option
+#% key: rules_tributes
+#% type: string
+#% description: Rule file for the reclassification of the landuse to associate for each landuse a tribute rate [currency/ha].
+#% required: no
+#% guisection: Compensation
+#%end
+#%option
+#% key: rules_stumpage
+#% type: string
+#% description: Rule file for the reclassification of the landuse to associate for each landuse a stumpage value [currency/ha].
+#% required: no
+#% guisection: Compensation
+#%end
+#%option
+#% key: rules_rotation
+#% type: string
+#% description: Rule file for the reclassification of the landuse to associate for each landuse a rotation value [year].
+#% required: no
+#% guisection: Compensation
+#%end
+#%option
+#% key: rules_age
+#% type: string
+#% description: Rule file for the reclassification of the landuse to associate for each landuse an age [year].
+#% required: no
+#% guisection: Compensation
+#%end
+
+#############################################################################
+# DEFINE EXCAVATION COSTS
+
+#%option
+#% key: width
+#% type: double
+#% description: With of the excavation works [m]
+#% required: no
+#% answer: 2.
+#% guisection: Excavation
+#%end
+#%option
+#% key: depth
+#% type: double
+#% description:Depth of the excavation works [m]
+#% required: no
+#% answer: 2.
+#% guisection: Excavation
+#%end
+#%option
+#% key: slope_limit
+#% type: double
+#% description: Slope limit, above this limit the cost will equal to the maximum [degree]
+#% required: no
+#% answer: 50.
+#% guisection: Excavation
+#%end
+
+# provide raster maps
+
+#%option G_OPT_R_INPUT
+#% key: min_exc
+#% description: Minimum excavation costs [currency/mc]
+#% required: no
+#% guisection: Excavation
+#%end
+#%option G_OPT_R_INPUT
+#% key: max_exc
+#% description: Maximum excavation costs [currency/mc]
+#% required: no
+#% guisection: Excavation
+#%end
+#%option G_OPT_R_INPUT
+#% key: slope
+#% description: Slope raster map
+#% required: yes
+#% guisection: Excavation
+#%end
+
+
+
+# RULES
+#%option
+#% key: rules_min_exc
+#% type: string
+#% description: Rule file for the reclassification of the landuse to associate for each landuse a minimum excavation cost [currency/mc].
+#% required: no
+#% guisection: Excavation
+#%end
+#%option
+#% key: rules_max_exc
+#% type: string
+#% description: Rule file for the reclassification of the landuse to associate for each landuse a maximum excavation cost [currency/mc].
+#% required: no
+#% guisection: Excavation
+#%end
+
+#############################################################################
+# DEFINE ELECTROMECHANICAL COSTS
+#%option
+#% key: alpha_em
+#% type: double
+#% description: Electro-mechanical parameter costs, default values taken from Aggidis et al. 2010
+#% required: no
+#% answer: 0.56
+#% guisection: Electro-mechanical
+#%end
+#%option
+#% key: beta_em
+#% type: double
+#% description: Electro-mechanical parameter costs, default values taken from Aggidis et al. 2010
+#% required: no
+#% answer: -0.112
+#% guisection: Electro-mechanical
+#%end
+#%option
+#% key: gamma_em
+#% type: double
+#% description: Electro-mechanical parameter costs, default values taken from Aggidis et al. 2010
+#% required: no
+#% answer: 15600.0
+#% guisection: Electro-mechanical
+#%end
+#%option
+#% key: const_em
+#% type: double
+#% description: Electro-mechanical parameter costs constant value, default values taken from Aggidis et al. 2010
+#% required: no
+#% answer: 0.
+#% guisection: Electro-mechanical
+#%end
+
+#############################################################################
+# DEFINE SUPPLY AND INSTALLATION COSTS
+#%option
+#% key: lc_pipe
+#% type: double
+#% description: Supply and installation linear cost for the pipeline [currency/m]
+#% answer: 310.
+#% guisection: Supply & Installation
+#%end
+#%option
+#% key: lc_electro
+#% type: double
+#% description: Supply and installation linear cost for the electroline [currency/m]
+#% answer: 250.
+#% guisection: Supply & Installation
+#%end
+#%option G_OPT_V_MAP
+#% key: electro
+#% description: Name of the vector map with the electric grid
+#% required: yes
+#%end
+#%option G_OPT_V_INPUT
+#% key: electro_layer
+#% description: Name of the vector map layer of other plants, with the following attributes: kind (water intake/turbine), discharge [m3/year], id point, id plant
+#% required: no
+#% answer: 1
+#%end
+
+#############################################################################
+# DEFINE POWER STATION COSTS
+#%option
+#% key: alpha_station
+#% type: double
+#% description: Power station costs are assessd as a fraction of the Electro-mechanical costs
+#% answer: 0.52
+#% guisection: Power station
+#%end
+
+#############################################################################
+# DEFINE INLET COSTS
+#%option
+#% key: alpha_inlet
+#% type: double
+#% description: Inlet costs are assessd as a fraction of the Electro-mechanical costs
+#% answer: 0.38
+#% guisection: Inlet
+#%end
+
+#############################################################################
+# DEFINE OTHER COSTS
+#%option
+#% key: grid
+#% type: double
+#% description: Cost for grid connection
+#% answer: 50000
+#% guisection: Other
+#%end
+#%option
+#% key: general
+#% type: double
+#% description: Factor for general expenses
+#% answer: 0.15
+#% guisection: Other
+#%end
+#%option
+#% key: hindrances
+#% type: double
+#% description: Factor for hindrances expenses
+#% answer: 0.1
+#% guisection: Other
+#%end
+# DEFINE MAINTENANCE COSTS
+#%option
+#% key: cost_maintenance_per_kw
+#% type: double
+#% description: Maintenace costs per kW
+#% answer: 7000.
+#% guisection: Maintenance
+#%end
+#%option
+#% key: alpha_maintenance
+#% type: double
+#% description: alpha coefficient to assess the maintenance costs
+#% answer: 0.05
+#% guisection: Maintenance
+#%end
+#%option
+#% key: beta_maintenance
+#% type: double
+#% description: beta coefficient to assess the maintenance costs
+#% answer: 0.45
+#% guisection: Maintenance
+#%end
+#%option
+#% key: const_maintenance
+#% type: double
+#% description: constant to assess the maintenance costs
+#% answer: 0.
+#% guisection: Maintenance
+#%end
+
+# DEFINE REVENUES PARAMETERS
+#%option
+#% key: energy_price
+#% type: double
+#% description: Energy price per kW [currency/kW]
+#% answer: 0.09
+#% guisection: Revenues
+#%end
+#%option
+#% key: eta
+#% type: double
+#% description: Energy price per kW [currency/kW]
+#% answer: 0.81
+#% guisection: Revenues
+#%end
+#%option
+#% key: operative_hours
+#% type: double
+#% description: The number of operative hours per year [hours/year]
+#% answer: 6500.
+#% guisection: Revenues
+#%end
+#%option
+#% key: alpha_revenue
+#% type: double
+#% description: coefficient to transform installed power to mean
+#% answer: 0.5
+#% guisection: Revenues
+#%end
+#%option
+#% key: const_revenue
+#% type: double
+#% description: constant to assess the revenues
+#% answer: 0.
+#% guisection: Revenues
+#%end
+
+#############################################################################
+# DEFINE OUTPUTS
+#%option G_OPT_V_OUTPUT
+#% key: output
+#% description: Name of the output vector map
+#% required: yes
+#%end
+
+## COSTS
+#%option G_OPT_R_OUTPUT
+#% key: compensation
+#% description: Raster map with the compensation values
+#% required: no
+#%end
+#%option G_OPT_R_OUTPUT
+#% key: excavation
+#% description: Raster map with the excavation costs
+#% required: no
+#%end
+
+## VALUES
+#%option G_OPT_R_OUTPUT
+#% key: upper
+#% description: Raster map with the compensation values
+#% required: no
+#%end
+
+#############################################################################
+#%rules
+#% required: landvalue, rules_landvalue
+#% requires: rules_landvalue, landuse
+#% required: tributes, rules_tributes
+#% requires: rules_tributes, landuse
+#% required: stumpage, rules_stumpage
+#% requires: rules_stumpage, landuse
+#% required: tributes, rules_rotation
+#% requires: rules_rotation, landuse
+#% required: tributes, rules_age
+#% requires: rules_age, landuse
+#% required: min_exc, rules_min_exc
+#% requires: rules_min_exc, landuse
+#% required: max_exc, rules_max_exc
+#% requires: rules_max_exc, landuse
+#%end
+#exclusive: at most one of the options may be given
+#required: at least one of the options must be given
+#requires: if the first option is given, at least one of the subsequent options must also be given
+#requires_all: if the first option is given, all of the subsequent options must also be given
+#excludes: if the first option is given, none of the subsequent options may be given
+#collective: all or nothing; if any option is given, all must be given
+from __future__ import print_function
+
+import os
+import atexit
+import random
+
+import numpy as np
+import numexpr as ne
+
+from grass.exceptions import ParameterError
+from grass.script.core import parser, run_command, overwrite
+from grass.pygrass.modules.shortcuts import raster as r
+from grass.pygrass.modules.shortcuts import vector as v
+
+from grass.pygrass.utils import set_path
+from grass.pygrass.raster import RasterRow
+from grass.pygrass.vector import VectorTopo, sql
+from grass.pygrass.vector.basic import Cats
+
+from grass.pygrass.messages import get_msgr
+
+# set python path to the shared r.green libraries
+set_path('r.green', 'libhydro', '..')
+set_path('r.green', 'libgreen', os.path.join('..', '..'))
+
+from libgreen.utils import cleanup
+from libgreen.checkparameter import check_required_columns, exception2error
+from libhydro.plant import read_plants, write_structures
+
+
+def rname(base):
+ return '%s_%03d' % (base, random.randint(0, 1000))
+
+
+def check_raster_or_landuse(opts, params):
+ """Return a list of raster name, i f necessary the rasters are generated
+ using the file with the rules to convert a landuse map to the raster
+ that it is needed.
+ """
+ msgr = get_msgr()
+ rasters = []
+ for par in params:
+ if opts[par]:
+ rasters.append(opts[par])
+ elif opts['rules_%s' % par] and opts['landuse']:
+ output = rname(par)
+ msg = 'Creating: {out} using the rules: {rul}'
+ msgr.verbose(msg.format(out=output, rul='rules_%s' % par))
+ r.reclass(input=opts['landuse'], output=output,
+ rules=opts['rules_%s' % par])
+ rasters.append(output)
+ else:
+ msg = '{par} or rule_{par} are required'
+ raise ParameterError(msg.format(par=par))
+ return rasters
+
+
+def upper_value(upper, stu, lan, rot, age, irate, overwrite=False):
+ """Compute the upper value of a land.
+ """
+ expr = ("{upper} = ({stu} + {lan}) / ((1 + {irate})^({rot} - {age}) )"
+ " - {lan}")
+ r.mapcalc(expr.format(upper=upper, stu=stu, lan=lan, irate=irate,
+ rot=rot, age=age), overwrite=overwrite)
+
+
+def compensation_cost(comp, lan, tri, upper,
+ irate, gamma, life, width, overwrite=False):
+ """Compute the compensation raster map costs"""
+ expr = ("{comp} = (({lan} + {tri} * "
+ "(1 + {irate})^{life}/({irate} * (1 + {irate}))) "
+ "* {gamma} + {upper}) * {width} * nsres() / 10000")
+ r.mapcalc(expr.format(comp=comp, lan=lan, tri=tri, upper=upper,
+ irate=irate, gamma=gamma, life=life, width=width),
+ overwrite=overwrite)
+
+
+def excavation_cost(exc, excmin, excmax, slope,
+ slim, width, depth, overwrite=False):
+ """Compute the excavation cost"""
+ expr = ("{exc} = if({slope} < {slim}, "
+ "({excmin} + ({excmax} - {excmin}) / {slim} * {slope})"
+ "* {width} * {depth} * nsres(),"
+ "{excmax} * {width} * {depth} * nsres())")
+ r.mapcalc(expr.format(exc=exc, slope=slope, excmin=excmin, excmax=excmax,
+ slim=slim, width=width, depth=depth),
+ overwrite=overwrite)
+
+def vmapcalc2(vmap, vlayer, cname, ctype, expr, overwrite=False):
+ v.db_addcolumn(map=vmap, layer=vlayer, columns=[(cname, ctype)],
+ overwrite=overwrite)
+ v.db_update(map=vmap, layer=vlayer, column=cname, query_column=expr,
+ overwrite=overwrite)
+
+
+def get_cnames(expr,
+ _names_cache=ne.utils.CacheDict(256),
+ _numexpr_cache=ne.utils.CacheDict(256), **kwargs):
+ if not isinstance(expr, (str, unicode)):
+ raise ValueError("must specify expression as a string")
+ # Get the names for this expression
+ context = ne.necompiler.getContext(kwargs, frame_depth=1)
+ expr_key = (expr, tuple(sorted(context.items())))
+ if expr_key not in _names_cache:
+ _names_cache[expr_key] = ne.necompiler.getExprNames(expr.strip(),
+ context)
+ names, ex_uses_vml = _names_cache[expr_key]
+ return names
+
+
+def vcolcalc(vname, vlayer, ctype, expr, condition=lambda x: x is None,
+ notfinitesubstitute=None, **kwargs):
+ equal = expr.index('=')
+ if equal < 0:
+ raise
+ cname = expr[:equal].strip()
+ expr = expr[(equal + 1):].strip()
+ cnames = get_cnames(expr)
+ with VectorTopo(vname, layer=vlayer, mode='r') as vect:
+ if vect.table is None:
+ msg = 'Vector: {vname} is without table.'
+ raise TypeError(msg.format(vname=vname))
+
+ cols = vect.table.columns
+ # check if the column in the expressions exist or not
+ for col in cnames:
+ if col not in cols:
+ msg = 'Vector: {vname} has not column: {col} in layer: {layer}'
+ raise TypeError(msg.format(vname=vname, col=col,
+ layer=vlayer))
+ if cname not in cols:
+ vect.table.columns.add(cname, ctype)
+ # extract value from the attribute table
+ if cnames[0] != vect.table.key:
+ cnames.insert(0, vect.table.key)
+ # TODO: find a more elegant way to do it
+ sql = vect.table.filters.select(', '.join(cnames)).where(' is not NULL AND '.join(cnames)).get_sql()[:-1] + ' is not NULL;'
+ data = np.array(list(vect.table.execute(sql)))
+ # create a dictionary with local variables
+ lvars = {col: array for col, array in zip(cnames, data.T)}
+ # compute the result with numexpr
+ res = ne.evaluate(expr, local_dict=lvars, **kwargs)
+ if notfinitesubstitute is not None:
+ res[~np.isfinite(res)] = notfinitesubstitute
+ # save the result to the vector point
+ cur = vect.table.conn.cursor()
+ upstr = 'UPDATE {tname} SET {cname}=? WHERE {key} == ?;'
+ cur.executemany(upstr.format(tname=vect.table.name, cname=cname,
+ key=vect.table.key),
+ zip(res, lvars[vect.table.key]))
+ # save changes
+ vect.table.conn.commit()
+
+
+def electromechanical_cost(vname, power, head,
+ gamma=15600., alpha=0.56, beta=0.112, const=0.,
+ vlayer=1, cname='em_cost',
+ ctype='double precision',
+ overwrite=False):
+ expr = "{cname} = {gamma} * {power}**{alpha} * {head}**{beta} + {const}"
+ vcolcalc(vname, vlayer, ctype, notfinitesubstitute=0.,
+ expr=expr.format(cname=cname, gamma=gamma, power=power,
+ alpha=alpha, head=head, beta=beta, const=const))
+
+
+def col_exist(vname, cname, ctype='double precision', vlayer=1, create=False):
+ with VectorTopo(vname, layer=vlayer, mode='r') as vect:
+ if vect.table is None:
+ msg = 'Vector: {vname} is without table.'
+ raise TypeError(msg.format(vname=vname))
+ res = cname in vect.table.columns
+ if res is False:
+ vect.table.columns.add(cname, ctype)
+ return res
+
+
+def linear_cost(vname, cname='lin_cost', alpha=310., length='length', vlayer=1,
+ ctype='double precision', overwrite=False):
+ # check if length it is alread in the db
+ if not col_exist(vname, 'length', create=True):
+ v.to_db(map=vname, type='line', layer=vlayer, option='length',
+ columns='length')
+ expr = "{cname} = {alpha} * {length}"
+ vcolcalc(vname, vlayer, ctype, notfinitesubstitute=0.,
+ expr=expr.format(cname=cname, alpha=alpha, length=length))
+
+
+def get_electro_length(opts):
+ # open vector plant
+ with VectorTopo(opts['plant'], layer=int(opts['plant_layer']),
+ mode='r') as vect:
+ kcol = opts['plant_kind_column']
+ ktype = opts['plant_kind_turbine']
+ # check if electro_length it is alredy in the table
+ if 'electro_length' not in vect.table.columns:
+ vect.table.columns.add('electro_length', 'double precision')
+ # open vector map with the existing electroline
+ with VectorTopo(opts['electro'], layer=int(opts['electro_layer']),
+ mode='r') as electro:
+ for line in vect:
+ if line.attrs[kcol] == ktype:
+ # the turbine is the last point of the penstock
+ turbine = line[-1]
+ # find the closest electro line
+ eline = electro.find['by_point'].geo(turbine, maxdist=1e6)
+ dist = eline.distance(turbine)
+ line.attrs['electro_length'] = dist.dist
+ else:
+ line.attrs['electro_length'] = 0.
+ vect.table.conn.commit()
+
+
+def get_gamma_NPV(r=0.03, y=30):
+ """gamma it is a coefficient define as:
+
+ $\gamma = 1 - \frac{1 - (1+r)^y}{r(1+r)^y}$
+
+ with $r$ as the interest rate (default value: 0.03) and
+ $y$ as the number of years of the plant [years] (default value: 30);
+ """
+ return 1 - (1 - (1 + r)**y) / (r * (1 + r)**y)
+
+
+def group_by(vinput, voutput, isolate=None, aggregate=None,
+ function='sum', vtype='lines',
+ where='', group_by=None, linput=1, loutput=1):
+ with VectorTopo(vinput, mode='r') as vin:
+ columns = ['cat', ]
+ vincols = vin.table.columns
+ types = ['PRIMARY KEY', ]
+ siso = ''
+ if isolate is not None:
+ columns += list(isolate)
+ types += [vincols[c] for c in isolate]
+ siso = ', '.join(isolate)
+ sagg = ''
+ if aggregate is not None:
+ ct = '{func}({col}) as {col}'
+ sagg = ', '.join([ct.format(func=function, col=col)
+ for col in aggregate])
+ columns += list(aggregate)
+ types += [vincols[c] for c in aggregate]
+
+ scols = "%s, %s" % (siso, sagg) if siso and sagg else siso + sagg
+ base = "SELECT {cols} FROM {tin}"
+ bwhere = " WHERE %s" % where if where else ''
+ bgroup = " GROUP BY %s" % ', '.join(group_by) if group_by else ''
+ grp = (base + bwhere + bgroup + ';').format(cols=scols,
+ tin=vin.table.name,
+ tout=voutput)
+ bqry = "SELECT {cat} FROM {tin} WHERE {cond};"
+ qry = bqry.format(cat=vin.table.key, tin=vin.table.name,
+ cond=' AND '.join(['%s=?' % g for g in group_by]))
+ selcols = columns[1:]
+ gindexs = [selcols.index(col) for col in group_by]
+ insrt = sql.INSERT.format(tname=voutput,
+ values=','.join(['?', ] * len(columns)))
+ with VectorTopo(voutput, mode='w', tab_name=voutput, layer=loutput,
+ tab_cols=list(zip(columns, types)), link_key='cat',
+ overwrite=True) as vout:
+ cur = vout.table.conn.cursor()
+ ncat = 1
+ #import ipdb; ipdb.set_trace()
+ # TODO: why do I need to use list(cur) and I can not iterate directly
+ for row in list(cur.execute(grp)):
+ # add the new line to table
+ print(row)
+ cur.execute(insrt, tuple([ncat, ] + list(row)))
+ for cat in cur.execute(qry, tuple([row[i] for i in gindexs])):
+ for line in vin.cat(cat[0], vtype):
+ # set the new category
+ category = Cats(line.c_cats)
+ category.reset()
+ category.set(ncat, loutput)
+ # write geometry
+ vout.write(line, set_cats=False)
+ ncat += 1
+ vout.table.conn.commit()
+
+
+def main(opts, flgs):
+ # check or generate raster map from rules files
+ ecovalues = ['landvalue', 'tributes', 'stumpage', 'rotation', 'age',
+ 'min_exc', 'max_exc']
+ (lan, tri, stu, rot, age,
+ excmin, excmax) = check_raster_or_landuse(opts, ecovalues)
+ upper = opts['upper'] if opts['upper'] else 'tmp_upper'
+ comp = opts['compensation'] if opts['compensation'] else 'tmp_compensation'
+ exc = opts['excavation'] if opts['excavation'] else 'tmp_excavation'
+ vlayer = int(opts['plant_layer'])
+
+ # read common scalar parameters
+ irate = float(opts['interest_rate'])
+ life = float(opts['life'])
+ width = float(opts['width'])
+ depth = float(opts['depth'])
+ slim = float(opts['slope_limit'])
+ overw = overwrite()
+
+ # RASTERS
+ # Start computing the raster map of the costs
+ # Upper value
+ upper_value(upper, stu, lan, rot, age, irate, overw)
+ # Compensation raster costs
+ compensation_cost(comp, lan, tri, upper,
+ irate, float(opts['gamma_comp']), life, width, overw)
+ # Excavation raster costs
+ excavation_cost(exc, excmin, excmax, opts['slope'],
+ slim, width, depth, overw)
+ # TODO: extra cost when crossing roads and rivers are missing
+
+ # VECTOR
+ # add columns with costs from rasters
+ # add compensation costs
+ v.rast_stats(map=opts['plant'], layer=vlayer, flags='c',
+ raster=comp, column_prefix='comp_cost', method='sum')
+ # add excavation costs
+ v.rast_stats(map=opts['plant'], layer=vlayer, flags='c',
+ raster=exc, column_prefix='exc_cost', method='sum')
+
+ # add elecro-mechanical costs
+ electromechanical_cost(opts['plant'],
+ power=opts['plant_power_column'],
+ head=opts['plant_head_column'],
+ gamma=float(opts['gamma_em']),
+ alpha=float(opts['alpha_em']),
+ beta=float(opts['beta_em']),
+ const=float(opts['const_em']),
+ vlayer=vlayer, cname='em_cost',
+ ctype='double precision',
+ overwrite=overw)
+
+ # add linear cost for pipeline
+ linear_cost(vname=opts['plant'], cname='lin_pipe_cost',
+ alpha=float(opts['lc_pipe']), vlayer=vlayer,
+ ctype='double precision', overwrite=overw)
+
+ # add linear for for electroline
+ get_electro_length(opts)
+ linear_cost(vname=opts['plant'], cname='lin_electro_cost',
+ alpha=float(opts['lc_electro']), length='electro_length',
+ vlayer=vlayer, ctype='double precision', overwrite=overw)
+
+ xcost = "{cname} = {alpha} * {em}"
+ # add power station costs
+ vcolcalc(vname=opts['plant'], vlayer=vlayer,
+ ctype='double precision',
+ expr=xcost.format(cname='station_cost', em='em_cost',
+ alpha=opts['alpha_station']))
+ # add inlet costs
+ vcolcalc(vname=opts['plant'], vlayer=vlayer,
+ ctype='double precision', notfinitesubstitute=0.,
+ expr=xcost.format(cname='inlet_cost', em='em_cost',
+ alpha=opts['alpha_inlet']))
+ # add total inlet costs
+ # TODO: to be check to avoid to count cost more than one time I have moltiplied by 0.5
+ tot = ('tot_cost = (comp_cost_sum + em_cost + '
+ 'lin_pipe_cost + lin_electro_cost + '
+ 'station_cost + inlet_cost + {grid}*0.5) * '
+ '(1 + {general} + {hindrances})')
+ vcolcalc(vname=opts['plant'], vlayer=vlayer,
+ ctype='double precision', notfinitesubstitute=0.,
+ expr=tot.format(grid=opts['grid'], general=opts['general'],
+ hindrances=opts['hindrances']))
+
+ # TODO: produce a new vector map (output), with the conduct + penstock in
+ # a unique line and sum the costs grouping by intake_id and side
+ # SELECT {key} FROM {tname}
+ group_by(opts['plant'], opts['output'],
+ isolate=['intake_id', 'plant_id', 'side', 'power',
+ 'head', 'discharge'],
+ aggregate=['tot_cost', ],
+ function='sum',
+ group_by=['intake_id', 'side'])
+
+ """
+ where these values (3871.2256 and -0.45) are coming from?
+ import numpy as np
+ from scipy import stats
+
+ power= np.array([50., 100., 200., 400., 600., 1000., 5000.])
+ maint = np.array([707., 443., 389., 261., 209., 163., 88.])
+
+ plt.plot(np.log(power), np.log(maint))
+ plt.show()
+ slope, intercept, r_value, p_value, std_err = stats.linregress(np.log(power), np.log(maint))
+ #slope -0.45000431409701719
+ #intercept 8.2613264284076049
+ np.exp(intercept) * power ** slope
+
+ """
+# maint = "{cname} = {alpha} * {power} ** (1 + {beta}) + {const}"
+ maint = "{cname} = {cost_per_kW} * {power} * {alpha}) + {const}"
+ # compute yearly maintenance costs
+ vcolcalc(vname=opts['output'], vlayer=vlayer,
+ ctype='double precision', notfinitesubstitute=0.,
+ expr=maint.format(cname='maintenance',
+ cost_per_kw=opts['cost_maintenance_per_kw'],
+ alpha=opts['alpha_maintenance'],
+ power=opts['plant_power_column'],
+ beta=opts['beta_maintenance'],
+ const=opts['const_maintenance']))
+
+ # compute yearly revenues
+ rev = "{cname} = {eta} * {power} * {eprice} * {ophours} * {alpha} + {const}"
+ vcolcalc(vname=opts['output'], vlayer=vlayer,
+ ctype='double precision', notfinitesubstitute=0.,
+ expr=rev.format(cname='revenue',
+ eta=opts['eta'],
+ power=opts['plant_power_column'],
+ eprice=opts['energy_price'],
+ ophours=opts['operative_hours'],
+ alpha=opts['alpha_revenue'],
+ const=opts['const_revenue']))
+
+ # compute the Net Present Value
+ npv = "{cname} = {gamma} * ({revenue} - {maintenance}) - {tot}"
+ gamma_npv = get_gamma_NPV(irate, life)
+ vcolcalc(vname=opts['output'], vlayer=vlayer,
+ ctype='double precision', notfinitesubstitute=0.,
+ expr=npv.format(cname='NPV',
+ gamma=gamma_npv,
+ revenue='revenue',
+ maintenance='maintenance',
+ tot='tot_cost'))
+
+
+if __name__ == "__main__":
+ main(*parser())
Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.economic/r.green.hydro.economic.py
___________________________________________________________________
Added: svn:mime-type
+ text/x-python
Added: svn:eol-style
+ native
Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.html
___________________________________________________________________
Added: svn:mime-type
+ text/html
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.legal/Makefile
===================================================================
--- grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.legal/Makefile (rev 0)
+++ grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.legal/Makefile 2015-04-27 06:36:52 UTC (rev 65146)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../../../..
+
+PGM = r.green.hydro.legal
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.legal/Makefile
___________________________________________________________________
Added: svn:mime-type
+ text/x-makefile
Added: svn:eol-style
+ native
Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.legal/r.green.hydro.legal.html
___________________________________________________________________
Added: svn:mime-type
+ text/html
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.legal/r.green.hydro.legal.py
===================================================================
--- grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.legal/r.green.hydro.legal.py (rev 0)
+++ grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.legal/r.green.hydro.legal.py 2015-04-27 06:36:52 UTC (rev 65146)
@@ -0,0 +1,202 @@
+#!/usr/bin/env python
+# -- coding: utf-8 --
+#
+############################################################################
+#
+# MODULE: r.green.hydro.legal
+# AUTHOR(S): Giulia Garegnani
+# PURPOSE: Calculate the optimal position of a plant along a river
+# following legal constrains
+# COPYRIGHT: (C) 2014 by the GRASS Development Team
+#
+# This program is free software under the GNU General Public
+# License (>=v2). Read the file COPYING that comes with GRASS
+# for details.
+#
+#############################################################################
+#
+
+#%module
+#% description: Hydropower energy potential with legal constrains
+#% keywords: raster
+#% overwrite: yes
+#%end
+#%option
+#% key: discharge
+#% type: string
+#% gisprompt: old,cell,raster
+#% key_desc: name
+#% description: Name of river discharge [l/s]
+#% required: yes
+#%end
+#%option
+#% key: mvf
+#% type: string
+#% gisprompt: old,cell,raster
+#% key_desc: name
+#% description: Name of river mvf discharge [l/s]
+#% required: yes
+#%end
+#%option
+#% key: river
+#% type: string
+#% gisprompt: old,vector,vector
+#% key_desc: name
+#% description: Name of vector map with interested segments of rivers
+#% required: yes
+#%end
+#%option
+#% key: elevation
+#% type: string
+#% gisprompt: old,cell,raster
+#% key_desc: name
+#% description: Name of dtm [m]
+#% required: yes
+#%end
+#%flag
+#% key: d
+#% description: Debug with intermediate maps
+#%end
+#%flag
+#% key: c
+#% description: Clean vector lines
+#%end
+#%option
+#% key: len_plant
+#% type: double
+#% key_desc: name
+#% description: maximum length plant [m]
+#% required: yes
+#% answer: 100000000
+#%end
+#%option
+#% key: len_min
+#% type: double
+#% key_desc: name
+#% description: minimum length plant [m]
+#% required: yes
+#% answer: 10
+#%end
+#%option
+#% key: distance
+#% type: double
+#% description: minimum distance among plants.
+#% required: yes
+#% answer: 0.5
+#%end
+#%option
+#% key: area
+#% type: string
+#% gisprompt: old,vector,vector
+#% key_desc: name
+#% description: areas to be excluded.
+#% required: no
+#%end
+#%option
+#% key: buff
+#% type: double
+#% description: buffer for areas to be excluded.
+#% required: no
+#%end
+#%option
+#% key: output_plant
+#% type: string
+#% key_desc: name
+#% description: Name of output vector with potential segments
+#% required: no
+#%end
+#%option
+#% key: output_point
+#% type: string
+#% key_desc: name
+#% description: Name of output vector with potential intakes and restitution
+#% required: yes
+#%END
+
+# import system libraries
+from __future__ import print_function
+import os
+import sys
+import atexit
+
+# import grass libraries
+from grass.script import core as gcore
+from grass.pygrass.utils import set_path
+from grass.pygrass.messages import get_msgr
+from grass.script import mapcalc
+
+# r.green lib
+set_path('r.green', 'libhydro', '..')
+set_path('r.green', 'libgreen', os.path.join('..', '..'))
+# finally import the module in the library
+from libgreen.utils import cleanup
+
+
+if "GISBASE" not in os.environ:
+ print("You must be in GRASS GIS to run this program.")
+ sys.exitd(1)
+
+
+def main(opts, flgs):
+ DEBUG = False
+ TMPRAST = []
+ atexit.register(cleanup, rast=TMPRAST, debug=DEBUG)
+ dtm = options['elevation']
+ river = options['river'] # raster
+ discharge = options['discharge'] # vec
+ mvf = options['mvf'] # vec
+ len_plant = float(options['len_plant'])
+ len_min = float(options['len_min'])
+ distance = float(options['distance'])
+ output_plant = options['output_plant']
+ output_point = options['output_point']
+ area = options['area']
+ buff = options['buff']
+ DEBUG = flags['d']
+ c = flags['c']
+ msgr = get_msgr()
+
+ TMPRAST = []
+ if not gcore.overwrite():
+ for m in TMPRAST:
+ if gcore.find_file(m)['name']:
+ msgr.fatal(_("Temporary raster map %s exists") % (m))
+ #FIXME:check if it works for vectors
+
+ if area:
+ if buff:
+ gcore.run_command('v.buffer',
+ input=area,
+ output='buff_area',
+ distance=buff)
+ area = 'buff_area'
+
+ gcore.run_command('v.overlay',
+ ainput=river,
+ binput=area,
+ operator='not',
+ output='new_river')
+ river = 'new_river'
+
+ command = ('q_leg= %s - %s') % (discharge, mvf)
+ mapcalc(command, overwrite=True)
+
+ gcore.run_command('r.green.hydro.optimal',
+ flags=c,
+ discharge=discharge,
+ river=river,
+ elevation=dtm,
+ len_plant=len_plant,
+ output_plant=output_plant,
+ output_point=output_point,
+ distance=distance,
+ len_min=len_min,
+ efficiency=1)
+
+ #TODO: add the possibility to exclude vector segments of rivers
+ # by using the attribute of the river input
+
+if __name__ == "__main__":
+ atexit.register(cleanup)
+ options, flags = gcore.parser()
+ sys.exit(main(options, flags))
Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.legal/r.green.hydro.legal.py
___________________________________________________________________
Added: svn:mime-type
+ text/x-python
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.optimal/Makefile
===================================================================
--- grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.optimal/Makefile (rev 0)
+++ grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.optimal/Makefile 2015-04-27 06:36:52 UTC (rev 65146)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../../../..
+
+PGM = r.green.hydro.optimal
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.optimal/Makefile
___________________________________________________________________
Added: svn:mime-type
+ text/x-makefile
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.optimal/availablestreams.jpg
===================================================================
(Binary files differ)
Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.optimal/availablestreams.jpg
___________________________________________________________________
Added: svn:mime-type
+ image/jpeg
Added: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.optimal/potentialplants.jpg
===================================================================
(Binary files differ)
Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.optimal/potentialplants.jpg
___________________________________________________________________
Added: svn:mime-type
+ image/jpeg
Added: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.optimal/r.green.hydro.optimal.html
===================================================================
--- grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.optimal/r.green.hydro.optimal.html (rev 0)
+++ grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.optimal/r.green.hydro.optimal.html 2015-04-27 06:36:52 UTC (rev 65146)
@@ -0,0 +1,55 @@
+<h2>DESCRIPTION</h2>
+<em>r.green.hydro.optimal</em> calculates the hydropower energy potential for each basin. Deciding the range of plant length and the distance between plants, the module returns a vector file with the intakes and the restitution of the river which maximize the power that can be produced.
+
+<h2>NOTES</h2>
+The power is defined as P=ηρgQΔh where :<br>
+η is the efficiency of the plant <br>
+ρ the density of water (1000 kg/m<sup>3</sup>) <br>
+g the gravity term (9,81 m/s<sup>2</sup>) <br>
+Q the discharge of the river<br>
+Δh the gross head of the considered segment<br> <br>
+
+The module maximizes the power over a given range by a brute-force search in order to examine all possible arrangements of Q and Δh. Thus, the potential segments can be shorter than the maximum length plant chosen because it depends on the maximization of the product Q * Δh. <br> <br>
+The three input files have to give : the river considered, the discharge for each point of this river and the dtm to calculate the gross head. <br> <br>
+For each potential segment, the potential power is given in kW in attribute. The categories are divided in range of segment lengths. Arms of river are distinguished by a number, as well as the plants on these arms.
+
+
+<h2>EXAMPLE</h2>
+These examples are based on the Pnam file which refers to the Gesso and Vermenagna Valley in Piedmont, Italy.<br><br>
+
+Here is the vector file availablestreams of the interested streams. The river segments already exploited by an existing plant do not appear in the file.<br>
+
+<center>
+<img src="availablestreams.jpg" alt="availablestreams"><br>
+Vector map availablestreams
+</center><br>
+
+<div class="code"><pre>r.green.hydro.optimal discharge=discharge river=availablestreams dtm=elevation len_plant=800 distance=200 output_plant=potentialsegments output_point=potentialplants efficiency=1<br>
+d.vect map= potentialplants color=red<br>
+d.vect map= potentialsegments color=blue<br></pre></div><br>
+
+This command calculates the energy potential for a plant length range from 10 to 800 m and a distance between plants of 200m.<br><br>
+
+<center>
+<img src="potentialplants.jpg" alt="potentialplants"><br>
+Superposition of the potential segments vector file (potentialsegments, in blue) and the potential intakes and restitution vector file (potentialplants, in red)
+</center>
+
+
+<h2>SEE ALSO</h2>
+<em>
+<a href="r.green.hydro.structure.html">r.green.hydro.structure</a><br>
+<a href="r.green.hydro.technical.html">r.green.hydro.technical</a>
+</em>
+
+<h2>REFERENCE</h2>
+
+
+<h2>AUTHORS</h2>
+
+
+
+</div>
+</body>
+</html>
+
Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.optimal/r.green.hydro.optimal.html
___________________________________________________________________
Added: svn:mime-type
+ text/html
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.optimal/r.green.hydro.optimal.py
===================================================================
--- grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.optimal/r.green.hydro.optimal.py (rev 0)
+++ grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.optimal/r.green.hydro.optimal.py 2015-04-27 06:36:52 UTC (rev 65146)
@@ -0,0 +1,188 @@
+#!/usr/bin/env python
+# -- coding: utf-8 --
+#
+############################################################################
+#
+# MODULE: r.green.hydro.optimal
+# AUTHOR(S): Giulia Garegnani
+# PURPOSE: Calculate the optimal position of a plant along a river
+# COPYRIGHT: (C) 2014 by the GRASS Development Team
+#
+# This program is free software under the GNU General Public
+# License (>=v2). Read the file COPYING that comes with GRASS
+# for details.
+#
+#############################################################################
+#
+#%module
+#% description: Calculate the hydropower energy potential for each basin
+#% keyword: raster
+#%end
+#%option
+#% key: discharge
+#% type: string
+#% gisprompt: old,cell,raster
+#% key_desc: name
+#% description: Name of river discharge [m^3/s]
+#% required: yes
+#%end
+#%option
+#% key: river
+#% type: string
+#% gisprompt: old,vector,vector
+#% key_desc: name
+#% description: Name of vector map with interested segments of rivers
+#% required: yes
+#%end
+#%option
+#% key: elevation
+#% type: string
+#% gisprompt: old,cell,raster
+#% key_desc: name
+#% description: Name of elevation raster map [m]
+#% required: yes
+#%end
+#%flag
+#% key: d
+#% description: Debug with intermediate maps
+#%end
+#%flag
+#% key: c
+#% description: Clean vector lines
+#%end
+#%option
+#% key: len_plant
+#% type: double
+#% description: Maximum length of the plant [m]
+#% required: yes
+#% answer: 10000
+#%end
+#%option
+#% key: len_min
+#% type: double
+#% description: Minimum length of the plant [m]
+#% required: yes
+#% answer: 10
+#%end
+#%option
+#% key: distance
+#% type: double
+#% description: Minimum distance among plants [m]
+#% required: yes
+#% answer: 0.5
+#%end
+#%option
+#% key: p_max
+#% type: double
+#% description: Max mean power [kW]
+#% required: no
+#%end
+#%option
+#% key: output_plant
+#% type: string
+#% key_desc: name
+#% description: Name of output vector with potential power segments [kW]
+#% required: no
+#%end
+#%option
+#% key: output_point
+#% type: string
+#% key_desc: name
+#% description: Name of output vector with potential power intakes and restitution [kW]
+#% required: yes
+#%end
+#%option
+#% key: efficiency
+#% type: double
+#% description: Efficiency [-]
+#% required: yes
+#% answer: 1
+#%END
+
+from __future__ import print_function
+
+# import system libraries
+import atexit
+import os
+import sys
+
+# import grass libraries
+from grass.script import core as gcore
+#from grass.script import mapcalc
+from grass.pygrass.messages import get_msgr
+
+#from grass.pygrass.raster.buffer import Buffer
+from grass.pygrass.utils import set_path
+
+set_path('r.green', 'libhydro', '..')
+set_path('r.green', 'libgreen', os.path.join('..', '..'))
+
+from libgreen.utils import cleanup
+from libgreen.utils import dissolve_lines
+from libhydro.optimal import find_segments
+from libhydro.optimal import write_plants
+from libhydro.optimal import write_points
+
+##################################################
+# optimization problem
+# the coordinate along the river is s
+# the delta (distance betwwen intake and restitution) is delta
+# the discharge is q
+# the equation is f=[h(s,0)-h(s,delta)]*q
+# x = [s,delta]
+# s e delta are integer (the index of the vector)
+#
+
+
+if "GISBASE" not in os.environ:
+ print("You must be in GRASS GIS to run this program.")
+ sys.exit(1)
+
+
+def main(options, flags):
+ TMPRAST, TMPVECT, DEBUG = [], [], False
+ atexit.register(cleanup, rast=TMPRAST, vect=TMPVECT, debug=DEBUG)
+ elevation = options['elevation']
+ river = options['river'] # raster
+ discharge = options['discharge'] # vec
+ len_plant = float(options['len_plant'])
+ len_min = float(options['len_min'])
+ distance = float(options['distance'])
+ efficiency = float(options['efficiency'])
+ output_plant = options['output_plant']
+ output_point = options['output_point']
+ if options['p_max']:
+ p_max = float(options['p_max'])
+ else:
+ p_max = None
+ DEBUG = flags['d']
+ c = flags['c']
+ msgr = get_msgr()
+
+ # pdb.set_trace()
+
+ TMPVEC = ['river_clean']
+ if not gcore.overwrite():
+ for m in TMPVEC:
+ if gcore.find_file(m)['name']:
+ msgr.fatal(_("Temporary vector %s exists") % (m))
+
+ if c:
+ msgr.message("\Clean rivers\n")
+ dissolve_lines(river, 'river_clean')
+ river = 'river_clean'
+ # number of cell of the river
+ # range for the solution
+ msgr.message("\Loop on the category of segments\n")
+ #pdb.set_trace()
+ range_plant = (len_min, len_plant)
+ plants = find_segments(river, discharge, elevation, range_plant, distance,
+ p_max)
+ if output_plant:
+ write_plants(plants, output_plant, efficiency)
+ write_points(plants, output_point, efficiency)
+# else:
+
+if __name__ == "__main__":
+ options, flags = gcore.parser()
+ sys.exit(main(options, flags))
Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.optimal/r.green.hydro.optimal.py
___________________________________________________________________
Added: svn:mime-type
+ text/x-python
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.recommended/Makefile
===================================================================
--- grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.recommended/Makefile (rev 0)
+++ grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.recommended/Makefile 2015-04-27 06:36:52 UTC (rev 65146)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../../../..
+
+PGM = r.green.hydro.recommended
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.recommended/Makefile
___________________________________________________________________
Added: svn:mime-type
+ text/x-makefile
Added: svn:eol-style
+ native
Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.recommended/r.green.hydro.recommended.html
___________________________________________________________________
Added: svn:mime-type
+ text/html
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.recommended/r.green.hydro.recommended.py
===================================================================
--- grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.recommended/r.green.hydro.recommended.py (rev 0)
+++ grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.recommended/r.green.hydro.recommended.py 2015-04-27 06:36:52 UTC (rev 65146)
@@ -0,0 +1,196 @@
+#!/usr/bin/env python
+# -- coding: utf-8 --
+#
+############################################################################
+#
+# MODULE: r.green.hydro.reccomended
+# AUTHOR(S): Giulia Garegnani
+# PURPOSE: Calculate the optimal position of a plant along a river
+# following user reccomandations
+# COPYRIGHT: (C) 2014 by the GRASS Development Team
+#
+# This program is free software under the GNU General Public
+# License (>=v2). Read the file COPYING that comes with GRASS
+# for details.
+#
+#############################################################################
+#
+
+#%module
+#% description: Hydropower energy potential with user reccomandations
+#% keywords: raster
+#% overwrite: yes
+#%end
+#%option
+#% key: discharge
+#% type: string
+#% gisprompt: old,cell,raster
+#% key_desc: name
+#% description: Name of river discharge [l/s]
+#% required: yes
+#%end
+#%option
+#% key: river
+#% type: string
+#% gisprompt: old,vector,vector
+#% key_desc: name
+#% description: Name of vector map with interested segments of rivers
+#% required: yes
+#%end
+#%option
+#% key: dtm
+#% type: string
+#% gisprompt: old,cell,raster
+#% key_desc: name
+#% description: Name of dtm [m]
+#% required: yes
+#%end
+#%flag
+#% key: d
+#% description: Debug with intermediate maps
+#%end
+#%flag
+#% key: c
+#% description: Clean vector lines
+#%end
+#%option
+#% key: efficiency
+#% type: double
+#% key_desc: name
+#% description: efficiency [-]]
+#% required: yes
+#% answer: 0.8
+#%end
+#%option
+#% key: len_plant
+#% type: double
+#% key_desc: name
+#% description: maximum length plant [m]
+#% required: yes
+#% answer: 100
+#%end
+#%option
+#% key: len_min
+#% type: double
+#% key_desc: name
+#% description: minimum length plant [m]
+#% required: yes
+#% answer: 10
+#%end
+#%option
+#% key: distance
+#% type: double
+#% description: minimum distance among plants.
+#% required: yes
+#% answer: 0.5
+#%end
+#%option
+#% key: area
+#% type: string
+#% gisprompt: old,vector,vector
+#% key_desc: name
+#% description: areas to be excluded.
+#% required: no
+#%end
+#%option
+#% key: buff
+#% type: double
+#% description: buffer for areas to be excluded.
+#% required: no
+#%end
+#%option
+#% key: output_plant
+#% type: string
+#% key_desc: name
+#% description: Name of output vector with potential segments
+#% required: no
+#%end
+#%option
+#% key: output_point
+#% type: string
+#% key_desc: name
+#% description: Name of output vector with potential intakes and restitution
+#% required: yes
+#%END
+
+# import system libraries
+from __future__ import print_function
+import os
+import sys
+import atexit
+
+# import grass libraries
+from grass.script import core as gcore
+from grass.pygrass.utils import set_path
+from grass.pygrass.messages import get_msgr
+
+# r.green lib
+set_path('r.green', 'libhydro', '..')
+set_path('r.green', 'libgreen', os.path.join('..', '..'))
+# finally import the module in the library
+from libgreen.utils import cleanup
+
+DEBUG = False
+TMPRAST = []
+
+if "GISBASE" not in os.environ:
+ print("You must be in GRASS GIS to run this program.")
+ sys.exit(1)
+
+
+def main(opts, flgs):
+ global DEBUG, TMPRAST
+ atexit.register(cleanup, rast=TMPRAST, debug=DEBUG)
+ dtm = options['dtm']
+ river = options['river'] # raster
+ discharge = options['discharge'] # vec
+ len_plant = options['len_plant']
+ len_min = options['len_min']
+ distance = options['distance']
+ output_plant = options['output_plant']
+ output_point = options['output_point']
+ area = options['area']
+ buff = options['buff']
+ efficiency = options['efficiency']
+ DEBUG = flags['d']
+ c = flags['c']
+ msgr = get_msgr()
+
+ TMPRAST = ['new_river', 'buff_area']
+ if not gcore.overwrite():
+ for m in TMPRAST:
+ if gcore.find_file(m)['name']:
+ msgr.fatal(_("Temporary raster map %s exists") % (m))
+ #FIXME:check if it works for vectors
+
+ if area:
+ if buff:
+ gcore.run_command('v.buffer',
+ input=area,
+ output='buff_area',
+ distance=buff)
+ area = 'buff_area'
+
+ gcore.run_command('v.overlay',
+ ainput=river,
+ binput=area,
+ operator='not',
+ output='new_river')
+ river = 'new_river'
+
+ gcore.run_command('r.green.hydro.optimal',
+ flags='c',
+ discharge=discharge,
+ river=river,
+ dtm=dtm,
+ len_plant=len_plant,
+ output_plant=output_plant,
+ output_point=output_point,
+ distance=distance,
+ len_min=len_min,
+ efficiency=efficiency)
+
+if __name__ == "__main__":
+ atexit.register(cleanup)
+ options, flags = gcore.parser()
+ sys.exit(main(options, flags))
Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.recommended/r.green.hydro.recommended.py
___________________________________________________________________
Added: svn:mime-type
+ text/x-python
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.structure/Makefile
===================================================================
--- grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.structure/Makefile (rev 0)
+++ grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.structure/Makefile 2015-04-27 06:36:52 UTC (rev 65146)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../../../..
+
+PGM = r.green.hydro.structure
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.structure/Makefile
___________________________________________________________________
Added: svn:mime-type
+ text/x-makefile
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.structure/plantstructure.png
===================================================================
(Binary files differ)
Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.structure/plantstructure.png
___________________________________________________________________
Added: svn:mime-type
+ image/png
Added: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.structure/potentialplants.jpg
===================================================================
(Binary files differ)
Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.structure/potentialplants.jpg
___________________________________________________________________
Added: svn:mime-type
+ image/jpeg
Added: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.structure/r.green.hydro.structure.html
===================================================================
--- grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.structure/r.green.hydro.structure.html (rev 0)
+++ grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.structure/r.green.hydro.structure.html 2015-04-27 06:36:52 UTC (rev 65146)
@@ -0,0 +1,44 @@
+<h2>DESCRIPTION</h2>
+<em>r.green.hydro.structure</em> computes the channels and the penstock between each intake and restitution. The input has to be a vector map with the intakes and restitution of the potential plants as the one computed by r.green.hydro.optimal.
+
+
+<h2>NOTES</h2>
+The structure of the plant is composed of a channel and a penstock. The channel conveys the water from the intake to the pentock along a very small slope and the penstock conveys the water to the turbine with the highest possible head. Indeed, the power is maximized for the highest head (the discharge is the same in the river and in the structure). Thus, the channel is computed along the same quote - even if in reality a slope is necessary, it's so small that it's neglected - until the point which maximizes the head of the pentock until the turbine. The structure is computed for both sides of the river in order to determine which one produces more power.
+
+
+<h2>EXAMPLE</h2>
+These examples are based on the Pnam file which refers to the Gesso and Vermenagna Valley in Piedmont, Italy.<br><br>
+
+Here is the input vector file potentialplants with the intakes and restitution (in red) computed by r.green.hydro.optimal. The vector map with the segments of river is also visibile in blue on this picture.<br><br>
+
+<center>
+<img src="potentialplants.jpg" alt="potentialplants"><br>
+Potential intakes and restitution
+</center><br>
+
+The following command computes the channel and the pentock for each potential segment of river and for each side of the river from the input file potentialplants :<br>
+<div class="code"><pre>r.green.hydro.structure elevation=elevation plant=potentialplants output=structplants</pre></div><br><br>
+
+The result is shown in the following vector map structplants.<br><br>
+
+<center>
+<img src="plantstructure.png" alt="plantstructure"><br>
+Structure of the potential plants
+</center>
+
+
+<h2>SEE ALSO</h2>
+<em>
+<a href="r.green.hydro.optimal.html">r.green.hydro.optimal</a><br>
+<a href="r.green.hydro.technical.html">r.green.hydro.technical</a>
+</em>
+
+<h2>REFERENCE</h2>
+
+<h2>AUTHORS</h2>
+
+
+
+</div>
+</body>
+</html>
Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.structure/r.green.hydro.structure.html
___________________________________________________________________
Added: svn:mime-type
+ text/html
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.structure/r.green.hydro.structure.py
===================================================================
--- grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.structure/r.green.hydro.structure.py (rev 0)
+++ grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.structure/r.green.hydro.structure.py 2015-04-27 06:36:52 UTC (rev 65146)
@@ -0,0 +1,155 @@
+#!/usr/bin/env python
+# -- coding: utf-8 --
+#
+############################################################################
+#
+# MODULE: r.green.hydro.structure
+# AUTHOR(S):
+# PURPOSE:
+# COPYRIGHT: (C) 2014 by the GRASS Development Team
+#
+# This program is free software under the GNU General Public
+# License (>=v2). Read the file COPYING that comes with GRASS
+# for details.
+#
+#############################################################################
+#
+#%Module
+#% description: Compute channels and penstocks
+#% overwrite: yes
+#%End
+#%option G_OPT_R_ELEV
+#% required: yes
+#%end
+#%option G_OPT_V_MAP
+#% key: plant
+#% description: Name of the vector map points with the plants
+#% required: yes
+#%end
+#%option G_OPT_V_INPUT
+#% key: plant_layer
+#% description: Name of the vector map layer of plants
+#% required: no
+#% answer: 1
+#%end
+#%option
+#% key: plant_column_plant_id
+#% type: string
+#% description: Column name with the plant id
+#% required: no
+#% answer: plant_id
+#% guisection: Input columns
+#%end
+#%option
+#% key: plant_column_point_id
+#% type: string
+#% description: Column name with the point id
+#% required: no
+#% answer: cat
+#% guisection: Input columns
+#%end
+#%option
+#% key: plant_column_elevation
+#% type: string
+#% description: Column name with the elevation values [m]
+#% required: no
+#% answer: elevation
+#% guisection: Input columns
+#%end
+#%option
+#% key: plant_column_discharge
+#% type: string
+#% description: Column name with the discharge values [m3/s]
+#% required: no
+#% answer: discharge
+#% guisection: Input columns
+#%end
+#%option
+#% key: plant_column_kind
+#% type: string
+#% description: Column name (string) with the kind type of the points
+#% required: no
+#% answer: kind_label
+#% guisection: Input columns
+#%end
+
+#%option
+#% key: plant_column_kind_intake
+#% type: string
+#% description: Value contained in the column: hydro_kind that indicates the plant is an intake.
+#% required: no
+#% answer: intake
+#% guisection: Input columns
+#%end
+#%option
+#% key: plant_column_kind_turbine
+#% type: string
+#% description: Value contained in the column: hydro_kind that indicates the plant is a restitution.
+#% required: no
+#% answer: restitution
+#% guisection: Input columns
+#%end
+#%option G_OPT_V_OUTPUT
+#% key: output
+#% required: yes
+#%end
+from __future__ import print_function
+
+import os
+import atexit
+
+from grass.exceptions import ParameterError
+from grass.script.core import parser, run_command, overwrite
+from grass.pygrass.utils import set_path
+from grass.pygrass.raster import RasterRow
+from grass.pygrass.vector import VectorTopo
+
+# set python path to the shared r.green libraries
+set_path('r.green', 'libhydro', '..')
+set_path('r.green', 'libgreen', os.path.join('..', '..'))
+
+from libgreen.utils import cleanup
+from libgreen.checkparameter import check_required_columns, exception2error
+from libhydro.plant import read_plants, write_structures
+
+
+def main(opts, flgs):
+ #TMPVECT = []
+ #DEBUG = True if flgs['d'] else False
+ #atexit.register(cleanup, vect=TMPVECT, debug=DEBUG)
+ # check input maps
+ plant = [opts['plant_column_kind'], opts['plant_column_discharge'],
+ opts['plant_column_point_id'], opts['plant_column_plant_id']]
+ ovwr = overwrite()
+ #import pdb; pdb.set_trace()
+ try:
+ plnt = check_required_columns(opts['plant'], int(opts['plant_layer']),
+ plant, 'plant')
+ except ParameterError as exc:
+ exception2error(exc)
+ return
+
+ el, mset = (opts['elevation'].split('@') if '@' in opts['elevation']
+ else (opts['elevation'], ''))
+
+ elev = RasterRow(name=el, mapset=mset)
+ elev.open('r')
+ plnt.open('r')
+ #import pdb; pdb.set_trace()
+ plants, skipped = read_plants(plnt, elev=elev,
+ restitution=opts['plant_column_kind_turbine'],
+ intake=opts['plant_column_kind_intake'],
+ ckind_label=opts['plant_column_kind'],
+ cdischarge=opts['plant_column_discharge'],
+ celevation=opts['plant_column_elevation'],
+ cid_point=opts['plant_column_point_id'],
+ cid_plant=opts['plant_column_plant_id'])
+ plnt.close()
+ #import ipdb; ipdb.set_trace()
+
+ write_structures(plants, opts['output'], elev, overwrite=ovwr)
+ elev.close()
+
+
+if __name__ == "__main__":
+ main(*parser())
Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.structure/r.green.hydro.structure.py
___________________________________________________________________
Added: svn:mime-type
+ text/x-python
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.theoretical/Makefile
===================================================================
--- grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.theoretical/Makefile (rev 0)
+++ grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.theoretical/Makefile 2015-04-27 06:36:52 UTC (rev 65146)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../../../..
+
+PGM = r.green.hydro.theoretical
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.theoretical/Makefile
___________________________________________________________________
Added: svn:mime-type
+ text/x-makefile
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.theoretical/r.green.hydro.theoretical.html
===================================================================
--- grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.theoretical/r.green.hydro.theoretical.html (rev 0)
+++ grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.theoretical/r.green.hydro.theoretical.html 2015-04-27 06:36:52 UTC (rev 65146)
@@ -0,0 +1,131 @@
+<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
+<html>
+ <head>
+ <title>GRASS GIS manual: r.green.hydro.theoretical</title>
+ <meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
+ <link rel="stylesheet" href="grassdocs.css" type="text/css">
+ </head>
+ <body bgcolor="white">
+ <div id="container">
+ <a href="index.html"><img src="grass_logo.png" alt="GRASS logo"></a>
+ <hr class="header">
+ <h2>NAME</h2>
+ <em><b>r.green.hydro.theoretical</b></em> - Calculate the hydropower energy potential for each basin starting from discharge and elevation data. If existing plants are available it computes the potential installed power in the available part of the rivers.
+ <h2>KEYWORDS</h2>
+ <a href="raster.html">raster</a>
+ <h2>SYNOPSIS</h2>
+ <div id="name"><b>r.green.hydro.theoretical</b><br></div>
+ <b>r.green.hydro.theoretical --help</b><br>
+ <div id="synopsis"><b>r.green.hydro.theoretical</b> [-<b>d</b>] <b>elevation</b>=<em>name</em> <b>discharge</b>=<em>name</em> [<b>rivers</b>=<em>name</em>] [<b>lakes</b>=<em>name</em>] [<b>input</b>=<em>name</em>] <b>threshold</b>=<em>string</em> [<b>basins</b>=<em>name</em>] [<b>stream</b>=<em>name</em>] <b>output</b>=<em>name</em> [<b>segments</b>=<em>name</em>] [<b>hydro</b>=<em>name</em>] [<b>hydro_layer</b>=<em>string</em>] [<b>hydro_kind_intake</b>=<em>string</em>] [<b>hydro_kind_turbine</b>=<em>string</em>] [<b>other</b>=<em>name</em>] [<b>other_layer</b>=<em>name</em>] [<b>other_kind_intake</b>=<em>string</em>] [<b>other_kind_turbine</b>=<em>string</em>] [<b>output_segm</b>=<em>name</em>] [<b>output_point</b>=<em>name</em>] [--<b>help</b>] [--<b>verbose</b>] [--<b>quiet</b>] [--<b>ui</b>] </div>
+ <div id="flags">
+ <h3>Flags:</h3>
+ <dl>
+ <dt><b>-d</b></dt>
+ <dd>Debug with intermediate maps</dd>
+ <dt><b>--help</b></dt>
+ <dd>Print usage summary</dd>
+ <dt><b>--verbose</b></dt>
+ <dd>Verbose module output</dd>
+ <dt><b>--quiet</b></dt>
+ <dd>Quiet module output</dd>
+ <dt><b>--ui</b></dt>
+ <dd>Force launching GUI dialog</dd>
+ </dl>
+ </div>
+ <div id="parameters">
+ <h3>Parameters:</h3>
+ <h4>Required:</h4>
+ <dl>
+ <dt><b>elevation</b>=<em>name</em> <b>[required]</b></dt>
+ <dd>Name of input elevation raster map</dd>
+ <dt><b>discharge</b>=<em>name</em> <b>[required]</b></dt>
+ <dd>Name of river discharge [l/s]</dd>
+ </dl>
+ <h4>Optional:</h4>
+ <dl>
+ <dt><b>rivers</b>=<em>name</em></dt>
+ <dd>Name of river network</dd>
+ <dt><b>lakes</b>=<em>name</em></dt>
+ <dd>Name of lakes</dd>
+ </dl>
+ <h4>Discharge:</h4>
+ <dl>
+ <dt><b>input</b>=<em>name</em></dt>
+ <dd>Name of vector points map with discharge data</dd>
+ <dd>Or data source for direct OGR access</dd>
+ </dl>
+ <h4>Basin Potential:</h4>
+ <dl>
+ <dt><b>threshold</b>=<em>string</em> <b>[required]</b></dt>
+ <dd>Minimum size of exterior watershed basin</dd>
+ <dd>Default: <em>0</em></dd>
+ </dl>
+ OR if basins and streams were generated by r.watershed
+ <dl>
+ <dt><b>basins</b>=<em>name</em></dt>
+ <dd>Name of basin map obtained by r.watershed</dd>
+ <dt><b>stream</b>=<em>name</em></dt>
+ <dd>Name of stream map obtained by r.watershed</dd>
+ <dt><b>output</b>=<em>name</em> <b>[required]</b></dt>
+ <dd>Name of vector map with basin potential</dd>
+ </dl>
+ <h4>Residual Theoretical Potential:</h4>
+ <dl>
+ <dt><b>segments</b>=<em>name</em></dt>
+ <dd>Name of vector map with river segment without plants</dd>
+ <dt><b>hydro</b>=<em>name</em></dt>
+ <dd>Name of input vector map</dd>
+ <dd>Name of the vector map with points of intakes and turbines.</dd>
+ <dt><b>hydro_layer</b>=<em>string</em></dt>
+ <dd>Layer number or name</dd>
+ <dd>Name of the vector map layer of the hydro plant, with the following attributes: kind (water intake/turbine), discharge [l/s], working hours [hours], altitute[m], id point, id plant</dd>
+ <dd>Default: <em>1</em></dd>
+ <dt><b>hydro_kind_intake</b>=<em>string</em></dt>
+ <dd>Value contained in the column: hydro_kind that indicate the plant is an intake.</dd>
+ <dd>Default: <em>intake</em></dd>
+ <dt><b>hydro_kind_turbine</b>=<em>string</em></dt>
+ <dd>Value contained in the column: hydro_kind that indicate the plant is a turbine.</dd>
+ <dd>Default: <em>turbine</em></dd>
+ <dt><b>other</b>=<em>name</em></dt>
+ <dd>Name of vector map</dd>
+ <dd>Name of the vector map point of other plants such as irrigation, acqueducts, etc.</dd>
+ <dt><b>other_layer</b>=<em>name</em></dt>
+ <dd>Name of input vector map</dd>
+ <dd>Name of the vector map layer of other plants, with the following attributes: kind (water intake/turbine), discharge [m3/year], id point, id plant</dd>
+ <dd>Default: <em>1</em></dd>
+ <dt><b>other_kind_intake</b>=<em>string</em></dt>
+ <dd>Value contained in the column: other_kind that indicate the plant is an intake.</dd>
+ <dd>Default: <em>intake</em></dd>
+ <dt><b>other_kind_turbine</b>=<em>string</em></dt>
+ <dd>Value contained in the column: other_kind that indicate the plant is a turbine.</dd>
+ <dd>Default: <em>turbine</em></dd>
+ <dt><b>output_segm</b>=<em>name</em></dt>
+ <dd>Name of output vector with potential segments</dd>
+ <dt><b>output_point</b>=<em>name</em></dt>
+ <dd>Name of output vector with potential intakes and restitution</dd>
+ </dl>
+ </div>
+ <h2>SEE ALSO</h2>
+ For the hydro moduls
+ <em>
+ <a href="r.green.hydro.technical.html">r.green.hydro.technical</a>,
+ <a href="r.green.hydro.legal.html">r.green.hydro.legal</a>,
+ <a href="r.green.hydro.reccomended.html">r.green.hydro.reccomended.</a>
+ </em>
+ For the r.green modules
+ <em>
+ <a href="r.green.impact">r.green.impact.</a>,
+ </em>
+ Others modules called by r.green.hydro.theoretical are
+ <em>
+ <a href="r.green.hydro.discharge.html">r.green.hydro.discharge.</a>,
+ </em>
+ <h2>AUTHORS</h2>
+ EURAC,<br>
+ Giulia Garegnani, Pietro Zambelli
+ <hr class="header">
+ <p><a href="index.html">Main index</a> | <a href="raster.html">Raster index</a> | <a href="topics.html">Topics index</a> | <a href="keywords.html">Keywords index</a> | <a href="full_index.html">Full index</a></p>
+ <p>© 2003-2015 <a href="http://grass.osgeo.org">GRASS Development Team</a>, GRASS GIS 7.1.svn Reference Manual</p>
+ </div>
+ </body>
+</html>
\ No newline at end of file
Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.theoretical/r.green.hydro.theoretical.html
___________________________________________________________________
Added: svn:mime-type
+ text/html
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.theoretical/r.green.hydro.theoretical.py
===================================================================
--- grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.theoretical/r.green.hydro.theoretical.py (rev 0)
+++ grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.theoretical/r.green.hydro.theoretical.py 2015-04-27 06:36:52 UTC (rev 65146)
@@ -0,0 +1,296 @@
+#!/usr/bin/env python
+# -- coding: utf-8 --
+#
+############################################################################
+#
+# MODULE: r.green.hydro.potential
+# AUTHOR(S): Giulia Garegnani, Pietro Zambelli
+# PURPOSE: Calculate the theorethical hydropower energy potential for each basin and segments of river
+# COPYRIGHT: (C) 2014 by the GRASS Development Team
+#
+# This program is free software under the GNU General Public
+# License (>=v2). Read the file COPYING that comes with GRASS
+# for details.
+#
+#############################################################################
+#
+
+#%module
+#% description: Calculate the hydropower energy potential for each basin starting from discharge and elevation data. If existing plants are available it computes the potential installed power in the available part of the rivers.
+#% keywords: raster
+#% keywords: hydropower
+#% keywords: renewable energy
+#%end
+#%option G_OPT_R_ELEV
+#% required: yes
+#%end
+#%option
+#% key: discharge
+#% type: string
+#% gisprompt: old,cell,raster
+#% key_desc: name
+#% description: Name of river discharge [m3/s]
+#% required: yes
+#%end
+#%option
+#% key: rivers
+#% type: string
+#% gisprompt: old,vector,vector
+#% key_desc: name
+#% description: Name of river network
+#% required: no
+#%end
+#%option
+#% key: lakes
+#% type: string
+#% gisprompt: old,vector,vector
+#% key_desc: name
+#% description: Name of lakes
+#% required: no
+#%end
+#%option
+#% key: threshold
+#% type: string
+#% description: Minimum size of exterior watershed basin
+#% required: yes
+#% answer: 0
+#% guisection: Basin Potential
+#%end
+#%option
+#% key: basins
+#% type: string
+#% gisprompt: old,cell,raster
+#% key_desc: name
+#% description: Name of basin map obtained by r.watershed
+#% required: no
+#% guisection: Basin Potential
+#%end
+#%option
+#% key: stream
+#% type: string
+#% gisprompt: old,cell,raster
+#% key_desc: name
+#% description: Name of stream map obtained by r.watershed
+#% required: no
+#% guisection: Basin Potential
+#%end
+#TODO: add flags
+#%flag
+#% key: d
+#% description: Debug with intermediate maps
+#%end
+#%option
+#% key: output
+#% type: string
+#% key_desc: name
+#% description: Name of vector map with basin potential
+#% required: yes
+#% guisection: Basin Potential
+#%END
+
+from __future__ import print_function
+
+# import system libraries
+import os
+import sys
+import atexit
+import pdb
+
+# import grass libraries
+from grass.script import core as gcore
+from grass.pygrass.messages import get_msgr
+from grass.pygrass.utils import set_path
+
+set_path('r.green', 'libhydro', '..')
+set_path('r.green', 'libgreen', os.path.join('..', '..'))
+# finally import the module in the library
+from libgreen.utils import cleanup
+from libhydro import basin
+from libgreen.utils import check_overlay_rr
+#from libgreen.utils import check_overlay_rv
+from libgreen.utils import raster2numpy
+from libgreen.utils import remove_pixel_from_raster
+
+
+if "GISBASE" not in os.environ:
+ print("You must be in GRASS GIS to run this program.")
+ sys.exit(1)
+
+
+def main(options, flags):
+ #############################################################
+ # inizialitation
+ #############################################################
+ TMPRAST, TMPVECT, DEBUG = [], [], False
+ atexit.register(cleanup, rast=TMPRAST, vect=TMPVECT, debug=DEBUG)
+ #TOD_add the possibilities to have q_points
+ # required
+ discharge = options['discharge']
+ dtm = options['elevation']
+ # basin potential
+ basins = options['basins']
+ stream = options['stream'] # raster
+ rivers = options['rivers'] # vec
+ lakes = options['lakes'] # vec
+ E = options['output']
+ threshold = options['threshold']
+ # existing plants
+# segments = options['segments']
+# output_segm = options['output_segm']
+# output_point = options['output_point']
+# hydro = options['hydro']
+# hydro_layer = options['hydro_layer']
+# hydro_kind_intake = options['hydro_kind_intake']
+# hydro_kind_turbine = options['hydro_kind_turbine']
+# other = options['other']
+# other_kind_intake = options['other_kind_intake']
+# other_kind_turbine = options['other_kind_turbine']
+
+ # optional
+ DEBUG = flags['d']
+
+ msgr = get_msgr()
+ # info = gcore.parse_command('g.region', flags='m')
+ # print info
+ #############################################################
+ # check temporary raster
+ #############################################################
+ TMPRAST = ['dtm_mean', 'neighbors',
+ 'closure', 'down', 'lakes', 'stream_thin',
+ 'debug', 'vec_rast']
+ TMPVECT = ['tmptmp', 'segments']
+
+ if not gcore.overwrite():
+ for m in TMPRAST:
+ if gcore.find_file(m)['name']:
+ msgr.fatal(_("Temporary raster map %s exists") % (m))
+
+ if rivers:
+ # cp the vector in the current mapset in order to clean it
+
+ to_copy = '%s,river_new' % rivers
+ gcore.run_command('g.copy', vector=to_copy)
+ rivers = 'river_new'
+ TMPVECT.append('river_new')
+ gcore.run_command('v.build', map=rivers)
+ pdb.set_trace()
+ dtm_corr, tmp_v, tmp_r = basin.dtm_corr(dtm, rivers, lakes)
+ TMPRAST = TMPRAST + tmp_r
+ TMPVECT = TMPVECT + tmp_v
+ basins, stream = basin.check_compute_basin_stream(basins, stream,
+ dtm_corr, threshold)
+ else:
+ basins, stream = basin.check_compute_basin_stream(basins, stream,
+ dtm, threshold)
+
+ perc_overlay = check_overlay_rr(discharge, stream)
+ #pdb.set_trace()
+ if float(perc_overlay) < 90:
+ warn = ("Discharge map doesn't overlay all the stream map."
+ "It covers only the %s %% of rivers") % (perc_overlay)
+ msgr.warning(warn)
+
+ msgr.message("\Init basins\n")
+ #pdb.set_trace()
+ #############################################################
+ # core
+ #############################################################
+ basins_tot, inputs = basin.init_basins(basins)
+
+ msgr.message("\nCompute mean elevation\n")
+
+ gcore.run_command('r.stats.zonal',
+ base=basins,
+ cover=dtm, flags='r',
+ output='dtm_mean',
+ method='average')
+
+ msgr.message("\nBuild the basin network\n")
+ #############################################################
+ # build_network(stream, dtm, basins_tot) build relationship among basins
+ # Identify the closure point with the r.neigbours module
+ # if the difference betwwen the stream and the neighbours
+ # is negative it means a new subsanin starts
+ #############################################################
+ basin.build_network(stream, dtm, basins_tot)
+ stream_n = raster2numpy(stream)
+ discharge_n = raster2numpy(discharge)
+ #pdb.set_trace()
+ basin.fill_basins(inputs, basins_tot, basins, dtm, discharge_n, stream_n)
+
+ ###################################################################
+ # check if lakes and delate stream in lakes optional
+ ###################################################################
+
+ if lakes:
+ remove_pixel_from_raster(lakes, stream)
+
+ ####################################################################
+ # write results
+ ####################################################################
+ # if not rivers or debug I write the result in the new vector stream
+ msgr.message("\nWrite results\n")
+
+ basin.write_results2newvec(stream, E, basins_tot, inputs)
+
+ ####################################################################
+ # try tobuilt a vector shp of rivers coerent with vector input
+ #TODO: and the new branches created by r.watershed, check how to add to
+ # the existing vector river
+ ####################################################################
+#==============================================================================
+# if rivers:
+# if DEBUG:
+# basin.write_results2newvec(stream, 'debug', basins_tot, inputs)
+# basin.add_results2shp(basins, rivers, basins_tot, E, inputs)
+# else:
+# basin.add_results2shp(basins, rivers, basins_tot, E, inputs)
+#==============================================================================
+
+ ####################################################################
+ # if there are plants, it deletes the piece of river
+ # pdb.set_trace()
+# if (segments and hydro):
+# warn = (("%s map could have been generated by %s."
+# "%s is used for next computations")
+# % (segments, hydro, segments))
+# msgr.warning(warn)
+#
+# if hydro and not(segments):
+# msgr.message("\nr.green.hydro.delplants\n")
+# gcore.run_command('r.green.hydro.delplants',
+# hydro=hydro,
+# hydro_layer=hydro_layer,
+# river=E,
+# plants='plant',
+# output='segments',
+# hydro_kind_intake=hydro_kind_intake,
+# hydro_kind_turbine=hydro_kind_turbine,
+# elevation=dtm,
+# other=other,
+# other_kind_intake=other_kind_intake,
+# other_kind_turbine=other_kind_turbine)
+# segments = 'segments'
+#
+# #FIXME: with temporary plants it don't work
+# if segments:
+# perc_overlay = check_overlay_rv(discharge, segments)
+# if float(perc_overlay) < 90:
+# warn = ("Discharge map doesn't overlay all the segments map."
+# "It covers only the %s %% of rivers") % (perc_overlay)
+# msgr.warning(warn)
+# gcore.run_command('r.green.hydro.optimal',
+# flags='c',
+# discharge=discharge,
+# river=segments,
+# dtm=dtm,
+# len_plant='10000000',
+# output_plant=output_segm,
+# output_point=output_point,
+# distance='0.5',
+# len_min='0.5')
+
+
+if __name__ == "__main__":
+ options, flags = gcore.parser()
+ sys.exit(main(options, flags))
Property changes on: grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.theoretical/r.green.hydro.theoretical.py
___________________________________________________________________
Added: svn:mime-type
+ text/x-python
Added: svn:eol-style
+ native
More information about the grass-commit
mailing list