[GRASS-SVN] r61020 - sandbox/annakrat/r3.flow
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri Jun 27 14:18:37 PDT 2014
Author: annakrat
Date: 2014-06-27 14:18:37 -0700 (Fri, 27 Jun 2014)
New Revision: 61020
Modified:
sandbox/annakrat/r3.flow/gradient.py
sandbox/annakrat/r3.flow/integrate.py
sandbox/annakrat/r3.flow/r3.flow.py
sandbox/annakrat/r3.flow/rast3d_functions.py
Log:
r3.flow: compute gradient on the fly, currently doesn't yet produce identical results when gradient is precomputed
Modified: sandbox/annakrat/r3.flow/gradient.py
===================================================================
--- sandbox/annakrat/r3.flow/gradient.py 2014-06-27 20:50:02 UTC (rev 61019)
+++ sandbox/annakrat/r3.flow/gradient.py 2014-06-27 21:18:37 UTC (rev 61020)
@@ -57,7 +57,7 @@
gradient_array_z[depth, row, col] = (depth_array[depth + 1] -
depth_array[depth - 1]) / (2 * step[2])
- return gradient_array_x, gradient_array_y, gradient_array_z
+ return gradient_array_x, -gradient_array_y, gradient_array_z
class TestGradient(unittest.TestCase):
@@ -72,7 +72,7 @@
gradient_x, gradient_y, gradient_z = gradient(array, step)
self.assertEqual(True, np.allclose(gradient_z, numpy_gradient_z, atol=1e-8))
- self.assertEqual(True, np.allclose(gradient_y, numpy_gradient_y, atol=1e-8))
+ self.assertEqual(True, np.allclose(gradient_y, -numpy_gradient_y, atol=1e-8))
self.assertEqual(True, np.allclose(gradient_x, numpy_gradient_x, atol=1e-8))
Modified: sandbox/annakrat/r3.flow/integrate.py
===================================================================
--- sandbox/annakrat/r3.flow/integrate.py 2014-06-27 20:50:02 UTC (rev 61019)
+++ sandbox/annakrat/r3.flow/integrate.py 2014-06-27 21:18:37 UTC (rev 61020)
@@ -18,11 +18,11 @@
# nice explanation of runge-kutta with picture:
# http://www.marekfiser.com/Projects/Vector-field-visualization-on-GPU-using-CUDA/
#4-Vector-field-integrators-for-stream-line-visualization
-def rk4_integrate_next(rast3d_region, arrays, point, next_point, delta_t):
+def rk4_integrate_next(rast3d_region, velocity_obj, point, next_point, delta_t):
"""Integrates new point of a flowline using 4th order Runge-Kutta method.
:param dict rast3d_region: information about region
- :param arrays: 3D arrays of velocity components
+ :param velocity_obj: velocity class to handle different sources of velocity field
:param list point: list of x, y, z coordinates from where to integrate
:param list next_point: list of x, y, z coordinates of the integrated point
(computed here and the list is modified)
@@ -35,27 +35,27 @@
k4 = [None] * 3
x, y, z = point[0], point[1], point[2]
# first
- velocity = get_velocity(rast3d_region, arrays, point)
+ velocity = get_velocity(rast3d_region, velocity_obj, point)
if velocity is None:
return False
for i in range(3):
k1[i] = delta_t * velocity[i]
# second
- velocity = get_velocity(rast3d_region, arrays,
+ velocity = get_velocity(rast3d_region, velocity_obj,
(x + k1[0] / 2, y + k1[1] / 2, z + k1[2] / 2))
if velocity is None:
return False
for i in range(3):
k2[i] = delta_t * velocity[i]
# third
- velocity = get_velocity(rast3d_region, arrays,
+ velocity = get_velocity(rast3d_region, velocity_obj,
(x + k2[0] / 2, y + k2[1] / 2, z + k2[2] / 2))
if velocity is None:
return False
for i in range(3):
k3[i] = delta_t * velocity[i]
# fourth
- velocity = get_velocity(rast3d_region, arrays, (x + k3[0], y + k3[1], z + k3[2]))
+ velocity = get_velocity(rast3d_region, velocity_obj, (x + k3[0], y + k3[1], z + k3[2]))
if velocity is None:
return False
for i in range(3):
@@ -87,12 +87,12 @@
-277. / 14336, 512. / 1771 - 1. / 4]
-def rk45_integrate_next(rast3d_region, arrays, point, next_point, delta_t, min_step, max_step):
+def rk45_integrate_next(rast3d_region, velocity_obj, point, next_point, delta_t, min_step, max_step):
"""Integrates new point of a flowline using Runge-Kutta method with adaptive step size.
Uses Cash-Karp coefficients. Tries to decrease the estimated error below max error.
:param dict rast3d_region: information about region
- :param arrays: 3D arrays of velocity components
+ :param velocity_obj: velocity class to handle different sources of velocity field
:param list point: list of x, y, z coordinates from where to integrate
:param list next_point: list of x, y, z coordinates of the integrated point
(computed here and the list is modified)
@@ -109,7 +109,7 @@
# try to iteratively decrease error to less than max error
while estimated_error > MAX_ERROR:
# compute next point and get estimated error
- if rk45_next(rast3d_region, arrays, point, next_point, delta_t, error):
+ if rk45_next(rast3d_region, velocity_obj, point, next_point, delta_t, error):
estimated_error = norm(error)
else:
return False
@@ -135,19 +135,19 @@
delta_t = tmp
# break when the adjustment was needed (not sure why)
if do_break:
- if not rk45_next(rast3d_region, arrays, point, next_point, delta_t, error):
+ if not rk45_next(rast3d_region, velocity_obj, point, next_point, delta_t, error):
return False
break
return True
-def rk45_next(rast3d_region, arrays, point, next_point, delta_t, error):
+def rk45_next(rast3d_region, velocity_obj, point, next_point, delta_t, error):
"""The actual method which integrates new point of a flowline
using Runge-Kutta method with Cash-Karp coefficients. Provides error estimate.
:param dict rast3d_region: information about region
- :param arrays: 3D arrays of velocity components
+ :param velocity_obj: velocity class to handle different sources of velocity field
:param list point: list of x, y, z coordinates from where to integrate
:param list next_point: list of x, y, z coordinates of the integrated point
(computed here and the list is modified)
@@ -158,7 +158,7 @@
# 3 is 3 dimensions, 6 is the number of k's
tmp1 = [[0] * 3] * 6
tmp_point = [0] * 3
- velocity = get_velocity(rast3d_region, arrays, point)
+ velocity = get_velocity(rast3d_region, velocity_obj, point)
if velocity is None:
return False
tmp1[0] = velocity
@@ -171,7 +171,7 @@
tmp_point[j] = point[j] + delta_t * sum_tmp
- velocity = get_velocity(rast3d_region, arrays, tmp_point)
+ velocity = get_velocity(rast3d_region, velocity_obj, tmp_point)
if velocity is None:
return False
tmp1[i] = velocity
Modified: sandbox/annakrat/r3.flow/r3.flow.py
===================================================================
--- sandbox/annakrat/r3.flow/r3.flow.py 2014-06-27 20:50:02 UTC (rev 61019)
+++ sandbox/annakrat/r3.flow/r3.flow.py 2014-06-27 21:18:37 UTC (rev 61020)
@@ -99,7 +99,7 @@
from grass.pygrass.vector.geometry import Line, Point
from integrate import rk4_integrate_next, rk45_integrate_next, get_time_step, MIN_STEP, MAX_STEP
-from rast3d_functions import get_velocity, norm, compute_gradient, rast3d_location2coord
+from rast3d_functions import get_velocity, norm, rast3d_location2coord
from voxel_traversal import traverse
EPSILON = 1e-8
@@ -115,27 +115,40 @@
self.flowaccum = flowaccum
+class Velocity:
+ def __init__(self):
+ self.compute_gradient = False
+ self.vector_arrays = None
+ self.scalar_array = None
+ self.last_neighbors_values = None
+ self.last_neighbors_extent = None
+ self.last_shift = None
+
+
def main():
options, flags = gcore.parser()
+ velocity_obj = Velocity()
if options['vector_field']:
try:
v_x, v_y, v_z = options['vector_field'].split(',')
except ValueError:
gcore.fatal(_("Please provide 3 input 3D raster maps representing components of vector field."))
# read input 3D maps
- velocity = [garray.array3d(), garray.array3d(), garray.array3d()]
- for inp, vel in zip((v_x, v_y, v_z), velocity):
+ velocity_arrays = [garray.array3d(), garray.array3d(), garray.array3d()]
+ for inp, vel in zip((v_x, v_y, v_z), velocity_arrays):
if vel.read(mapname=inp) != 0:
gcore.fatal(_("Error when reading input 3D raster maps"))
map_info = grast3d.raster3d_info(v_x)
+ velocity_obj.vector_arrays = velocity_arrays
+ velocity_obj.compute_gradient = False
elif options['input']:
- gcore.info(_("Computing gradient field..."))
- array = garray.array3d()
- array.read(mapname=options['input'])
- velocity = compute_gradient(array)
+ scalar_array = garray.array3d()
+ scalar_array.read(mapname=options['input'])
map_info = grast3d.raster3d_info(options['input'])
+ velocity_obj.compute_gradient = True
+ velocity_obj.scalar_array = scalar_array
else:
# this should be handled by parser (in the future)
gcore.fatal("Either input or vector_field options must be specified.")
@@ -181,7 +194,7 @@
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
+ 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)
@@ -208,7 +221,6 @@
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()
@@ -228,7 +240,7 @@
last_coords = (None, None, None)
new_point = [None, None, None]
while count <= limit:
- velocity_vector = get_velocity(map_info, velocity, point)
+ velocity_vector = get_velocity(map_info, velocity_obj, point)
if velocity_vector is None:
break # outside region
velocity_norm = norm(velocity_vector)
@@ -239,12 +251,12 @@
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, point, new_point,
+ 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, point, new_point,
+ ret = rk45_integrate_next(map_info, velocity_obj, point, new_point,
delta_T * direction, min_step, max_step)
if not ret:
break
Modified: sandbox/annakrat/r3.flow/rast3d_functions.py
===================================================================
--- sandbox/annakrat/r3.flow/rast3d_functions.py 2014-06-27 20:50:02 UTC (rev 61019)
+++ sandbox/annakrat/r3.flow/rast3d_functions.py 2014-06-27 21:18:37 UTC (rev 61020)
@@ -5,16 +5,97 @@
import math
import numpy
+from gradient import gradient
#from grass.script import raster3d as grast3d
-def get_velocity(rast3d_region, velocity_arrays, point):
- """Returns tuple of interpolated velocity values
- 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])
+def get_velocity(rast3d_region, velocity_obj, point):
+ if velocity_obj.compute_gradient:
+ return get_gradient(rast3d_region, velocity_obj, point)
+ else:
+ return rast3d_trilinear_interpolation(rast3d_region, velocity_obj.vector_arrays,
+ point[1], point[0], point[2])
+#def get_velocity(rast3d_region, velocity_arrays, point):
+# """Returns tuple of interpolated velocity values
+# 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])
+
+
+def get_gradient(rast3d_region, velocity_obj, point):
+ nearest = find_nearest_voxels(rast3d_region, point[1], point[0], point[2])
+
+ minx, maxx = nearest[0][0], nearest[7][0]
+ miny, maxy = nearest[7][1], nearest[0][1]
+ minz, maxz = nearest[0][2], nearest[7][2]
+ # min, max of x, y, z neighboring voxels
+ neighbors_extent = [minx, maxx,
+ miny, maxy,
+ minz, maxz]
+ if not velocity_obj.last_neighbors_extent or neighbors_extent != velocity_obj.last_neighbors_extent:
+ velocity_obj.last_neighbors_extent = neighbors_extent
+
+ if minx < 0 or maxx >= rast3d_region['cols'] or \
+ miny < 0 or maxy >= rast3d_region['rows'] or \
+ minz < 0 or maxz >= rast3d_region['depths']:
+ return None
+
+ if minx == 0:
+ maxx = minx + 3
+ xshift = 0
+ elif maxx == rast3d_region['cols'] - 1:
+ minx = maxx - 3
+ xshift = 1
+ else:
+ minx -= 1
+ maxx += 1
+ xshift = 2
+
+ if miny == 0:
+ maxy = miny + 3
+ yshift = 0
+ elif maxy == rast3d_region['rows'] - 1:
+ miny = maxy - 3
+ yshift = 1
+ else:
+ miny -= 1
+ maxy += 1
+ yshift = 2
+
+ if minz == 0:
+ maxz = minz + 3
+ zshift = 0
+ elif maxz == rast3d_region['depths'] - 1:
+ minz = maxz - 3
+ zshift = 1
+ else:
+ minz -= 1
+ maxz += 1
+ zshift = 2
+
+ array = velocity_obj.scalar_array[minz:maxz + 1, miny:maxy + 1, minx:maxx + 1]
+# print (minx, maxx), (miny, maxy), (minz, maxz)
+ xyz_gradients = gradient(array, step=[rast3d_region['ewres'], rast3d_region['nsres'],
+ rast3d_region['tbres']])
+
+ array_values = [[] for i in range(3)]
+ for i in range(3):
+ for z in range(2):
+ for y in range(2):
+ for x in range(2):
+ array_values[i].append(xyz_gradients[i][z + zshift][y + yshift][x + xshift])
+
+ velocity_obj.last_neighbors_values = array_values
+ else:
+ array_values = velocity_obj.last_neighbors_values
+
+ x, y, z = get_relative_coords_for_interp(rast3d_region, point[1], point[0], point[2])
+
+ return trilinear_interpolation(array_values, x, y, z)
+
+
def norm(vector):
return math.sqrt(vector[0] * vector[0] + vector[1] * vector[1] + vector[2] * vector[2])
@@ -59,6 +140,27 @@
return array[depth][row][col]
+def find_nearest_voxels(rast3d_region, north, east, top):
+ """Finds 8 nearest voxel to given point and returns their indices."""
+ reg = rast3d_region
+ n_minus, n_plus = north - reg['nsres'] / 2., north + reg['nsres'] / 2.
+ e_minus, e_plus = east - reg['ewres'] / 2., east + reg['ewres'] / 2.
+ t_minus, t_plus = top - reg['tbres'] / 2., top + reg['tbres'] / 2.
+ # find nearest cells
+ points = [(n_minus, e_minus, t_minus),
+ (n_minus, e_plus, t_minus),
+ (n_plus, e_minus, t_minus),
+ (n_plus, e_plus, t_minus),
+ (n_minus, e_minus, t_plus),
+ (n_minus, e_plus, t_plus),
+ (n_plus, e_minus, t_plus),
+ (n_plus, e_plus, t_plus)]
+ coordinates = []
+ for point in points:
+ coordinates.append(rast3d_location2coord(reg, *point))
+
+ return coordinates
+
#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)
@@ -68,20 +170,10 @@
"""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
- points = [(north - reg['nsres'] / 2., east - reg['ewres'] / 2., top - reg['tbres'] / 2.),
- (north - reg['nsres'] / 2., east + reg['ewres'] / 2., top - reg['tbres'] / 2.),
- (north + reg['nsres'] / 2., east - reg['ewres'] / 2., top - reg['tbres'] / 2.),
- (north + reg['nsres'] / 2., east + reg['ewres'] / 2., top - reg['tbres'] / 2.),
- (north - reg['nsres'] / 2., east - reg['ewres'] / 2., top + reg['tbres'] / 2.),
- (north - reg['nsres'] / 2., east + reg['ewres'] / 2., top + reg['tbres'] / 2.),
- (north + reg['nsres'] / 2., east - reg['ewres'] / 2., top + reg['tbres'] / 2.),
- (north + reg['nsres'] / 2., east + reg['ewres'] / 2., top + reg['tbres'] / 2.)]
# get values of the nearest cells
array_values = [[] for i in range(len(arrays))]
- for point in points:
- col, row, depth = rast3d_location2coord(reg, *point)
+ for col, row, depth in find_nearest_voxels(reg, north, east, top):
for i in range(len(arrays)):
value = rast3d_get_value_region(reg, arrays[i], col, row, depth)
if value is None:
@@ -95,18 +187,29 @@
if rast3d_get_value_region(reg, arrays[0], col, row, depth) is None:
return None
- # hm, is there some more elegant way?
- # x
+ x, y, z = get_relative_coords_for_interp(rast3d_region, north, east, top)
+ return trilinear_interpolation(array_values, x, y, z)
+
+
+def get_relative_coords_for_interp(rast3d_region, north, east, top):
+ reg = rast3d_region
+ col, row, depth = rast3d_location2coord(rast3d_region, north, east, top)
+ # x
temp = east - reg['west'] - col * reg['ewres']
x = (temp - reg['ewres'] / 2. if temp > reg['ewres'] / 2. else temp + reg['ewres'] / 2.) / reg['ewres']
- rx = 1 - x
# y
temp = north - reg['south'] - (reg['rows'] - row - 1) * reg['nsres']
y = (temp - reg['nsres'] / 2. if temp > reg['nsres'] / 2. else temp + reg['nsres'] / 2.) / reg['nsres']
- ry = 1 - y
# z
temp = top - reg['bottom'] - depth * reg['tbres']
z = (temp - reg['tbres'] / 2. if temp > reg['tbres'] / 2. else temp + reg['tbres'] / 2.) / reg['tbres']
+
+ return x, y, z
+
+
+def trilinear_interpolation(array_values, x, y, z):
+ rx = 1 - x
+ ry = 1 - y
rz = 1 - z
weights = [rx * ry * rz,
x * ry * rz,
More information about the grass-commit
mailing list