[GRASS-SVN] r60824 - sandbox/annakrat/r3.flow

[svn_grass at osgeo.org]

Author: annakrat
Date: 2014-06-13 14:29:09 -0700 (Fri, 13 Jun 2014)
New Revision: 60824

Added:
   sandbox/annakrat/r3.flow/voxel_traversal.py
Modified:
   sandbox/annakrat/r3.flow/r3.flow.py
Log:
r3.flow: implemented voxel traversal algorithm for precise flow accumulation computation

Modified: sandbox/annakrat/r3.flow/r3.flow.py
===================================================================
--- sandbox/annakrat/r3.flow/r3.flow.py	2014-06-13 20:52:10 UTC (rev 60823)
+++ sandbox/annakrat/r3.flow/r3.flow.py	2014-06-13 21:29:09 UTC (rev 60824)
@@ -66,7 +66,7 @@
 #% type: double
 #% required: no
 #% label: Integration step in selected unit
-#% description: Default step is 0.5 cell
+#% description: Default step is 0.25 cell
 #% gisprompt: Integration
 #%end
 #%option
@@ -254,7 +254,5 @@
     if options['flowaccumulation']:
         flowacc.write(mapname=options['flowaccumulation'])
 
-
-
 if __name__ == '__main__':
     main()

Added: sandbox/annakrat/r3.flow/voxel_traversal.py
===================================================================
--- sandbox/annakrat/r3.flow/voxel_traversal.py	                        (rev 0)
+++ sandbox/annakrat/r3.flow/voxel_traversal.py	2014-06-13 21:29:09 UTC (rev 60824)
@@ -0,0 +1,106 @@
+# -*- coding: utf-8 -*-
+"""
+ at brief voxel traversal algorithm for precise flowaccumulation computing
+ at author: Anna Petrasova
+"""
+
+import numpy as np
+import math
+
+
+def traverse(start, end, region_info):
+    """!Computes which voxel cells are traversed by a line.
+    To each transected cell we add 1.
+
+    :param start: start point (3 coordinates)
+    :param end: end point (3 coordinates)
+    :param dict region_info: information about region
+    :return: coordinates of transected voxel cells
+    """
+    # initialize
+    dx = end[0] - start[0]
+    dy = end[1] - start[1]
+    dz = end[2] - start[2]
+
+    X = int(math.floor((start[0] - region_info['west']) / region_info['ewres']))
+    Y = int(math.floor((start[1] - region_info['south']) / region_info['nsres']))
+    Z = int(math.floor((start[2] - region_info['bottom']) / region_info['tbres']))
+    X_end = int(math.floor((end[0] - region_info['west']) / region_info['ewres']))
+    Y_end = int(math.floor((end[1] - region_info['south']) / region_info['nsres']))
+    Z_end = int(math.floor((end[2] - region_info['bottom']) / region_info['tbres']))
+
+    tDeltaX = region_info['ewres'] / float(dx) if dx != 0 else 1e6
+    tDeltaY = region_info['nsres'] / float(dy) if dy != 0 else 1e6
+    tDeltaZ = region_info['tbres'] / float(dz) if dz != 0 else 1e6
+
+    xtmp = start[0] / region_info['ewres']
+    ytmp = start[1] / region_info['nsres']
+    ztmp = start[2] / region_info['tbres']
+    xtmp = xtmp if xtmp > 0 else -xtmp
+    ytmp = ytmp if ytmp > 0 else -ytmp
+    ztmp = ztmp if ztmp > 0 else -ztmp
+    tMaxX = tDeltaX * (1.0 - (xtmp - int(xtmp)))
+    tMaxY = tDeltaY * (1.0 - (ytmp - int(ytmp)))
+    tMaxZ = tDeltaZ * (1.0 - (ztmp - int(ztmp)))
+
+    stepX = 1 if start[0] < end[0] else -1
+    stepY = 1 if start[1] < end[1] else -1
+    stepZ = 1 if start[2] < end[2] else -1
+
+    coordinates = []
+    while True:
+        if abs(tMaxX) < abs(tMaxY):
+            if abs(tMaxX) < abs(tMaxZ):
+                tMaxX = tMaxX + tDeltaX
+                X = X + stepX
+                if (X - X_end) * stepX > 0:
+                    break
+            else:
+                tMaxZ = tMaxZ + tDeltaZ
+                Z = Z + stepZ
+                if (Z - Z_end) * stepZ > 0:
+                    break
+        else:
+            if abs(tMaxY) < abs(tMaxZ):
+                tMaxY = tMaxY + tDeltaY
+                Y = Y + stepY
+                if (Y - Y_end) * stepY > 0:
+                    break
+            else:
+                tMaxZ = tMaxZ + tDeltaZ
+                Z = Z + stepZ
+                if (Z - Z_end) * stepZ > 0:
+                    break
+
+        coordinates.append((X, Y, Z))
+    return coordinates
+
+if __name__ == '__main__':
+    info = {}
+    info['cols'] = 10
+    info['rows'] = 10
+    info['depths'] = 10
+    info['nsres'] = 4
+    info['ewres'] = 8
+    info['tbres'] = 1
+    info['west'] = -80
+    info['south'] = -40
+    info['bottom'] = -10
+    start = [-66.2,-38.9,-9.5]
+    end = [-25,-24.5,-9.5]
+    X = int(math.floor((start[0] - info['west']) / info['ewres']))
+    Y = int(math.floor((start[1] - info['south']) / info['nsres']))
+    Z = int(math.floor((start[2] - info['bottom']) / info['tbres']))
+    X_end = int(math.floor((end[0] - info['west']) / info['ewres']))
+    Y_end = int(math.floor((end[1] - info['south']) / info['nsres']))
+    Z_end = int(math.floor((end[2] - info['bottom']) / info['tbres']))
+#    start = [1, 1, 0.5]
+#    end = [74, 42, 0.5]
+
+    array = np.zeros((info['depths'], info['rows'], info['cols']), dtype=np.int8)
+    coords = traverse(start, end, info)
+    array[Z, array.shape[1] - Y - 1, X] = 1
+    array[Z_end, array.shape[1] - Y_end - 1, X_end] = 3
+    for coord in coords:
+        array[coord[2], array.shape[1] - coord[1] - 1, coord[0]] = 2
+    q = array[0]