[GRASS-SVN] r60813 - sandbox/annakrat/r3.flow
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Jun 12 20:52:47 PDT 2014
Author: annakrat
Date: 2014-06-12 20:52:47 -0700 (Thu, 12 Jun 2014)
New Revision: 60813
Modified:
sandbox/annakrat/r3.flow/r3.flow.py
sandbox/annakrat/r3.flow/rast3d_functions.py
Log:
r3.flow: added more input options; simple flow accumulation implemented
Modified: sandbox/annakrat/r3.flow/r3.flow.py
===================================================================
--- sandbox/annakrat/r3.flow/r3.flow.py 2014-06-12 22:26:32 UTC (rev 60812)
+++ sandbox/annakrat/r3.flow/r3.flow.py 2014-06-13 03:52:47 UTC (rev 60813)
@@ -36,12 +36,20 @@
#%end
#%option G_OPT_V_INPUT
#% key: seed_points
-#% description:
+#% required: no
+#% description: If no map is provided, flow lines are generated from each cell of the input 3D raster
#% label: Name of vector map with points from which flow lines are generated
#%end
#%option G_OPT_V_OUTPUT
+#% key: flowline
+#% required: no
#% description: Name of vector map of flow lines
#%end
+#%option G_OPT_R3_OUTPUT
+#% key: flowaccumulation
+#% required: no
+#% description: Name of the output flow accumulation 3D raster map
+#%end
#%option
#% key: unit
#% type: string
@@ -69,12 +77,20 @@
#% description: Maximum number of steps
#% gisprompt: Integration
#%end
+#%option
+#% key: skip
+#% type: integer
+#% required: no
+#% multiple: yes
+#% description: Number of cells between flow lines in x, y and z direction
+#%end
#%flag
#% key: u
#% description: Compute upstream flowlines instead of default downstream flowlines
#%end
import math
+import numpy as np
from grass.script import core as gcore
from grass.script import array as garray
@@ -83,11 +99,20 @@
from grass.pygrass.vector.geometry import Line, Point
from integrate import integrate_next, get_time_step
-from rast3d_functions import get_velocity, norm, compute_gradient
+from rast3d_functions import get_velocity, norm, compute_gradient, rast3d_location2coord
-EPSILON = 1e-7
+EPSILON = 1e-8
+class Seed:
+ def __init__(self, x, y, z, flowline, flowaccum):
+ self.x = x
+ self.y = y
+ self.z = z
+ self.flowline = flowline
+ self.flowaccum = flowaccum
+
+
def main():
options, flags = gcore.parser()
@@ -106,12 +131,15 @@
elif options['input']:
array = garray.array3d()
array.read(mapname=options['input'])
- velocity = compute_gradient(array)
+ velocity = compute_gradient(array)
map_info = grast3d.raster3d_info(options['input'])
else:
# this should be handled by parser (in the future)
gcore.fatal("Either input or vector_field options must be specified.")
+ if not (options['flowaccumulation'] or options['flowline']):
+ gcore.fatal("At least one of output options flowline or flowaccumulation must be specified.")
+
# integrating forward means upstream, backward means downstream
if flags['u']:
direction = 1
@@ -130,26 +158,69 @@
# cell size is the diagonal
cell_size = math.sqrt(map_info['nsres'] ** 2 + map_info['ewres'] ** 2 + map_info['tbres'] ** 2)
- # read seed points
+ if options['seed_points'] and options['skip']:
+ gcore.warning(_("Option skip is ignored because seed point map was provided."))
+
+ # read or create seed points
seeds = []
- inp_seeds = VectorTopo(options['seed_points'])
- inp_seeds.open()
- for point in inp_seeds:
- point.is2D = False # just a workaround
- seeds.append((point.x, point.y, point.z))
- inp_seeds.close()
+ if options['flowaccumulation'] or (options['flowline'] and not options['seed_points']):
+ if options['skip']:
+ try:
+ skipx, skipy, skipz = [int(each) for each in options['skip'].split(',')]
+ except ValueError:
+ gcore.fatal(_("Please provide 3 values for skip option"))
+ else:
+ skipx = max(1, int(map_info['cols'] / 10))
+ skipy = max(1, int(map_info['rows'] / 10))
+ skipz = max(1, int(map_info['depths'] / 10))
+ for r in range(map_info['rows'], 0, -1):
+ for c in range(0, map_info['cols']):
+ for d in range(0, map_info['depths']):
+ x = c * map_info['ewres'] + map_info['ewres'] / 2
+ y = r * map_info['nsres'] + map_info['nsres'] / 2
+ z = d * map_info['tbres'] + map_info['tbres'] / 2
+
+ seed = Seed(x, y, z, False, False)
+ if options['flowaccumulation']:
+ seed.flowaccum = True
+ if options['flowline'] and c % skipx == 0 and r % skipy == 0 and d % skipz == 0:
+ seed.flowline = True
+ if seed.flowaccum or seed.flowline:
+ seeds.append(seed)
+
+ else:
+ # read seed points
+ inp_seeds = VectorTopo(options['seed_points'])
+ inp_seeds.open()
+ for point in inp_seeds:
+ point.is2D = False # just a workaround
+ seeds.append(Seed(point.x, point.y, point.z, flowaccum=False, flowline=True))
+ inp_seeds.close()
+
# create new vector map to store flow lines
- flowline_vector = VectorTopo(options['output'])
- if flowline_vector.exist() and not gcore.overwrite:
- gcore.fatal(_("Vector map <{flowlines}> already exists.").format(flowlines=flowline_vector))
- flowline_vector.open('w')
+ if options['flowline']:
+ flowline_vector = VectorTopo(options['flowline'])
+ if flowline_vector.exist() and not gcore.overwrite:
+ gcore.fatal(_("Vector map <{flowline}> already exists.").format(flowline=flowline_vector))
+ flowline_vector.open('w')
+ # create the flowaccumulation 3D raster
+ if options['flowaccumulation']:
+ flowacc = garray.array3d()
+ for depth in range(map_info['depths']):
+ flowacc[depth] = np.zeros((map_info['rows'], map_info['cols']))
+
# go through all points and integrate flowline for each
+ seed_count = 0
for seed in seeds:
- line = Line([Point(seed[0], seed[1], seed[2])])
- point = seed
count = 0
+ gcore.percent(seed_count, len(seeds), 1)
+ seed_count += 1
+ point = seed.x, seed.y, seed.z
+ if seed.flowline:
+ line = Line([Point(point[0], point[1], point[2])])
+ last_coords = (None, None, None)
while count <= limit:
velocity_vector = get_velocity(map_info, velocity, point)
if velocity_vector is None:
@@ -163,17 +234,27 @@
point = integrate_next(map_info, velocity, point, delta_T * direction)
if point is None:
break
- line.append(Point(point[0], point[1], point[2]))
+ if seed.flowline:
+ line.append(Point(point[0], point[1], point[2]))
+ if seed.flowaccum:
+ col, row, depth = rast3d_location2coord(map_info, point[1], point[0], point[2])
+ if last_coords != (col, row, depth):
+ flowacc[depth][row][col] += 1
+ last_coords = (col, row, depth)
count += 1
- if len(line) > 1:
+ if seed.flowline and len(line) > 1:
# workaround for pygrass bug
line.is2D = False
flowline_vector.write(line)
- gcore.message(_("Flowline ended after %s steps.") % (count - 1), flag='v')
+ gcore.verbose(_("Flowline ended after %s steps") % (count - 1))
- flowline_vector.close()
+ if options['flowline']:
+ flowline_vector.close()
+ if options['flowaccumulation']:
+ flowacc.write(mapname=options['flowaccumulation'])
+
if __name__ == '__main__':
main()
Modified: sandbox/annakrat/r3.flow/rast3d_functions.py
===================================================================
--- sandbox/annakrat/r3.flow/rast3d_functions.py 2014-06-12 22:26:32 UTC (rev 60812)
+++ sandbox/annakrat/r3.flow/rast3d_functions.py 2014-06-13 03:52:47 UTC (rev 60813)
@@ -34,9 +34,10 @@
def rast3d_location2coord(rast3d_region, north, east, top):
"""!Returns column, row and depth of a cell where the input
coordinates (north, east, top) lie."""
- col = int((east - rast3d_region['west']) / rast3d_region['ewres'])
- row = int((rast3d_region['north'] - north) / rast3d_region['nsres'])
- depth = int((top - rast3d_region['bottom']) / rast3d_region['tbres'])
+ # math.floor to get because int(-0.8) gives 0
+ col = int(math.floor((east - rast3d_region['west']) / rast3d_region['ewres']))
+ row = int(math.floor((rast3d_region['north'] - north) / rast3d_region['nsres']))
+ depth = int(math.floor((top - rast3d_region['bottom']) / rast3d_region['tbres']))
return col, row, depth
More information about the grass-commit
mailing list