[GRASS-SVN] r61143 - sandbox/annakrat/r3.flow
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Jul 3 15:48:19 PDT 2014
Author: annakrat
Date: 2014-07-03 15:48:19 -0700 (Thu, 03 Jul 2014)
New Revision: 61143
Modified:
sandbox/annakrat/r3.flow/r3.flow.py
Log:
r3.flow: reorganize main to avoid creating array of seed points
Modified: sandbox/annakrat/r3.flow/r3.flow.py
===================================================================
--- sandbox/annakrat/r3.flow/r3.flow.py 2014-07-03 15:02:22 UTC (rev 61142)
+++ sandbox/annakrat/r3.flow/r3.flow.py 2014-07-03 22:48:19 UTC (rev 61143)
@@ -94,7 +94,6 @@
from grass.script import core as gcore
from grass.script import array as garray
-from grass.script import raster3d as grast3d
from grass.pygrass.vector import VectorTopo
from grass.pygrass.vector.geometry import Line, Point
@@ -124,10 +123,69 @@
self.neighbors_pos = None
+class Integrate:
+ def __init__(self):
+ self.direction = None
+ self.unit = None
+ self.step = None
+ self.cell_size = None
+ self.limit = None
+
+
+def compute_flowline(map_info, seed, velocity_obj, integration, flowline_vector, flowacc):
+ count = 0
+ point = seed.x, seed.y, seed.z
+ if seed.flowline:
+ line = Line([Point(seed.x, seed.y, seed.z)])
+ last_coords = (None, None, None)
+ new_point = [None, None, None]
+ while count <= integration.limit:
+ velocity_vector = get_velocity(map_info, velocity_obj, point)
+ if velocity_vector is None:
+ break # outside region
+ velocity_norm = norm(velocity_vector)
+
+ if velocity_norm <= EPSILON:
+ break # zero velocity means end of propagation
+ # convert to time
+ delta_T = get_time_step(integration.unit, integration.step, velocity_norm, integration.cell_size)
+ # decide which integration method to choose
+ if not RK45:
+ ret = rk4_integrate_next(map_info, velocity_obj, point, new_point,
+ delta_T * integration.direction)
+ else:
+ min_step = get_time_step('cell', MIN_STEP, velocity_norm, integration.cell_size)
+ max_step = get_time_step('cell', MAX_STEP, velocity_norm, integration.cell_size)
+ ret = rk45_integrate_next(map_info, velocity_obj, point, new_point,
+ delta_T * integration.direction, min_step, max_step)
+ if not ret:
+ break
+
+ if seed.flowline:
+ line.append(Point(new_point[0], new_point[1], new_point[2]))
+ if seed.flowaccum:
+ col, row, depth = rast3d_location2coord(map_info, new_point[1], new_point[0], new_point[2])
+ if last_coords != (col, row, depth):
+ flowacc[depth][row][col] += 1
+ if last_coords[0] is not None and \
+ sum([abs(coord[0] - coord[1]) for coord in zip(last_coords, (col, row, depth))]) > 1:
+ additional = traverse(point, new_point, map_info)
+ for p in additional:
+ flowacc[p[2]][p[1]][p[0]] += 1
+ last_coords = (col, row, depth)
+ point = new_point
+ count += 1
+
+ if seed.flowline and len(line) > 1:
+ flowline_vector.write(line)
+ gcore.verbose(_("Flowline ended after %s steps") % (count - 1))
+
+
def main():
options, flags = gcore.parser()
velocity_obj = Velocity()
+ integration = Integrate()
if options['vector_field']:
try:
v_x, v_y, v_z = options['vector_field'].split(',')
@@ -154,132 +212,109 @@
gcore.fatal("At least one of output options flowline or flowaccumulation must be specified.")
map_info = get_region_info()
+
+ 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))
+
+ # initialize integration variables
+
# integrating forward means upstream, backward means downstream
if flags['u']:
- direction = 1
+ integration.direction = 1
else:
- direction = -1
+ integration.direction = -1
- limit = int(options['limit'])
+ integration.limit = int(options['limit'])
# see which units to use
step = options['step']
if step:
- step = float(step)
- unit = options['unit']
+ integration.step = float(step)
+ integration.unit = options['unit']
else:
- unit = 'cell'
- step = 0.25
+ integration.unit = 'cell'
+ integration.step = 0.25
# cell size is the diagonal
- cell_size = math.sqrt(map_info['nsres'] ** 2 + map_info['ewres'] ** 2 + map_info['tbres'] ** 2)
+ integration.cell_size = math.sqrt(map_info['nsres'] ** 2 + map_info['ewres'] ** 2 +
+ map_info['tbres'] ** 2)
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 = []
+ # create new vector map to store flow lines
+ 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(mode='w', with_z=True)
+
+ # 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']))
+
+ # try open seed points map and get number of points
+ if options['seed_points']:
+ inp_seeds = VectorTopo(options['seed_points'])
+ inp_seeds.open()
+ if not inp_seeds.is_3D():
+ inp_seeds.close()
+ gcore.fatal(_("Input vector map of seed points must be 3D."))
+ num_points = inp_seeds.num_primitive_of('point')
+ if num_points <= 0:
+ inp_seeds.close()
+ gcore.fatal(_("No points found in seed vector map."))
+
+ # compute number of total seeds to show percent
+ num_seeds = 0
+ if options['seed_points']:
+ num_seeds += num_points
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"))
+ if options['flowaccumulation']:
+ num_seeds += map_info['cols'] * map_info['rows'] * map_info['depths']
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))
+ num_seeds += (math.ceil(map_info['cols'] / float(skipx)) *
+ math.ceil(map_info['rows'] / float(skipy)) *
+ math.ceil(map_info['depths'] / float(skipz)))
+ gcore.info(_("Flow lines for {count} seed points will be computed.").format(count=num_seeds))
+ # compute flowlines from vector seed points
+ seed_count = 0
+ if options['seed_points']:
+ for point in inp_seeds:
+ gcore.percent(seed_count, num_seeds, 1)
+ seed = Seed(point.x, point.y, point.z, flowaccum=False, flowline=True)
+ compute_flowline(map_info, seed, velocity_obj, integration, flowline_vector, flowacc)
+ seed_count += 1
+
+ inp_seeds.close()
+
+ # compute flowlines from seeds created on grid
+ if options['flowaccumulation'] or (options['flowline'] and not options['seed_points']):
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 = map_info['west'] + c * map_info['ewres'] + map_info['ewres'] / 2
y = map_info['south'] + r * map_info['nsres'] - map_info['nsres'] / 2
z = map_info['bottom'] + 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)
+ gcore.percent(seed_count, num_seeds, 1)
+ compute_flowline(map_info, seed, velocity_obj, integration, flowline_vector, flowacc)
+ seed_count += 1
- else:
- # read seed points
- inp_seeds = VectorTopo(options['seed_points'])
- inp_seeds.open()
- if not inp_seeds.is_3D():
- gcore.fatal(_("Input vector map of seed points must be 3D."))
- for point in inp_seeds:
- seeds.append(Seed(point.x, point.y, point.z, flowaccum=False, flowline=True))
- inp_seeds.close()
-
- # create new vector map to store flow lines
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(mode='w', with_z=True)
- # 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
- gcore.info(_("Flow lines for {count} seed points will be computed.").format(count=len(seeds)))
- for seed in seeds:
- 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)
- new_point = [None, None, None]
- while count <= limit:
- velocity_vector = get_velocity(map_info, velocity_obj, point)
- if velocity_vector is None:
- break # outside region
- velocity_norm = norm(velocity_vector)
-
- if velocity_norm <= EPSILON:
- break # zero velocity means end of propagation
- # convert to time
- delta_T = get_time_step(unit, step, velocity_norm, cell_size)
- # decide which integration method to choose
- if not RK45:
- ret = rk4_integrate_next(map_info, velocity_obj, point, new_point,
- delta_T * direction)
- else:
- min_step = get_time_step('cell', MIN_STEP, velocity_norm, cell_size)
- max_step = get_time_step('cell', MAX_STEP, velocity_norm, cell_size)
- ret = rk45_integrate_next(map_info, velocity_obj, point, new_point,
- delta_T * direction, min_step, max_step)
- if not ret:
- break
-
- if seed.flowline:
- line.append(Point(new_point[0], new_point[1], new_point[2]))
- if seed.flowaccum:
- col, row, depth = rast3d_location2coord(map_info, new_point[1], new_point[0], new_point[2])
- if last_coords != (col, row, depth):
- flowacc[depth][row][col] += 1
- if last_coords[0] is not None and \
- sum([abs(coord[0] - coord[1]) for coord in zip(last_coords, (col, row, depth))]) > 1:
- additional = traverse(point, new_point, map_info)
- for p in additional:
- flowacc[p[2]][p[1]][p[0]] += 1
- last_coords = (col, row, depth)
- point = new_point
- count += 1
-
- if seed.flowline and len(line) > 1:
- # workaround for pygrass bug
- flowline_vector.write(line)
- gcore.verbose(_("Flowline ended after %s steps") % (count - 1))
-
- if options['flowline']:
flowline_vector.close()
if options['flowaccumulation']:
flowacc.write(mapname=options['flowaccumulation'])
More information about the grass-commit
mailing list