r.tri added distance weighting to calculation and optional tiled calc…

…ulation (#317)

grass7/raster/r.tri/r.tri.html

@@ -1,28 +1,41 @@
 <h2>DESCRIPTION</h2>
 
-<p><i>r.tri </i>calculates the Terrain Ruggedness Index (TRI) of Riley et al. (1999). The index represents the mean
-    change in elevation between a grid cell and its neighbours, over a user-specified moving window size. The
-    original calculation in Riley et al. (1999) used only a 3x3 neighbourhood and represented the sum of the absolute
-    deviations between the center pixel and its immediate 8 neighbours. In r.tri, this calculation is modified so that
-    the calculation can be extended over any scale by taking the mean of the absolute deviations.</p>
+<p><i>r.tri </i>calculates the Terrain Ruggedness Index (TRI) of Riley et al.
+    (1999). The index represents the mean change in elevation between a grid
+    cell and its neighbours, over a user-specified moving window size. The
+    original calculation in Riley et al. (1999) used only a 3x3 neighbourhood
+    and represented the sum of the absolute deviations between the center pixel
+    and its immediate 8 neighbours. In r.tri, this calculation is modified so
+    that the calculation can be extended over any scale by taking the mean of
+    the absolute deviations.</p>
 
 <h2>NOTES</h2>
 
-<p><i>r.tri</i> produces fairly similar results to the average deviation of elevation values, apart from the center
-    pixel is used in place of the mean. In practice, this produces a slightly less smoothed result that can potentially
-    better highlight finer-scale terrain features.
+<p>
+    <i>r.tri</i> produces fairly similar results to the average deviation of
+    elevation values, apart from the center pixel is used in place of the mean.
+    In practice, this produces a slightly less smoothed result that can
+    potentially highlight finer-scale terrain features. However, because the
+    terrain ruggedness index does not consider the distance of each cell from
+    the centre cell in it's calculation, the TRI results can become noisy with
+    large window sizes. To avoid this, weighting each cell by the inverse of
+    its distance can be used by setting the <i>exponent</i> parameter to &gt 0.
 </p>
 
 <p>
-    Similar to many other GRASS GIS algorithms, cell padding is not performed automatically which will leave null values
-    at the boundaries of the output raster relative to the size of the input raster. To minimize this effect the DEM
-    needs to be buffered/grown prior to using r.tri.
+    Similar to many other GRASS GIS algorithms, cell padding is not performed
+    automatically which will leave null values at the boundaries of the output
+    raster relative to the size of the input raster. To minimize this effect
+    the DEM needs to be buffered/grown prior to using r.tri.
 </p>
 
 <p>
-    To reduce computational times for large raster datasets, setting <em>n_jobs</em> &gt 1 will use a parallelized and
-    tiled calculations that is spread across multiple processing cores. This option requires the
-    GRASS addon <a href="r.mapcalc.tiled.html">r.mapcalc.tiled</a> to be installed.
+    Currently, <i>r.tri</i> is implemented using a <i>r.mapcalc</i> focal
+    function. This becomes slow for large window sizes. To reduce computational
+    times for large raster datasets, setting <em>processes</em> parameter to 
+    &gt 1 will use a parallelized and tiled calculations that is spread across
+    multiple processing cores. This option requires the GRASS addon
+    <a href="r.mapcalc.tiled.html">r.mapcalc.tiled</a> to be installed.
 </p>
 
 <h2>EXAMPLE</h2>

grass7/raster/r.tri/r.tri.py

@@ -32,50 +32,52 @@
 #%option
 #% key: size
 #% type: integer
-#% description: Size of neighbourhood in cells
+#% description: Size of neighbourhood in cells (> 2 and <= 51)
 #% required: yes
 #% answer: 3
 #%end
 
+#%option
+#% key: exponent
+#% type: double
+#% description: Distance weighting exponent (>= 0 and <= 4.0)
+#% required: no
+#% answer: 0.0
+#%end
+
+#%option
+#% key: processes
+#% type: integer
+#% label: Number of processing cores for tiled calculation
+#% description: Number of processing cores for tiled calculation (negative numbers are all cpus -1, -2 etc.)
+#% required: no
+#% answer: 1
+#%end
+
 #%flag
 #% key: c
 #% description: Use circular neighborhood
 #%end
 
+import math
+import multiprocessing as mp
 import sys
-import atexit
-import random
-import string
-import grass.script as gs
-
-TMP_RAST = []
-
-
-def cleanup():
-    gs.message("Deleting intermediate files...")
-
-    for f in TMP_RAST:
-        gs.run_command("g.remove", type="raster", name=f, flags="f", quiet=True)
-
-
-def temp_map(name):
-    tmp = name + "".join(
-        [random.choice(string.ascii_letters + string.digits) for n in range(4)]
-    )
-    TMP_RAST.append(tmp)
 
-    return tmp
+import grass.script as gs
+import numpy as np
+from grass.pygrass.gis.region import Region
+from grass.pygrass.modules.shortcuts import raster as gr
 
 
-def focal_expr(radius, window_square=False):
+def focal_expr(radius, circular=False):
     """Returns array offsets relative to centre pixel (0,0) for a matrix of
     size radius
 
     Args
     ----
     radius : int
         Radius of the focal function
-    window_square : bool. Optional (default is False)
+    circular : bool. Optional (default is False)
         Boolean to use a circular or square window
 
     Returns
@@ -84,70 +86,217 @@ def focal_expr(radius, window_square=False):
         List of pixel positions (row, col) relative to the center pixel
         ( 1, -1)  ( 1, 0)  ( 1, 1)
         ( 0, -1)  ( 0, 0)  ( 0, 1)
-        (-1, -1)  (-1, 0)  (-1, 1)"""
-
-    offsets = []
+        (-1, -1)  (-1, 0)  (-1, 1)
+    """
 
     # generate a list of spatial neighbourhood offsets for the chosen radius
     # ignoring the centre cell
-    if window_square:
+    size = radius * 2 + 1
+    centre = int(size / 2)
 
-        for i in range(-radius, radius + 1):
-            for j in range(-radius, radius + 1):
-                if (i, j) != (0, 0):
-                    offsets.append((i, j))
+    rows, cols = np.ogrid[-radius : radius + 1, -radius : radius + 1]
+    row_offsets, col_offsets = np.meshgrid(rows, cols)
 
+    # create circular mask (also masking centre)
+    if circular:
+        mask = distance_from_centre(radius) <= radius
     else:
+        mask = np.ones((size, size), dtype=np.bool)
+    mask[centre, centre] = False
 
-        for i in range(-radius, radius + 1):
-            for j in range(-radius, radius + 1):
-                row = i + radius
-                col = j + radius
+    # mask and flatten the offsets
+    row_offsets = row_offsets[mask]
+    col_offsets = col_offsets[mask]
 
-                if pow(row - radius, 2) + pow(col - radius, 2) <= pow(radius, 2) and (
-                    i,
-                    j,
-                ) != (0, 0):
-                    offsets.append((j, i))
+    # create a list of offsets
+    offsets = list()
+    for i, j in zip(row_offsets, col_offsets):
+        offsets.append((i, j))
 
     return offsets
 
 
+def distance_from_centre(radius):
+    """Create a square matrix filled with the euclidean distance from the 
+    centre cell
+
+    Parameters
+    ----------
+    radius : int
+        Radius of the square matrix in cells.
+    
+    Returns
+    -------
+    dist_from_centre : 2d ndarray
+        Square matrix with each cell filled with the distance from the centre
+        cell.
+    
+    """
+    size = radius * 2 + 1
+    centre = int(size / 2)
+
+    Y, X = np.ogrid[:size, :size]
+    dist_from_centre = np.sqrt((X - centre) ** 2 + (Y - centre) ** 2)
+
+    return dist_from_centre
+
+
+def idw_weights(radius, p, circular=False):
+    """Create square matrix of inverse distance weights
+
+    The weights are normalized (sum to 1) and are flattened as a list
+
+    Parameters
+    ----------
+    radius : int
+        Radius of the square matrix in cells.
+    
+    p : float
+        Distance weighting power. p=0 gives equal weights.
+    
+    circular : bool (opt). Default is False
+        Optionally use a circular mask.
+    
+    Returns
+    -------
+    W : list
+        Returns a square matrix of weights that excludes the centre cell, or
+        other cells if a circular mask is used. The matrix is flattened and
+        returned as a list.
+    """
+    size = radius * 2 + 1
+    centre = int(size / 2)
+
+    # create distance matrix
+    dist = distance_from_centre(radius)
+
+    # create inverse distance weights
+    W = dist.copy()
+    W[centre, centre] = np.inf
+    W = 1 / (W ** p)
+
+    # normalize weights to sum to 1 (excluding centre)
+    W[centre, centre] = np.inf
+    W = W / W[np.isfinite(W)].sum()
+
+    # circular mask
+    if circular:
+        mask = dist <= radius
+    else:
+        mask = np.isfinite(W)
+    mask[centre, centre] = False
+
+    return list(W[mask])
+
+
+def tile_shape(region, n_jobs):
+    """Calculates the number of tiles required for one tile per cpu
+
+    Parameters
+    ----------
+    region : pygrass.gis.region.Region
+        The computational region object.
+    
+    n_jobs : int
+        The number of processing cores.
+    
+    Returns
+    -------
+    width, height : tuple
+        The width and height of each tile.
+    """
+    n = math.sqrt(n_jobs)
+    width = int(math.ceil(region.cols / n))
+    height = int(math.ceil(region.rows / n))
+
+    return width, height
+
+
 def main():
     dem = options["input"]
     tri = options["output"]
     size = int(options["size"])
+    exponent = float(options["exponent"])
+    processes = int(options["processes"])
     circular = flags["c"]
     radius = int((size - 1) / 2)
-
+
+    region = Region()
+
+    # Some checks
     if "@" in tri:
-    	tri = tri.split("@")[0]
+        tri = tri.split("@")[0]
+
+    if processes == 0:
+        gs.fatal("Number of processing cores for parallel computation must not equal 0")
 
-    # calculate TRI based on map calc statements
+    if processes < 0:
+        system_cores = mp.cpu_count()
+        processes = system_cores + processes + 1
+
+    if processes > 1:
+        if gs.find_program("r.mapcalc.tiled") is False:
+            gs.fatal(
+                "The GRASS addon r.mapcalc.tiled must also be installed if n_jobs != 1. Run 'g.extension r.mapcalc.tiled'"
+            )
+
+    if size <= 2 or size > 51:
+        gs.fatal("size must be > 2 and <= 51")
+
+    if size % 2 != 1:
+        gs.fatal("size must be an odd number")
+
+    if exponent < 0 or exponent > 4.0:
+        gs.fatal("exponent must be >= 0 and <= 4.0")
+
+    # Calculate TRI based on map calc statements
     gs.message("Calculating the Topographic Ruggedness Index...")
 
-    # generate a list of spatial neighbourhood offsets for the chosen radius
+    # Generate a list of spatial neighbourhood offsets for the chosen radius
     # ignoring the center cell
-    offsets = focal_expr(radius=radius, window_square=not circular)
+    offsets = focal_expr(radius, circular)
+    weights = idw_weights(radius, exponent, circular)
     terms = []
-    for d in offsets:
+    for d, w in zip(offsets, weights):
         d_str = ",".join(map(str, d))
-        terms.append("({dem}[{d}]-{dem})^2".format(dem=dem, d=d_str))
+        terms.append("{w}*abs({dem}[{d}]-{dem})".format(dem=dem, d=d_str, w=w))
 
-    # define the calculation expression
-    ncells = len(offsets) + 1
-    terms = " + ".join(terms)
+    # Define the calculation expression
+    terms = "+".join(terms)
 
-    # perform the r.mapcalc calculation with the moving window
-    expr = "{tri} = float(sqrt(({terms}) / {ncells}))".format(
-        tri=tri, ncells=ncells, terms=terms
+    # Perform the r.mapcalc calculation with the moving window
+    expr = "{tri} = float({terms})".format(tri=tri, terms=terms)
+
+    width, height = tile_shape(region, processes)
+
+    if width < region.cols and height < region.rows and processes > 1:
+        gr.mapcalc_tiled(
+            expression=expr,
+            width=width,
+            height=height,
+            processes=processes,
+            overlap=int((size + 1) / 2),
+            output=tri,
+        )
+    else:
+        gs.mapcalc(expr)
+
+    # update metadata
+    opts = ""
+    if circular:
+        opts = "-c"
+
+    gr.support(
+        map=tri,
+        title="Terrain Ruggedness Index",
+        description="Generated with r.tri input={dem} output={tri} size={size} exponent={exponent} {flags}".format(
+            dem=dem, tri=tri, size=size, exponent=exponent, flags=opts
+        ),
     )
-    gs.mapcalc(expr)
 
     return 0
 
 
 if __name__ == "__main__":
     options, flags = gs.parser()
-    atexit.register(cleanup)
     sys.exit(main())