[GRASS-SVN] r60620 - in sandbox/annakrat: . r3.flow
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri May 30 19:11:14 PDT 2014
Author: annakrat
Date: 2014-05-30 19:11:14 -0700 (Fri, 30 May 2014)
New Revision: 60620
Added:
sandbox/annakrat/r3.flow/
sandbox/annakrat/r3.flow/integrate.py
sandbox/annakrat/r3.flow/r3.flow.html
sandbox/annakrat/r3.flow/r3.flow.py
sandbox/annakrat/r3.flow/rast3d_functions.py
sandbox/annakrat/r3.flow/test_interpolation.py
Log:
initial r3.flow GSoC code added to sandbox
Added: sandbox/annakrat/r3.flow/integrate.py
===================================================================
--- sandbox/annakrat/r3.flow/integrate.py (rev 0)
+++ sandbox/annakrat/r3.flow/integrate.py 2014-05-31 02:11:14 UTC (rev 60620)
@@ -0,0 +1,50 @@
+# -*- coding: utf-8 -*-
+"""
+ at author: Anna Petrasova
+"""
+from rast3d_functions import get_velocity
+
+
+# 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 integrate_next(rast3d_region, arrays, point, delta_t):
+ k1 = [None] * 3
+ k2 = [None] * 3
+ k3 = [None] * 3
+ 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)
+ 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)
+ 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)
+ 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])
+ 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
Added: sandbox/annakrat/r3.flow/r3.flow.html
===================================================================
--- sandbox/annakrat/r3.flow/r3.flow.html (rev 0)
+++ sandbox/annakrat/r3.flow/r3.flow.html 2014-05-31 02:11:14 UTC (rev 60620)
@@ -0,0 +1,27 @@
+<h2>DESCRIPTION</h2>
+
+
+
+
+<h2>NOTES</h2>
+
+
+
+<h2>EXAMPLE</h2>
+
+
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="r3.flow.html">r3.flow</a>
+</em>
+
+
+<h2>AUTHOR</h2>
+
+Anna Petrasova<br>
+<i><a href="http://gis.ncsu.edu/osgeorel/">NCSU OSGeoREL</a></i>
+
+<p>
+<i>Last changed: $Date: 2013-02-15 23:08:41 +0100 (Fri, 15 Feb 2013) $</i>
Added: sandbox/annakrat/r3.flow/r3.flow.py
===================================================================
--- sandbox/annakrat/r3.flow/r3.flow.py (rev 0)
+++ sandbox/annakrat/r3.flow/r3.flow.py 2014-05-31 02:11:14 UTC (rev 60620)
@@ -0,0 +1,95 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+############################################################################
+#
+# MODULE: r3.flow
+# AUTHOR: Anna Petrasova
+# PURPOSE:
+#
+#
+#
+# COPYRIGHT: (c) 2014 the GRASS Development Team
+# This program is free software under the GNU General Public
+# License (>=v2). Read the file COPYING that comes with GRASS
+# for details.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+#############################################################################
+#%module
+#% description: Generate 3D flow lines from 3D raster representing velocity field
+#% keywords: raster3d
+#% keywords: vector
+#%end
+#%option G_OPT_R3_INPUT
+#% description: Names of three 3D raster maps describing x, y, z components of vector field
+#%end
+#%option G_OPT_V_INPUT
+#% key: seed_points
+#% description:
+#% label: Name of vector map with points from which flow lines are generated
+#%end
+#%option G_OPT_V_OUTPUT
+#% description: Name of vector map of flow lines
+#%end
+
+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
+
+
+def main():
+ options, flags = gcore.parser()
+
+ try:
+ v_x, v_y, v_z = options['input'].split(',')
+ except ValueError:
+ gcore.fatal(_("Please provide 3 input 3D raster maps."))
+
+ # read input 3D maps
+ velocity = [garray.array3d(), garray.array3d(), garray.array3d()]
+ for inp, vel in zip((v_x, v_y, v_z), velocity):
+ if vel.read(mapname=inp) != 0:
+ 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']
+
+ # read 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()
+
+ # 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')
+ 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)
+ if point is None:
+ if len(line) > 1:
+ flowline_vector.write(line)
+ break
+ line.append(Point(point[0], point[1], point[2]))
+
+ flowline_vector.close()
+
+
+if __name__ == '__main__':
+ main()
Added: sandbox/annakrat/r3.flow/rast3d_functions.py
===================================================================
--- sandbox/annakrat/r3.flow/rast3d_functions.py (rev 0)
+++ sandbox/annakrat/r3.flow/rast3d_functions.py 2014-05-31 02:11:14 UTC (rev 60620)
@@ -0,0 +1,106 @@
+# -*- coding: utf-8 -*-
+"""
+ at author: Anna Petrasova
+"""
+
+#from grass.script import raster3d as grast3d
+
+
+def get_velocity(rast3d_region, v_x, v_y, v_z, x, y, z):
+ """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)
+
+ if vel_x is None or vel_y is None or vel_z is None:
+ return None
+ return vel_x, vel_y, vel_z
+
+
+# following functions are already in raster3d library, this is simpler Python version
+# except for trilinear interpolation which could be added there.
+
+
+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'])
+
+ return col, row, depth
+
+
+def rast3d_get_value_region(rast3d_region, array, col, row, depth):
+ """Returns value of 3D raster array in particular column, row and depth.
+ Returns None if it is outside of array."""
+ if col < 0 or col >= rast3d_region['cols'] or \
+ row < 0 or row >= rast3d_region['rows'] or \
+ depth < 0 or depth >= rast3d_region['depths']:
+ return None
+
+ 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_trilinear_interpolation(rast3d_region, array, north, east, top):
+ """Linearly interpolates value 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
+ values = []
+ 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)
+
+ # compute weights
+ col, row, depth = rast3d_location2coord(reg, north, east, top)
+ if rast3d_get_value_region(reg, array, col, row, depth) is None:
+ return None
+
+ # hm, is there some more elegant way?
+ # 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']
+ rz = 1 - z
+ weights = [rx * ry * rz,
+ x * ry * rz,
+ rx * y * rz,
+ x * y * rz,
+ rx * ry * z,
+ x * ry * z,
+ rx * y * z,
+ x * y * z]
+
+ # weighted of surrounding values
+ value = 0
+ for i in range(len(weights)):
+ value += weights[i] * values[i]
+ return value
Added: sandbox/annakrat/r3.flow/test_interpolation.py
===================================================================
--- sandbox/annakrat/r3.flow/test_interpolation.py (rev 0)
+++ sandbox/annakrat/r3.flow/test_interpolation.py 2014-05-31 02:11:14 UTC (rev 60620)
@@ -0,0 +1,51 @@
+# -*- coding: utf-8 -*-
+"""
+ at brief test trilinear interpolation by comparing to scipy implementation
+ at author: Anna Petrasova
+"""
+import unittest
+from random import uniform
+from grass.script import core as gcore
+from grass.script import array as garray
+from grass.script import raster3d as grast3d
+
+from rast3d_functions import rast3d_trilinear_interpolation
+
+
+class TestTrivariateInterpolation(unittest.TestCase):
+
+ def setUp(self):
+ from scipy import ndimage
+ gcore.use_temp_region()
+ self.name = 'random_voxel'
+ self.minim = -25.
+ self.maxim = 123.
+ self.res = 5.
+
+ def test_interpolation(self):
+ from scipy import ndimage
+ gcore.run_command('g.region', n=self.maxim, s=self.minim, e=self.maxim, w=self.minim,
+ t=self.maxim, b=self.minim, res3=self.res)
+ gcore.run_command('r3.mapcalc', expression=self.name + ' = rand(-20,500)')
+ rast3d = garray.array3d()
+ rast3d.read(mapname=self.name)
+ reg = grast3d.raster3d_info(self.name)
+
+ for i in range(100):
+ # avoid edges because numpy gives always 0
+ 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])
+ 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)
+
+ def tearDown(self):
+ gcore.run_command('g.remove', rast3d=self.name)
+ gcore.del_temp_region()
+
+if __name__ == '__main__':
+ unittest.main()
More information about the grass-commit
mailing list