[GRASS-SVN] r65958 - in grass-addons/grass7/raster/r.green/r.green.hydro: libhydro r.green.hydro.structure
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Aug 18 03:03:50 PDT 2015
Author: zarch
Date: 2015-08-18 03:03:50 -0700 (Tue, 18 Aug 2015)
New Revision: 65958
Modified:
grass-addons/grass7/raster/r.green/r.green.hydro/libhydro/plant.py
grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.structure/r.green.hydro.structure.py
Log:
r.green: Improve performances for the identification of power hydro structures and add ndigits and resolution parameters
Modified: grass-addons/grass7/raster/r.green/r.green.hydro/libhydro/plant.py
===================================================================
--- grass-addons/grass7/raster/r.green/r.green.hydro/libhydro/plant.py 2015-08-18 06:16:29 UTC (rev 65957)
+++ grass-addons/grass7/raster/r.green/r.green.hydro/libhydro/plant.py 2015-08-18 10:03:50 UTC (rev 65958)
@@ -9,6 +9,7 @@
from grass.pygrass.vector import VectorTopo
from grass.pygrass.vector.geometry import Line
from grass.pygrass.vector.table import Link
+from grass.pygrass import utils
COLS = [(u'cat', 'INTEGER PRIMARY KEY'),
@@ -33,6 +34,31 @@
['intake', 'conduct', 'penstock', 'side'])
+def closest(number, ndigits=0, resolution=None):
+ """Round a number defining the number of precision decimal digits and
+ the resolution.
+
+ Examples
+ --------
+
+ >>> closest(103.66778,)
+ 104.0
+ >>> closest(103.66778, ndigits=2)
+ 103.67
+ >>> closest(103.66778, ndigits=2, resolution=5)
+ 105.0
+ >>> closest(103.66778, ndigits=2, resolution=0.5)
+ 103.5
+ >>> closest(103.66778, ndigits=2, resolution=0.25)
+ 103.75
+ """
+ num = round(number, ndigits)
+ return (num if resolution is None else
+ round((num // resolution * resolution +
+ round((num % resolution) / float(resolution), 0) *
+ resolution), ndigits))
+
+
def isinverted(line, elev, region):
el0, el1 = elev.get_value(line[0]), elev.get_value(line[-1])
return False if el0 > el1 else True
@@ -125,16 +151,17 @@
out.table.conn.commit()
-def write_structures(plants, output, elev, stream=None, overwrite=False):
+def write_structures(plants, output, elev, stream=None,
+ ndigits=0, resolution=None, contour='',
+ 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))
+ 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,
@@ -149,6 +176,7 @@
(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')
@@ -156,17 +184,52 @@
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)
+ # check if contour vector map is provide by the user
+ if contour:
+ cname, cmset = (contour.split('@') if '@' in contour
+ else (contour, ''))
+ # check if the map already exist
+ if bool(utils.get_mapset_vector(cname, cmset)) and overwrite:
+ compute_contour = True
+ remove = False
+ else:
+ # create a random name
+ contour = 'tmpvect%02d' % random.randint(0, 99)
+ compute_contour = True
+ remove = True
+ if compute_contour:
+ # compute the levels of the contour lines map
+ levels = []
+ for p in plants.values():
+ for itk in p.intakes:
+ levels.append(closest(itk.elevation, ndigits=ndigits,
+ resolution=resolution))
+ levels = sorted(set(levels))
+
+ # generate the contur line that pass to the point
+ r.contour(input='%s@%s' % (elev.name, elev.mapset),
+ output=contour, step=0, levels=levels, overwrite=True)
+
+ # open the contur lines
+ with VectorTopo(contour, mode='r') as cnt:
+ for plant in plants.values():
+ print(plant.id)
+ for options in plant.structures(elev, stream=stream,
+ ndigits=ndigits,
+ resolution=resolution,
+ contour=cnt):
+ for hydro in options:
+ print('writing: ', hydro.intake)
+ write_hydrostruct(out, hydro, plant)
+
+ if remove:
+ cnt.remove()
+
+
class AbstractPoint(object):
def __init__(self, id_point, point, elevation, id_plants=None,
id_stream=None, discharge=1.):
@@ -356,7 +419,8 @@
ids.append(l_id)
return lines, ids
- def structures(self, elev, stream=None):
+ def structures(self, elev, stream=None,
+ ndigits=0, resolution=None, contour=None):
"""Return a tuple with lines structres options of a hypotetical plant.
::
@@ -435,37 +499,52 @@
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)
+ if contour is None:
+ levels = sorted(set([closest(itk.elevation,
+ ndigits=ndigits, resolution=resolution)
+ for itk in self.intakes]))
+
# 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
+ contour_tmp = 'tmpvect%04d' % random.randint(1000, 9999)
r.contour(input='%s@%s' % (elev.name, elev.mapset),
- output=tmpvect, step=0, levels=itk.elevation,
+ output=contour_tmp, step=0, levels=levels,
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 = VectorTopo(contour_tmp)
+ cnt.open()
+ else:
+ cnt = contour
+
+ for itk in self.intakes:
+ # find the closest contur line
+ contur_res = cnt.find['by_point'].geo(self.restitution.point,
+ maxdist=100000.0)
+
+ # 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."
+ "Therefore %d will be used.")
+ print(msg % (contur_itk.id, contur_res.id, cnt.name,
+ contur_itk.id))
+
+ # check contour
+ contur = not_overlaped(contur_itk)
+ contur = sort_by_west2east(contur)
+ result.append(get_all_structs(contur, itk, self.restitution))
+
+ # remove temporary vector map
+ if contour is None:
+ cnt.close()
cnt.remove()
return result
Modified: 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 2015-08-18 06:16:29 UTC (rev 65957)
+++ grass-addons/grass7/raster/r.green/r.green.hydro/r.green.hydro.structure/r.green.hydro.structure.py 2015-08-18 10:03:50 UTC (rev 65958)
@@ -89,7 +89,28 @@
#% answer: pot_power
#% guisection: Input columns
#%end
+#%option
+#% key: ndigits
+#% type: integer
+#% description: Number of digits to use for the elevation in the contour line vector map
+#% required: no
+#% answer: 0
+#% guisection: Contour
+#%end
+#%option
+#% key: resolution
+#% type: double
+#% description: Resolution use for the contour line vector map, if 0.25 approximate 703.31 tp 703.25
+#% required: no
+#% guisection: Contour
+#%end
#%option G_OPT_V_OUTPUT
+#% key: contour
+#% description: Name of the contour line vector map
+#% required: no
+#% guisection: Contour
+#%end
+#%option G_OPT_V_OUTPUT
#% key: output_point
#% label: Name of output vector map with potential intakes and restitution
#% required: no
@@ -172,7 +193,11 @@
plnt.close()
- write_structures(plants, opts['output_struct'], elev, overwrite=ovwr)
+ # contour options
+ resolution = float(opts['resolution']) if opts['resolution'] else None
+ write_structures(plants, opts['output_struct'], elev,
+ ndigits=int(opts['ndigits']), resolution=resolution,
+ contour=opts['contour'], overwrite=ovwr)
elev.close()
More information about the grass-commit
mailing list