I. natcap#1104 Test multiple decay filtering

Added euclidean distance transform methods for filtering (decaying)
threat rasters to test against the convolution implementation. Also
tested an exponential decay kernel with an add 2.99 constant as
presented in the User's Guide.

src/natcap/invest/habitat_quality.py

@@ -3,6 +3,7 @@
 import collections
 import logging
 import os
+import math
 
 import numpy
 import pygeoprocessing
@@ -539,32 +540,71 @@ def execute(args):
                 kernel_dir, f'kernel_{threat}{lulc_key}{file_suffix}.tif')
 
             decay_type = threat_data['decay']
-
-            create_kernel_task = task_graph.add_task(
-                func=_create_decay_kernel,
-                args=((threat_raster_path, 1), kernel_path, decay_type,
-                      threat_data['max_dist']),
-                target_path_list=[kernel_path],
-                dependent_task_list=[align_task],
-                task_name=f'decay_kernel_{decay_type}{lulc_key}_{threat}')
-
-            filtered_threat_raster_path = os.path.join(
-                intermediate_output_dir,
-                f'filtered_{threat}{lulc_key}{file_suffix}.tif')
-
-            convolve_task = task_graph.add_task(
-                func=pygeoprocessing.convolve_2d,
-                args=((threat_raster_path, 1), (kernel_path, 1),
-                      filtered_threat_raster_path),
-                kwargs={
-                    'target_nodata': _OUT_NODATA,
-                    'ignore_nodata_and_edges': False,
-                    'mask_nodata': False
-                },
-                target_path_list=[filtered_threat_raster_path],
-                dependent_task_list=[create_kernel_task],
-                task_name=f'convolve_{decay_type}{lulc_key}_{threat}')
-            threat_convolve_task_list.append(convolve_task)
+            if 'euclidean' in decay_type:
+                distance_raster_path = os.path.join(
+                    kernel_dir, f'edt_{threat}{lulc_key}{file_suffix}.tif')
+
+                dist_edt_task = task_graph.add_task(
+                    func=pygeoprocessing.distance_transform_edt,
+                    args=((threat_raster_path, 1), distance_raster_path),
+                    target_path_list=[distance_raster_path],
+                    dependent_task_list=[align_task],
+                    task_name=f'distance edt {lulc_key} {threat}')
+
+                if decay_type == 'euclidean-exponential':
+                    ### Decay using distance transform w/ exponential scalar ###
+                    filtered_threat_raster_path = os.path.join(
+                        intermediate_output_dir,
+                        f'filtered_{decay_type}_{threat}{lulc_key}{file_suffix}.tif')
+
+                    dist_exp_task = task_graph.add_task(
+                        func=dist_edt_exp,
+                        args=(
+                            distance_raster_path, 2.99, threat_data['max_dist'],
+                            filtered_threat_raster_path),
+                        target_path_list=[filtered_threat_raster_path],
+                        dependent_task_list=[dist_edt_task],
+                        task_name=f'distance exp {lulc_key} {threat}')
+                    threat_convolve_task_list.append(dist_exp_task)
+                elif decay_type == 'euclidean-linear':
+                    ### Decay using distance transform w/ exponential scalar ###
+                    filtered_threat_raster_path = os.path.join(
+                        intermediate_output_dir,
+                        f'filtered_{decay_type}_{threat}{lulc_key}{file_suffix}.tif')
+
+                    dist_exp_task = task_graph.add_task(
+                        func=dist_edt_linear,
+                        args=(distance_raster_path, threat_data['max_dist'], filtered_threat_raster_path),
+                        target_path_list=[filtered_threat_raster_path],
+                        dependent_task_list=[dist_edt_task],
+                        task_name='distance exp')
+                    threat_convolve_task_list.append(dist_exp_task)
+            else:
+                create_kernel_task = task_graph.add_task(
+                    func=_create_decay_kernel,
+                    args=((threat_raster_path, 1), kernel_path, decay_type,
+                          threat_data['max_dist']),
+                    target_path_list=[kernel_path],
+                    dependent_task_list=[align_task],
+                    task_name=f'decay_kernel_{decay_type}{lulc_key}_{threat}')
+
+                filtered_threat_raster_path = os.path.join(
+                    intermediate_output_dir,
+                    f'filtered_{threat}{lulc_key}{file_suffix}.tif')
+
+                convolve_task = task_graph.add_task(
+                    func=pygeoprocessing.convolve_2d,
+                    args=((threat_raster_path, 1), (kernel_path, 1),
+                          filtered_threat_raster_path),
+                    kwargs={
+                        'target_nodata': _OUT_NODATA,
+                        'ignore_nodata_and_edges': False,
+                        'mask_nodata': False
+                    },
+                    target_path_list=[filtered_threat_raster_path],
+                    dependent_task_list=[create_kernel_task],
+                    task_name=f'convolve_{decay_type}{lulc_key}_{threat}')
+                threat_convolve_task_list.append(convolve_task)
 
             # create sensitivity raster based on threat
             sens_raster_path = os.path.join(
@@ -894,6 +934,8 @@ def _create_decay_kernel(raster_path_band, kernel_path, decay_type, max_dist):
         decay_func = _make_linear_decay_kernel_path
     elif decay_type == 'exponential':
         decay_func = utils.exponential_decay_kernel_raster
+    elif decay_type == 'exponential-scalar':
+        decay_func = _exponential_decay_scalar_kernel_raster
     else:
         raise ValueError(
             "Unknown type of decay in biophysical table, should be"
@@ -1163,7 +1205,8 @@ def validate(args, limit_to=None):
 
     return validation_warnings
 
-def _exponential_decay_kernel_raster(expected_distance, kernel_filepath):
+
+def _exponential_decay_scalar_kernel_raster(expected_distance, kernel_filepath):
     """Create a raster-based exponential decay kernel.
 
     The raster created will be a tiled GeoTiff, with 256x256 memory blocks.
@@ -1234,7 +1277,9 @@ def _exponential_decay_kernel_raster(expected_distance, kernel_filepath):
                 row_indices, col_indices)
             kernel = numpy.where(
                 kernel_index_distances > max_distance, 0.0,
-                numpy.exp(-kernel_index_distances / expected_distance))
+                numpy.exp(
+                    (-kernel_index_distances * 2.99) / expected_distance))
+            integration += numpy.sum(kernel)
 
             kernel_band.WriteArray(kernel, xoff=col_offset,
                                    yoff=row_offset)
@@ -1256,96 +1301,54 @@ def _exponential_decay_kernel_raster(expected_distance, kernel_filepath):
     kernel_band = None
     kernel_dataset = None
 
-def _exponential_decay_scalar_kernel_raster(expected_distance, kernel_filepath):
-    """Create a raster-based exponential decay kernel.
 
-    The raster created will be a tiled GeoTiff, with 256x256 memory blocks.
-
-    Args:
-        expected_distance (int or float): The distance (in pixels) of the
-            kernel's radius, the distance at which the value of the decay
-            function is equal to `1/e`.
-        kernel_filepath (string): The path to the file on disk where this
-            kernel should be stored.  If this file exists, it will be
-            overwritten.
-
-    Returns:
-        None
-    """
-    max_distance = expected_distance * 5
-    kernel_size = int(numpy.round(max_distance * 2 + 1))
-
-    driver = gdal.GetDriverByName('GTiff')
-    kernel_dataset = driver.Create(
-        kernel_filepath.encode('utf-8'), kernel_size, kernel_size, 1,
-        gdal.GDT_Float32, options=[
-            'BIGTIFF=IF_SAFER', 'TILED=YES', 'BLOCKXSIZE=256',
-            'BLOCKYSIZE=256'])
-
-    # Make some kind of geotransform, it doesn't matter what but
-    # will make GIS libraries behave better if it's all defined
-    kernel_dataset.SetGeoTransform([0, 1, 0, 0, 0, -1])
-    srs = osr.SpatialReference()
-    srs.SetWellKnownGeogCS('WGS84')
-    kernel_dataset.SetProjection(srs.ExportToWkt())
-
-    kernel_band = kernel_dataset.GetRasterBand(1)
-    kernel_band.SetNoDataValue(-9999)
+def dist_edt_exp(dist_raster_path, scalar, max_dist, target_path):
+    # need the pixel size for the raster so we can create an appropriate
+    # kernel for convolution
+    threat_pixel_size = pygeoprocessing.get_raster_info(
+        dist_raster_path)['pixel_size']
 
-    cols_per_block, rows_per_block = kernel_band.GetBlockSize()
+    # convert max distance (given in KM) to meters
+    max_dist_m = max_dist * 1000.0
 
-    n_cols = kernel_dataset.RasterXSize
-    n_rows = kernel_dataset.RasterYSize
+    # convert max distance from meters to the number of pixels that
+    # represents on the raster
+    max_dist_pixel = max_dist_m / abs(threat_pixel_size[0])
+    LOGGER.debug(f'Max distance in pixels: {max_dist_pixel}')
 
-    n_col_blocks = int(math.ceil(n_cols / float(cols_per_block)))
-    n_row_blocks = int(math.ceil(n_rows / float(rows_per_block)))
+    def exp_op(dist):
+        out_array = numpy.empty_like(dist)
 
-    integration = 0.0
-    for row_block_index in range(n_row_blocks):
-        row_offset = row_block_index * rows_per_block
-        row_block_width = n_rows - row_offset
-        if row_block_width > rows_per_block:
-            row_block_width = rows_per_block
+        out_array = numpy.where(
+            dist > max_dist_pixel, 0.0,
+            numpy.exp((-dist * scalar) / max_dist_pixel))
+        return out_array
 
-        for col_block_index in range(n_col_blocks):
-            col_offset = col_block_index * cols_per_block
-            col_block_width = n_cols - col_offset
-            if col_block_width > cols_per_block:
-                col_block_width = cols_per_block
+    pygeoprocessing.raster_calculator(
+        [(dist_raster_path, 1)], exp_op, target_path, gdal.GDT_Float32, -1)
 
-            # Numpy creates index rasters as ints by default, which sometimes
-            # creates problems on 32-bit builds when we try to add Int32
-            # matrices to float64 matrices.
-            row_indices, col_indices = numpy.indices((row_block_width,
-                                                      col_block_width),
-                                                     dtype=float)
 
-            row_indices += float(row_offset - max_distance)
-            col_indices += float(col_offset - max_distance)
+def dist_edt_linear(dist_raster_path, max_dist, target_path):
+    # need the pixel size for the raster so we can create an appropriate
+    # kernel for convolution
+    threat_pixel_size = pygeoprocessing.get_raster_info(
+        dist_raster_path)['pixel_size']
 
-            kernel_index_distances = numpy.hypot(
-                row_indices, col_indices)
-            kernel = numpy.where(
-                kernel_index_distances > max_distance, 0.0,
-                numpy.exp(
-                    (-kernel_index_distances * 2.99) / expected_distance))
+    # convert max distance (given in KM) to meters
+    max_dist_m = max_dist * 1000.0
 
-            kernel_band.WriteArray(kernel, xoff=col_offset,
-                                   yoff=row_offset)
+    # convert max distance from meters to the number of pixels that
+    # represents on the raster
+    max_dist_pixel = max_dist_m / abs(threat_pixel_size[0])
+    LOGGER.debug(f'Max distance in pixels: {max_dist_pixel}')
 
-    # Need to flush the kernel's cache to disk before opening up a new Dataset
-    # object in interblocks()
-    kernel_band.FlushCache()
-    kernel_dataset.FlushCache()
+    def linear_op(dist):
+        out_array = numpy.empty_like(dist)
 
-    for block_data in pygeoprocessing.iterblocks(
-            (kernel_filepath, 1), offset_only=True):
-        kernel_block = kernel_band.ReadAsArray(**block_data)
-        kernel_block /= integration
-        kernel_band.WriteArray(kernel_block, xoff=block_data['xoff'],
-                               yoff=block_data['yoff'])
+        out_array = numpy.where(
+            dist > max_dist_pixel, 0.0,
+            (max_dist_pixel - dist) / max_dist_pixel)
+        return out_array
 
-    kernel_band.FlushCache()
-    kernel_dataset.FlushCache()
-    kernel_band = None
-    kernel_dataset = None
+    pygeoprocessing.raster_calculator(
+        [(dist_raster_path, 1)], linear_op, target_path, gdal.GDT_Float32, -1)