[GRASS-SVN] r60673 - sandbox/annakrat/r3.flow
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Jun 2 11:40:18 PDT 2014
Author: annakrat
Date: 2014-06-02 11:40:18 -0700 (Mon, 02 Jun 2014)
New Revision: 60673
Modified:
sandbox/annakrat/r3.flow/integrate.py
sandbox/annakrat/r3.flow/r3.flow.py
sandbox/annakrat/r3.flow/rast3d_functions.py
sandbox/annakrat/r3.flow/test_interpolation.py
Log:
r3.flow: first working version, support for different integration units
Modified: sandbox/annakrat/r3.flow/integrate.py
===================================================================
--- sandbox/annakrat/r3.flow/integrate.py 2014-06-02 17:03:50 UTC (rev 60672)
+++ sandbox/annakrat/r3.flow/integrate.py 2014-06-02 18:40:18 UTC (rev 60673)
@@ -15,36 +15,44 @@
k4 = [None] * 3
x, y, z = point[0], point[1], point[2]
# first
- velocity = get_velocity(rast3d_region, arrays[0], arrays[1], arrays[2], x, y, z)
+ velocity = get_velocity(rast3d_region, arrays, point)
if velocity is None:
return None
for i in range(3):
k1[i] = delta_t * velocity[i]
# second
- velocity = get_velocity(rast3d_region, arrays[0], arrays[1], arrays[2],
- x + k1[0] / 2, y + k1[1] / 2, z + k1[2] / 2)
+ velocity = get_velocity(rast3d_region, arrays,
+ (x + k1[0] / 2, y + k1[1] / 2, z + k1[2] / 2))
if velocity is None:
return None
for i in range(3):
k2[i] = delta_t * velocity[i]
# third
- velocity = get_velocity(rast3d_region, arrays[0], arrays[1], arrays[2],
- x + k2[0] / 2, y + k2[1] / 2, z + k2[2] / 2)
+ velocity = get_velocity(rast3d_region, arrays,
+ (x + k2[0] / 2, y + k2[1] / 2, z + k2[2] / 2))
if velocity is None:
return None
for i in range(3):
k3[i] = delta_t * velocity[i]
# fourth
- velocity = get_velocity(rast3d_region, arrays[0], arrays[1], arrays[2],
- x + k3[0], y + k3[1], z + k3[2])
+ velocity = get_velocity(rast3d_region, arrays, (x + k3[0], y + k3[1], z + k3[2]))
if velocity is None:
return None
for i in range(3):
k4[i] = delta_t * velocity[i]
- print k1, k2, k3, k4
+
# next point
next_point = [None] * 3
for i in range(3):
next_point[i] = point[i] + k1[i] / 6 + k2[i] / 3 + k3[i] / 3 + k4[i] / 6
return next_point
+
+
+def get_time_step(unit, step, velocity, cell_size):
+ if unit == 'time':
+ return step
+ elif unit == 'length':
+ return step / velocity
+ elif unit == 'cell':
+ return (step * cell_size) / velocity
Modified: sandbox/annakrat/r3.flow/r3.flow.py
===================================================================
--- sandbox/annakrat/r3.flow/r3.flow.py 2014-06-02 17:03:50 UTC (rev 60672)
+++ sandbox/annakrat/r3.flow/r3.flow.py 2014-06-02 18:40:18 UTC (rev 60673)
@@ -35,16 +35,52 @@
#%option G_OPT_V_OUTPUT
#% description: Name of vector map of flow lines
#%end
+#%option
+#% key: unit
+#% type: string
+#% required: no
+#% answer: cell
+#% options: time,length,cell
+#% descriptions: elapsed time,length in map units,length in cells (voxels)
+#% label: Unit of integration step
+#% description: Default unit is cell
+#% gisprompt: Integration
+#%end
+#%option
+#% key: step
+#% type: double
+#% required: no
+#% label: Integration step in selected unit
+#% description: Default step is 0.5 cell
+#% gisprompt: Integration
+#%end
+#%option
+#% key: limit
+#% type: integer
+#% required: no
+#% answer: 2000
+#% description: Maximum number of steps
+#% gisprompt: Integration
+#%end
+#%flag
+#% key: u
+#% description: Compute upstream flowlines instead of default downstream flowlines
+#%end
+import math
+
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
-from interpolate import integrate_next
+from integrate import integrate_next, get_time_step
+from rast3d_functions import get_velocity, norm
+EPSILON = 1e-7
+
def main():
options, flags = gcore.parser()
@@ -60,9 +96,24 @@
gcore.fatal(_("Error when reading input 3D raster maps"))
map_info = grast3d.raster3d_info(v_x)
- # how to estimate default delta t?
- delta_T = map_info['nsres'] / map_info['max']
+ # integrating forward means upstream, backward means downstream
+ if flags['u']:
+ direction = 1
+ else:
+ direction = -1
+ limit = int(options['limit'])
+ # see which units to use
+ step = options['step']
+ if step:
+ step = float(step)
+ unit = options['unit']
+ else:
+ unit = 'cell'
+ step = 0.25
+ # cell size is the diagonal
+ cell_size = math.sqrt(map_info['nsres'] ** 2 + map_info['ewres'] ** 2 + map_info['tbres'] ** 2)
+
# read seed points
seeds = []
inp_seeds = VectorTopo(options['seed_points'])
@@ -77,17 +128,32 @@
if flowline_vector.exist() and not gcore.overwrite:
gcore.fatal(_("Vector map <{flowlines}> already exists.").format(flowlines=flowline_vector))
flowline_vector.open('w')
+
+ # go through all points and integrate flowline for each
for seed in seeds:
line = Line([Point(seed[0], seed[1], seed[2])])
point = seed
- while True:
- point = integrate_next(map_info, velocity, point, delta_T)
+ count = 0
+ while count <= limit:
+ velocity_vector = get_velocity(map_info, velocity, 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)
+ point = integrate_next(map_info, velocity, point, delta_T * direction)
if point is None:
- if len(line) > 1:
- flowline_vector.write(line)
break
line.append(Point(point[0], point[1], point[2]))
+ count += 1
+ if len(line) > 1:
+ flowline_vector.write(line)
+ gcore.message(_("Flowline ended after %s steps.") % (count - 1), flag='v')
+
flowline_vector.close()
Modified: sandbox/annakrat/r3.flow/rast3d_functions.py
===================================================================
--- sandbox/annakrat/r3.flow/rast3d_functions.py 2014-06-02 17:03:50 UTC (rev 60672)
+++ sandbox/annakrat/r3.flow/rast3d_functions.py 2014-06-02 18:40:18 UTC (rev 60673)
@@ -3,21 +3,21 @@
@author: Anna Petrasova
"""
+import math
#from grass.script import raster3d as grast3d
-def get_velocity(rast3d_region, v_x, v_y, v_z, x, y, z):
+def get_velocity(rast3d_region, velocity_arrays, point):
"""Returns tuple of interpolated velocity values
- in xyz components for a given point"""
- vel_x = rast3d_get_window_value(rast3d_region, v_x, y, x, z)
- vel_y = rast3d_get_window_value(rast3d_region, v_y, y, x, z)
- vel_z = rast3d_get_window_value(rast3d_region, v_z, y, x, z)
+ in xyz components for a given point, or None if outside of region"""
+ return rast3d_trilinear_interpolation(rast3d_region, velocity_arrays,
+ point[1], point[0], point[2])
- if vel_x is None or vel_y is None or vel_z is None:
- return None
- return vel_x, vel_y, vel_z
+def norm(vector):
+ return math.sqrt(vector[0] * vector[0] + vector[1] * vector[1] + vector[2] * vector[2])
+
# following functions are already in raster3d library, this is simpler Python version
# except for trilinear interpolation which could be added there.
@@ -43,13 +43,13 @@
return array[depth][row][col]
-def rast3d_get_window_value(rast3d_region, array, north, east, top):
- """Returns interpolated value in a point defined by geogr. coordinates"""
- return rast3d_trilinear_interpolation(rast3d_region, array, north, east, top)
+#def rast3d_get_window_value(rast3d_region, array, north, east, top):
+# """Returns interpolated value in a point defined by geogr. coordinates"""
+# return rast3d_trilinear_interpolation(rast3d_region, array, north, east, top)
-def rast3d_trilinear_interpolation(rast3d_region, array, north, east, top):
- """Linearly interpolates value in a given point from the 8 closest cells
+def rast3d_trilinear_interpolation(rast3d_region, arrays, north, east, top):
+ """Linearly interpolates value of multiple 3D maps in a given point from the 8 closest cells
based on the centers of cells."""
reg = rast3d_region
# find nearest cells
@@ -63,20 +63,22 @@
(north + reg['nsres'] / 2., east + reg['ewres'] / 2., top + reg['tbres'] / 2.)]
# get values of the nearest cells
- values = []
+ array_values = [[] for i in range(len(arrays))]
for point in points:
col, row, depth = rast3d_location2coord(reg, *point)
- value = rast3d_get_value_region(reg, array, col, row, depth)
- if value is None:
- values.append(0)
- else:
- values.append(value)
+ for i in range(len(arrays)):
+ value = rast3d_get_value_region(reg, arrays[i], col, row, depth)
+ if value is None:
+ array_values[i].append(0)
+ else:
+ array_values[i].append(value)
# compute weights
col, row, depth = rast3d_location2coord(reg, north, east, top)
- if rast3d_get_value_region(reg, array, col, row, depth) is None:
+ # check if we are out of region, any array should work
+ if rast3d_get_value_region(reg, arrays[0], col, row, depth) is None:
return None
-
+
# hm, is there some more elegant way?
# x
temp = east - reg['west'] - col * reg['ewres']
@@ -100,7 +102,10 @@
x * y * z]
# weighted of surrounding values
- value = 0
- for i in range(len(weights)):
- value += weights[i] * values[i]
- return value
+ interpolated_values = []
+ for values in array_values:
+ value = 0
+ for i in range(len(weights)):
+ value += weights[i] * values[i]
+ interpolated_values.append(value)
+ return interpolated_values
Modified: sandbox/annakrat/r3.flow/test_interpolation.py
===================================================================
--- sandbox/annakrat/r3.flow/test_interpolation.py 2014-06-02 17:03:50 UTC (rev 60672)
+++ sandbox/annakrat/r3.flow/test_interpolation.py 2014-06-02 18:40:18 UTC (rev 60673)
@@ -36,12 +36,14 @@
point = [uniform(self.minim + reg['ewres'] / 2., self.maxim - reg['ewres'] / 2.),
uniform(self.minim + reg['nsres'] / 2., self.maxim - reg['nsres'] / 2.),
uniform(self.minim + reg['tbres'] / 2., self.maxim - reg['tbres'] / 2.)]
- tested = rast3d_trilinear_interpolation(reg, rast3d, point[1], point[0], point[2])
+ tested = rast3d_trilinear_interpolation(reg, [rast3d, rast3d, rast3d],
+ point[1], point[0], point[2])
coordinates = [[(point[2] - 0.5 * reg['tbres'] - self.minim) / reg['tbres']],
[(reg['north'] - point[1] - 0.5 * reg['nsres']) / reg['nsres']],
[(point[0] - 0.5 * reg['ewres'] - self.minim) / reg['ewres']]]
reference = ndimage.map_coordinates(rast3d, order=1, coordinates=coordinates)
- self.assertAlmostEqual(tested, reference, places=4)
+ self.assertAlmostEqual(tested[0], reference, places=4)
+ self.assertAlmostEqual(tested[1], reference, places=4)
def tearDown(self):
gcore.run_command('g.remove', rast3d=self.name)
More information about the grass-commit
mailing list