Repository: natcap/invest
Pull request #730 (2021-11-20 00:13:27)

Urban Nature Access - Kernels and Convolutions
This PR is centered around kernels and convolutions, and not strictly limited to the Urban Nature Access model.

Changes to UNA in this PR include:

* Added 2 new kernels:
  * "Instantaneous"/"dichotomous" kernel - what we've called a 'disk' in other models.  It's a boolean kernel within a circular search radius.
  * "density" kernel - looks like a sphere if you were to look at it head-on.
* Added a new bit of functionality from the model, the "2-Step Floating Catchment Area" calculation, implemented as 2 convolutions with non-normalized kernels and a raster_calculator call.
* The LULC is now warped to ensure that we're using square pixels (since we're convolving based on a search radius, square pixels are needed)

Changes outside of UNA include:

* ``utils`` now has a gaussian decay kernel function.
* I've refactored GLOBIO and the Scenario Generator:Proximity to use the new gaussian decay kernel function.
* @esoth I updated a number of membership  assertions in the `ARGS_SPEC` test to make it much easier to debug when a key is missing.  Previously `self.assertTrue(a in b)` would just raise an exception saying `"False is not True"`.  Now, using `self.assertIn`, the exception will say something like `"10 could not be found in {'foo': 1}"` where `a=10` and `b={'foo: 1}`, making debugging errors a whole lot easier.


Patch:
diff --git a/src/natcap/invest/globio.py b/src/natcap/invest/globio.py
index f884b1637..0060e9156 100644
--- a/src/natcap/invest/globio.py
+++ b/src/natcap/invest/globio.py
@@ -1,22 +1,21 @@
 """GLOBIO InVEST Model."""
-import os
-import logging
 import collections
+import logging
+import os
 import tempfile
 
-from osgeo import gdal
-from osgeo import ogr
-from osgeo import osr
 import numpy
 import pygeoprocessing
 import taskgraph
+from osgeo import gdal
+from osgeo import ogr
+from osgeo import osr
 
-from . import utils
+from . import MODEL_METADATA
 from . import spec_utils
-from .spec_utils import u
+from . import utils
 from . import validation
-from . import MODEL_METADATA
-
+from .spec_utils import u
 
 LOGGER = logging.getLogger(__name__)
 
@@ -268,8 +267,8 @@ def execute(args):
     gaussian_kernel_path = os.path.join(
         tmp_dir, 'gaussian_kernel%s.tif' % file_suffix)
     make_gaussian_kernel_task = task_graph.add_task(
-        func=make_gaussian_kernel_path,
-        args=(SIGMA, gaussian_kernel_path),
+        func=utils.gaussian_decay_kernel_raster,
+        args=(SIGMA, gaussian_kernel_path, 5.0),
         target_path_list=[gaussian_kernel_path],
         task_name='gaussian_kernel')
 
@@ -681,45 +680,6 @@ def _msa_calculation(
         msa_op, msa_path, gdal.GDT_Float32, msa_nodata)
 
 
-def make_gaussian_kernel_path(sigma, kernel_path):
-    """Create a gaussian kernel raster."""
-    max_distance = sigma * 5
-    kernel_size = int(numpy.round(max_distance * 2 + 1))
-
-    driver = gdal.GetDriverByName('GTiff')
-    kernel_dataset = driver.Create(
-        kernel_path.encode('utf-8'), kernel_size, kernel_size, 1,
-        gdal.GDT_Float32, options=['BIGTIFF=IF_SAFER'])
-
-    # 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([444720, 30, 0, 3751320, 0, -30])
-    srs = osr.SpatialReference()
-    srs.SetUTM(11, 1)
-    srs.SetWellKnownGeogCS('NAD27')
-    kernel_dataset.SetProjection(srs.ExportToWkt())
-
-    kernel_band = kernel_dataset.GetRasterBand(1)
-    kernel_band.SetNoDataValue(-9999)
-
-    col_index = numpy.array(range(kernel_size))
-    integration = 0.0
-    for row_index in range(kernel_size):
-        kernel = numpy.exp(
-            -((row_index - max_distance)**2 +
-              (col_index - max_distance) ** 2)/(2.0*sigma**2)).reshape(
-                  1, kernel_size)
-
-        integration += numpy.sum(kernel)
-        kernel_band.WriteArray(kernel, xoff=0, yoff=row_index)
-
-    for row_index in range(kernel_size):
-        kernel_row = kernel_band.ReadAsArray(
-            xoff=0, yoff=row_index, win_xsize=kernel_size, win_ysize=1)
-        kernel_row /= integration
-        kernel_band.WriteArray(kernel_row, 0, row_index)
-
-
 def load_msa_parameter_table(
         msa_parameter_table_filename, intensification_fraction):
     """Load parameter table to a dict that to define the MSA ranges.
diff --git a/src/natcap/invest/scenario_gen_proximity.py b/src/natcap/invest/scenario_gen_proximity.py
index c0c74691d..daaf0cca9 100644
--- a/src/natcap/invest/scenario_gen_proximity.py
+++ b/src/natcap/invest/scenario_gen_proximity.py
@@ -1,26 +1,26 @@
 """Scenario Generation: Proximity Based."""
+import collections
+import heapq
+import logging
 import math
-import shutil
 import os
-import logging
-import tempfile
+import shutil
 import struct
-import heapq
+import tempfile
 import time
-import collections
 
 import numpy
-from osgeo import osr
-from osgeo import gdal
-import scipy
 import pygeoprocessing
+import scipy
 import taskgraph
+from osgeo import gdal
+from osgeo import osr
 
-from . import utils
+from . import MODEL_METADATA
 from . import spec_utils
-from .spec_utils import u
+from . import utils
 from . import validation
-from . import MODEL_METADATA
+from .spec_utils import u
 
 LOGGER = logging.getLogger(__name__)
 
@@ -369,7 +369,8 @@ def _convert_landscape(
     }
     # a sigma of 1.0 gives nice visual results to smooth pixel level artifacts
     # since a pixel is the 1.0 unit
-    _make_gaussian_kernel_path(1.0, tmp_file_registry['gaussian_kernel'])
+    utils.gaussian_decay_kernel_raster(
+        1.0, tmp_file_registry['gaussian_kernel'])
 
     # create the output raster first as a copy of the base landcover so it can
     # be looped on for each step
@@ -796,62 +797,6 @@ def _convert_by_score(
         stats_cache)
 
 
-def _make_gaussian_kernel_path(sigma, kernel_path):
-    """Create a 2D Gaussian kernel.
-
-    Args:
-        sigma (float): the sigma as in the classic Gaussian function
-        kernel_path (string): path to raster on disk to write the gaussian
-            kernel.
-
-    Returns:
-        None.
-
-    """
-    # going 3.0 times out from the sigma gives you over 99% of area under
-    # the guassian curve
-    max_distance = sigma * 3.0
-    kernel_size = int(numpy.round(max_distance * 2 + 1))
-
-    driver = gdal.GetDriverByName('GTiff')
-    kernel_dataset = driver.Create(
-        kernel_path.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)
-
-    col_index = numpy.array(range(kernel_size))
-    running_sum = 0.0
-    for row_index in range(kernel_size):
-        distance_kernel_row = numpy.sqrt(
-            (row_index - max_distance) ** 2 +
-            (col_index - max_distance) ** 2).reshape(1, kernel_size)
-        kernel = numpy.where(
-            distance_kernel_row > max_distance, 0.0,
-            (1 / (2.0 * numpy.pi * sigma ** 2) *
-             numpy.exp(-distance_kernel_row**2 / (2 * sigma ** 2))))
-        running_sum += numpy.sum(kernel)
-        kernel_band.WriteArray(kernel, xoff=0, yoff=row_index)
-
-    kernel_dataset.FlushCache()
-    for kernel_data, kernel_block in pygeoprocessing.iterblocks(
-            (kernel_path, 1)):
-        # divide by sum to normalize
-        kernel_block /= running_sum
-        kernel_band.WriteArray(
-            kernel_block, xoff=kernel_data['xoff'], yoff=kernel_data['yoff'])
-
-
 @validation.invest_validator
 def validate(args, limit_to=None):
     """Validate args to ensure they conform to `execute`'s contract.
diff --git a/src/natcap/invest/urban_nature_access.py b/src/natcap/invest/urban_nature_access.py
index 64cc4248e..e8a6f6618 100644
--- a/src/natcap/invest/urban_nature_access.py
+++ b/src/natcap/invest/urban_nature_access.py
@@ -1,4 +1,5 @@
 import logging
+import math
 import os
 import shutil
 import tempfile
@@ -17,7 +18,8 @@ from .spec_utils import u
 
 LOGGER = logging.getLogger(__name__)
 UINT32_NODATA = int(numpy.iinfo(numpy.uint32).max)
-FLOAT32_NODATA = float(numpy.finfo(numpy.float32).max)
+FLOAT32_NODATA = float(numpy.finfo(numpy.float32).min)
+BYTE_NODATA = 255
 ARGS_SPEC = {
     'model_name': MODEL_METADATA['urban_nature_access'].model_title,
     'pyname': MODEL_METADATA['urban_nature_access'].pyname,
@@ -97,14 +99,33 @@ ARGS_SPEC = {
             'units': u.m,
             'expression': "value > 0",
             'about': "",  # TODO, will know more about this when I implement.
+        },
+        'decay_function': {
+            'name': 'decay function',
+            'type': 'option_string',
+            'required': False,
+            'options': [
+                'dichotomy',
+                'exponential',
+                'gaussian',
+                'density',
+            ],
+            'about': '',  # TODO
         }
     }
 }
 
 
-_OUTPUT_BASE_FILES = {}
+_OUTPUT_BASE_FILES = {
+    'greenspace_supply': 'greenspace_supply.tif',
+}
 _INTERMEDIATE_BASE_FILES = {
-    'resampled_population': 'population_resampled.tif',
+    'aligned_population': 'aligned_population.tif',
+    'aligned_lulc': 'aligned_lulc.tif',
+    'greenspace_area': 'greenspace_area.tif',
+    'kernel': 'kernel.tif',
+    'greenspace_population_ratio': 'greenspace_population_ratio.tif',
+    'convolved_population': 'convolved_population.tif',
 }
 
 
@@ -139,6 +160,9 @@ def execute(args):
             number indicating the required greenspace, in m² per capita.
         args['search_radius'] (number): (required) A positive, nonzero number
             indicating the maximum distance that people travel for recreation.
+        args['kernel_type'] (string): (optional) The selected kernel type.
+            Must be one of the keys in ``KERNEL_TYPES``.  If not provided, the
+            ``'dichotomy'`` kernel will be used.
 
     Returns:
         ``None``
@@ -165,29 +189,292 @@ def execute(args):
         n_workers = -1  # Synchronous execution
     graph = taskgraph.TaskGraph(work_token_dir, n_workers)
 
+    kernel_types = {
+        'dichotomy': dichotomous_decay_kernel_raster,
+        # "exponential" is more consistent with other InVEST models'
+        # terminology.  "Power function" is used in the design doc.
+        'exponential': utils.exponential_decay_kernel_raster,
+        'gaussian': utils.gaussian_decay_kernel_raster,
+        'density': density_decay_kernel_raster,
+    }
+    # Since we have these keys defined in two places, I want to be super sure
+    # that the labels match.
+    assert sorted(kernel_types.keys()) == (
+        sorted(ARGS_SPEC['args']['decay_function']['options']))
+
+    if 'kernel_type' not in args:
+        kernel_type = 'dichotomy'
+    elif args['kernel_type'] not in kernel_types:
+        raise ValueError(
+            f'Kernel type "{args["kernel_type"]}" is not recognized. '
+            f"Must be one of {', '.join(kernel_types.keys())}")
+    else:
+        kernel_type = args['kernel_type']
+
     # Align the population raster to the LULC.
     lulc_raster_info = pygeoprocessing.get_raster_info(
         args['lulc_raster_path'])
 
+    squared_lulc_pixel_size = _square_off_pixels(args['lulc_raster_path'])
+    lulc_alignment_task = graph.add_task(
+        pygeoprocessing.warp_raster,
+        kwargs={
+            'base_raster_path': args['lulc_raster_path'],
+            'target_pixel_size': squared_lulc_pixel_size,
+            'target_bb': lulc_raster_info['bounding_box'],
+            'target_raster_path': file_registry['aligned_lulc'],
+            'resample_method': 'nearest',
+        },
+        target_path_list=[file_registry['aligned_lulc']],
+        task_name='Resample LULC to have square pixels'
+    )
+
     population_alignment_task = graph.add_task(
         _resample_population_raster,
         kwargs={
             'source_population_raster_path': args['population_raster_path'],
             'target_population_raster_path': file_registry[
-                'resampled_population'],
-            'lulc_pixel_size': lulc_raster_info['pixel_size'],
+                'aligned_population'],
+            'lulc_pixel_size': squared_lulc_pixel_size,
             'lulc_bb': lulc_raster_info['bounding_box'],
             'lulc_projection_wkt': lulc_raster_info['projection_wkt'],
             'working_dir': intermediate_dir,
         },
-        target_path_list=[file_registry['resampled_population']],
+        target_path_list=[file_registry['aligned_population']],
         task_name='Resample population to LULC resolution')
 
+    greenspace_lulc_lookup = utils.build_lookup_from_csv(
+        args['lulc_attribute_table'], 'lucode')
+    squared_pixel_area = abs(numpy.multiply(*squared_lulc_pixel_size))
+    greenspace_area_map = {
+        lucode: int(attributes['greenspace']) * squared_pixel_area
+        for lucode, attributes in greenspace_lulc_lookup.items()
+    }
+    greenspace_area_map[lulc_raster_info['nodata'][0]] = FLOAT32_NODATA
+    greenspace_reclassification_task = graph.add_task(
+        utils.reclassify_raster,
+        kwargs={
+            'raster_path_band': (file_registry['aligned_lulc'], 1),
+            'value_map': greenspace_area_map,
+            'target_raster_path': file_registry['greenspace_area'],
+            'target_datatype': gdal.GDT_Float32,
+            'target_nodata': FLOAT32_NODATA,
+            'error_details': {
+                'raster_name': ARGS_SPEC['args']['lulc_raster_path']['name'],
+                'column_name': 'greenspace',
+                'table_name': ARGS_SPEC['args'][
+                    'lulc_attribute_table']['name'],
+            },
+        },
+        target_path_list=[file_registry['greenspace_area']],
+        task_name='Identify the area of greenspace in pixels',
+        dependent_task_list=[lulc_alignment_task]
+    )
+
+    search_radius_in_pixels = abs(
+        float(args['search_radius']) / squared_lulc_pixel_size[0])
+    kernel_creation_task = graph.add_task(
+        # All kernel creation types have the same function signature
+        kernel_types[kernel_type],
+        args=(search_radius_in_pixels, file_registry['kernel']),
+        kwargs={'normalize': False},  # Model math calls for un-normalized
+        task_name=f'2SFCA - Create {kernel_type} kernel',
+        target_path_list=[file_registry['kernel']]
+    )
+
+    convolved_population_task = graph.add_task(
+        _convolve_and_set_lower_bounds_for_population,
+        kwargs={
+            'signal_path_band': (file_registry['aligned_population'], 1),
+            'kernel_path_band': (file_registry['kernel'], 1),
+            'target_path': file_registry['convolved_population'],
+            'working_dir': intermediate_dir,
+        },
+        task_name='2SFCA - Convolve population',
+        target_path_list=[file_registry['convolved_population']],
+        dependent_task_list=[
+            kernel_creation_task,
+            population_alignment_task,
+        ])
+
+    greenspace_nodata = pygeoprocessing.get_raster_info(
+        file_registry['greenspace_area'])['nodata'][0]
+    population_nodata = pygeoprocessing.get_raster_info(
+        file_registry['convolved_population'])['nodata'][0]
+    greenspace_population_ratio_task = graph.add_task(
+        pygeoprocessing.raster_calculator,
+        kwargs={
+            'base_raster_path_band_const_list': [
+                (file_registry['greenspace_area'], 1),
+                (file_registry['convolved_population'], 1),
+                (greenspace_nodata, 'raw'),
+                (population_nodata, 'raw')],
+            'local_op': _greenspace_population_ratio,
+            'target_raster_path': file_registry['greenspace_population_ratio'],
+            'datatype_target': gdal.GDT_Float32,
+            'nodata_target': FLOAT32_NODATA
+        },
+        task_name='2SFCA: Calculate R_j greenspace/population ratio',
+        target_path_list=[file_registry['greenspace_population_ratio']],
+        dependent_task_list=[
+            convolved_population_task,
+            greenspace_reclassification_task,
+        ])
+
+    convolved_greenspace_population_ratio_task = graph.add_task(
+        pygeoprocessing.convolve_2d,
+        kwargs={
+            'signal_path_band': (
+                file_registry['greenspace_population_ratio'], 1),
+            'kernel_path_band': (file_registry['kernel'], 1),
+            'target_path': file_registry['greenspace_supply'],
+            'working_dir': intermediate_dir,
+        },
+        task_name='2SFCA - greenspace supply',
+        target_path_list=[file_registry['greenspace_supply']],
+        dependent_task_list=[
+            kernel_creation_task,
+            greenspace_population_ratio_task,
+        ])
+
     graph.close()
     graph.join()
     LOGGER.info('Finished Urban Nature Access Model')
 
 
+def _greenspace_population_ratio(
+        greenspace_area, convolved_population, greenspace_nodata,
+        population_nodata):
+    """Calculate the greenspace-population ratio R_j.
+
+    Args:
+        greenspace_area (numpy.array): A numpy array representing the area of
+            greenspace in the pixel.  Pixel values will be ``0`` if there is no
+            greenspace.  Pixel values may also match ``greenspace_nodata``.
+        convolved_population (numpy.array): A numpy array where each pixel
+            represents the total number of people within a search radius of
+            each pixel, perhaps adjusted by a search kernel.
+        greenspace_nodata (int or float): The nodata value of the
+            ``greenspace_area`` matrix.
+        population_nodata (int or float): The nodata value of the
+            ``convolved_population`` matrix.
+
+    Returns:
+        A numpy array with the ratio ``R_j`` representing the
+        greenspace-population ratio with the following constraints:
+
+            * ``convolved_population`` pixels that are numerically close to
+              ``0`` are snapped to ``0`` to avoid unrealistically small
+              denominators in the final ratio.
+            * Any non-greenspace pixels will have a value of ``0.0`` in the
+              output matrix.
+    """
+    # ASSUMPTION: population nodata value is not close to 0.
+    #  Shouldn't be if we're coming from convolution.
+    out_array = numpy.full(
+        greenspace_area.shape, FLOAT32_NODATA, dtype=numpy.float32)
+
+    # Small negative values should already have been filtered out in another
+    # function after the convolution.
+    # This avoids divide-by-zero errors when taking the ratio.
+    valid_pixels = (convolved_population > 0)
+
+    # R_j is a ratio only calculated for the greenspace pixels.
+    greenspace_pixels = ~numpy.isclose(greenspace_area, 0.0)
+    valid_pixels &= greenspace_pixels
+    if population_nodata is not None:
+        valid_pixels &= ~numpy.isclose(
+            convolved_population, population_nodata)
+
+    if greenspace_nodata is not None:
+        valid_pixels &= ~numpy.isclose(greenspace_area, greenspace_nodata)
+
+    # If the population in the search radius is numerically 0, the model
+    # specifies that the ratio should be set to the greenspace area.
+    population_close_to_zero = numpy.isclose(convolved_population, 0.0)
+    out_array[population_close_to_zero] = (
+        greenspace_pixels[population_close_to_zero])
+    out_array[~greenspace_pixels] = 0.0  # set non-greenspace non-nodata to 0
+    out_array[valid_pixels] = (
+        greenspace_area[valid_pixels] / convolved_population[valid_pixels])
+
+    return out_array
+
+
+def _convolve_and_set_lower_bounds_for_population(
+        signal_path_band, kernel_path_band, target_path, working_dir):
+    """Convolve a raster and set all values below 0 to 0.
+
+    Args:
+        signal_path_band (tuple): A 2-tuple of (signal_raster_path, band_index)
+            to use as the signal raster in the convolution.
+        kernel_path_band (tuple): A 2-tuple of (kernel_raster_path, band_index)
+            to use as the kernel raster in the convolution.
+        target_path (string): Where the target raster should be written.
+        working_dir (string): The working directory that
+            ``pygeoprocessing.convolve_2d`` may use for its intermediate files.
+
+    Returns:
+        ``None``
+    """
+    pygeoprocessing.convolve_2d(
+        signal_path_band=signal_path_band,
+        kernel_path_band=kernel_path_band,
+        target_path=target_path,
+        working_dir=working_dir)
+
+    # Sometimes there are negative values that should have been clamped to 0 in
+    # the convolution but weren't, so let's clamp them to avoid support issues
+    # later on.
+    target_raster = gdal.OpenEx(target_path, gdal.GA_Update)
+    target_band = target_raster.GetRasterBand(1)
+    target_nodata = target_band.GetNoDataValue()
+    for block_data in pygeoprocessing.iterblocks(
+            (target_path, 1), offset_only=True):
+        block = target_band.ReadAsArray(**block_data)
+        valid_pixels = slice(None)
+        if target_nodata is not None:
+            valid_pixels = ~numpy.isclose(block, target_nodata)
+        block[(block < 0.0) & valid_pixels] = 0.0
+        target_band.WriteArray(
+            block, xoff=block_data['xoff'], yoff=block_data['yoff'])
+
+    target_band = None
+    target_raster = None
+
+
+def _square_off_pixels(raster_path):
+    """Create square pixels from the provided raster.
+
+    The pixel dimensions produced will respect the sign of the original pixel
+    dimensions and will be the mean of the absolute source pixel dimensions.
+
+    Args:
+        raster_path (string): The path to a raster on disk.
+
+    Returns:
+        A 2-tuple of ``(pixel_width, pixel_height)``, in projected units.
+    """
+    raster_info = pygeoprocessing.get_raster_info(raster_path)
+    pixel_width, pixel_height = raster_info['pixel_size']
+
+    if abs(pixel_width) == abs(pixel_height):
+        return (pixel_width, pixel_height)
+
+    pixel_tuple = ()
+    average_absolute_size = (abs(pixel_width) + abs(pixel_height)) / 2
+    for pixel_dimension_size in (pixel_width, pixel_height):
+        # This loop allows either or both pixel dimension(s) to be negative
+        sign_factor = 1
+        if pixel_dimension_size < 0:
+            sign_factor = -1
+
+        pixel_tuple += (average_absolute_size * sign_factor,)
+
+    return pixel_tuple
+
+
+# TODO: refactor this into raster_calculator and align_and_resize...
 def _resample_population_raster(
         source_population_raster_path, target_population_raster_path,
         lulc_pixel_size, lulc_bb, lulc_projection_wkt, working_dir):
@@ -316,6 +603,186 @@ def _resample_population_raster(
     shutil.rmtree(tmp_working_dir, ignore_errors=True)
 
 
+def dichotomous_decay_kernel_raster(expected_distance, kernel_filepath,
+        normalize=False):
+    """Create a raster-based, discontinuous decay kernel based on a dichotomy.
+
+    This kernel has a value of ``1`` for all pixels within
+    ``expected_distance`` from the center of the kernel.  All values outside of
+    this distance are ``0``.
+
+    Args:
+        expected_distance (int or float): The distance (in pixels) after which
+            the kernel becomes 0.
+        kernel_filepath (string): The string path on disk to where this kernel
+            should be stored.
+        normalize=False (bool): Whether to divide the kernel values by the sum
+            of all values in the kernel.
+
+    Returns:
+        ``None``
+    """
+    pixel_radius = math.ceil(expected_distance)
+    kernel_size = pixel_radius * 2 + 1  # allow for a center pixel
+    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_nodata = FLOAT32_NODATA
+    kernel_band.SetNoDataValue(kernel_nodata)
+
+    kernel_band = None
+    kernel_dataset = None
+
+    kernel_raster = gdal.OpenEx(kernel_filepath, gdal.GA_Update)
+    kernel_band = kernel_raster.GetRasterBand(1)
+    band_x_size = kernel_band.XSize
+    band_y_size = kernel_band.YSize
+    running_sum = 0.0
+    for block_data in pygeoprocessing.iterblocks(
+            (kernel_filepath, 1), offset_only=True):
+        array_xmin = block_data['xoff'] - pixel_radius
+        array_xmax = min(
+            array_xmin + block_data['win_xsize'],
+            band_x_size - pixel_radius)
+        array_ymin = block_data['yoff'] - pixel_radius
+        array_ymax = min(
+            array_ymin + block_data['win_ysize'],
+            band_y_size - pixel_radius)
+
+        pixel_dist_from_center = numpy.hypot(
+            *numpy.mgrid[
+                array_ymin:array_ymax,
+                array_xmin:array_xmax])
+        search_kernel = numpy.array(
+            pixel_dist_from_center <= expected_distance, dtype=numpy.uint8)
+        running_sum += search_kernel.sum()
+        kernel_band.WriteArray(
+            search_kernel,
+            yoff=block_data['yoff'],
+            xoff=block_data['xoff'])
+
+    kernel_raster.FlushCache()
+    kernel_band = None
+    kernel_raster = None
+
+    if normalize:
+        kernel_raster = gdal.OpenEx(kernel_filepath, gdal.GA_Update)
+        kernel_band = kernel_raster.GetRasterBand(1)
+        for block_data, kernel_block in pygeoprocessing.iterblocks(
+                (kernel_filepath, 1)):
+            # divide by sum to normalize
+            kernel_block /= running_sum
+            kernel_band.WriteArray(
+                kernel_block, xoff=block_data['xoff'], yoff=block_data['yoff'])
+
+        kernel_raster.FlushCache()
+        kernel_band = None
+        kernel_raster = None
+
+
+def density_decay_kernel_raster(expected_distance, kernel_filepath,
+        normalize=False):
+    """Create a raster-based density decay kernel.
+
+    Args:
+        expected_distance (int or float): The distance (in pixels) after which
+            the kernel becomes 0.
+        kernel_filepath (string): The string path on disk to where this kernel
+            should be stored.
+        normalize=False (bool): Whether to divide the kernel values by the sum
+            of all values in the kernel.
+
+    Returns:
+        ``None``
+    """
+    pixel_radius = math.ceil(expected_distance)
+    kernel_size = pixel_radius * 2 + 1  # allow for a center pixel
+    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_nodata = float(numpy.finfo(numpy.float32).min)
+    kernel_band.SetNoDataValue(kernel_nodata)
+
+    kernel_band = None
+    kernel_dataset = None
+
+    kernel_raster = gdal.OpenEx(kernel_filepath, gdal.GA_Update)
+    kernel_band = kernel_raster.GetRasterBand(1)
+    band_x_size = kernel_band.XSize
+    band_y_size = kernel_band.YSize
+    running_sum = 0.0
+    for block_data in pygeoprocessing.iterblocks(
+            (kernel_filepath, 1), offset_only=True):
+        array_xmin = block_data['xoff'] - pixel_radius
+        array_xmax = min(
+            array_xmin + block_data['win_xsize'],
+            band_x_size - pixel_radius)
+        array_ymin = block_data['yoff'] - pixel_radius
+        array_ymax = min(
+            array_ymin + block_data['win_ysize'],
+            band_y_size - pixel_radius)
+
+        pixel_dist_from_center = numpy.hypot(
+            *numpy.mgrid[
+                array_ymin:array_ymax,
+                array_xmin:array_xmax])
+
+        density = numpy.zeros(
+            pixel_dist_from_center.shape, dtype=numpy.float32)
+        pixels_in_radius = (pixel_dist_from_center <= expected_distance)
+        density[pixels_in_radius] = (
+            0.75 * (1 - (pixel_dist_from_center[
+                pixels_in_radius] / expected_distance) ** 2))
+        running_sum += density.sum()
+
+        kernel_band.WriteArray(
+            density,
+            yoff=block_data['yoff'],
+            xoff=block_data['xoff'])
+
+    kernel_raster.FlushCache()
+    kernel_band = None
+    kernel_raster = None
+
+    if normalize:
+        kernel_raster = gdal.OpenEx(kernel_filepath, gdal.GA_Update)
+        kernel_band = kernel_raster.GetRasterBand(1)
+        for block_data, kernel_block in pygeoprocessing.iterblocks(
+                (kernel_filepath, 1)):
+            # divide by sum to normalize
+            kernel_block /= running_sum
+            kernel_band.WriteArray(
+                kernel_block, xoff=block_data['xoff'], yoff=block_data['yoff'])
+
+        kernel_raster.FlushCache()
+        kernel_band = None
+        kernel_raster = None
+
+
 def validate(args, limit_to=None):
     return validation.validate(
         args, ARGS_SPEC['args'], ARGS_SPEC['args_with_spatial_overlap'])
diff --git a/src/natcap/invest/utils.py b/src/natcap/invest/utils.py
index 4dcd08d23..70700beea 100644
--- a/src/natcap/invest/utils.py
+++ b/src/natcap/invest/utils.py
@@ -1,21 +1,20 @@
 """InVEST specific code utils."""
 import codecs
-import math
-import os
 import contextlib
 import logging
-import tempfile
+import math
+import os
 import shutil
-from datetime import datetime
+import tempfile
 import time
+from datetime import datetime
 
-import pandas
 import numpy
-from shapely.wkt import loads
+import pandas
+import pygeoprocessing
 from osgeo import gdal
 from osgeo import osr
-import pygeoprocessing
-
+from shapely.wkt import loads
 
 LOGGER = logging.getLogger(__name__)
 _OSGEO_LOGGER = logging.getLogger('osgeo')
@@ -326,7 +325,8 @@ def make_suffix_string(args, suffix_key):
     return file_suffix
 
 
-def exponential_decay_kernel_raster(expected_distance, kernel_filepath):
+def exponential_decay_kernel_raster(expected_distance, kernel_filepath,
+        normalize=True):
     """Create a raster-based exponential decay kernel.
 
     The raster created will be a tiled GeoTiff, with 256x256 memory blocks.
@@ -338,6 +338,8 @@ def exponential_decay_kernel_raster(expected_distance, kernel_filepath):
         kernel_filepath (string): The path to the file on disk where this
             kernel should be stored.  If this file exists, it will be
             overwritten.
+        normalize=True (bool): Whether to divide the kernel values by the sum
+            of all values in the kernel.
 
     Returns:
         None
@@ -407,19 +409,108 @@ def exponential_decay_kernel_raster(expected_distance, kernel_filepath):
     # object in interblocks()
     kernel_band.FlushCache()
     kernel_dataset.FlushCache()
+    kernel_band = None
+    kernel_dataset = None
 
-    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'])
+    if normalize:
+        kernel_dataset = gdal.OpenEx(kernel_filepath, gdal.GA_Update)
+        kernel_band = kernel_dataset.GetRasterBand(1)
+        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'])
+
+        kernel_band.FlushCache()
+        kernel_dataset.FlushCache()
+        kernel_band = None
+        kernel_dataset = None
+
+
+def gaussian_decay_kernel_raster(
+        sigma, kernel_filepath, n_std_dev=3.0, normalize=True):
+    """Create a raster-based gaussian decay kernel.
+
+    The raster will be a tiled GeoTIFF, with 256x256 memory blocks.
+
+    While the ``sigma`` parameter represents the width of a standard deviation
+    in pixels, the ``n_std_dev`` parameter defines how many standard deviations
+    should be included in the resulting kernel.  The resulting kernel raster
+    will be square in shape, with a width of ``(sigma * n_std_dev * 2) + 1``
+    pixels.
+
+    Args:
+        sigma (int or float): The distance (in pixels) of the standard
+            deviation from the center of the raster.
+        kernel_filepath (string): The path to the file on disk where this
+            kernel should be stored. If a file exists at this path, it will be
+            overwritten.
+        n_std_dev=3.0 (int or float): The number of times sigma should be
+            multiplied in order to get the pixel radius of the resulting
+            kernel.  The default of 3 standard deviations will cover 99.7% of
+            the area under the gaussian curve.
+        normalize=True (bool): Whether to divide the kernel values by the sum
+            of all values in the kernel.
+
+    Returns:
+        ``None``
+    """
+    # going 3.0 times out from the sigma gives you over 99% of area under
+    # the gaussian curve
+    max_distance = sigma * n_std_dev
+    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_nodata = -9999
+    kernel_band.SetNoDataValue(kernel_nodata)
+
+    col_index = numpy.array(range(kernel_size))
+    running_sum = 0.0
+    for row_index in range(kernel_size):
+        distance_kernel_row = numpy.sqrt(
+            (row_index - max_distance) ** 2 +
+            (col_index - max_distance) ** 2).reshape(1, kernel_size)
+        kernel = numpy.where(
+            distance_kernel_row > max_distance, 0.0,
+            (1 / (2.0 * numpy.pi * sigma ** 2) *
+             numpy.exp(-distance_kernel_row**2 / (2 * sigma ** 2))))
+        running_sum += numpy.sum(kernel)
+        kernel_band.WriteArray(kernel, xoff=0, yoff=row_index)
 
-    kernel_band.FlushCache()
     kernel_dataset.FlushCache()
     kernel_band = None
     kernel_dataset = None
 
+    if normalize:
+        kernel_dataset = gdal.OpenEx(kernel_filepath, gdal.GA_Update)
+        kernel_band = kernel_dataset.GetRasterBand(1)
+        for kernel_data, kernel_block in pygeoprocessing.iterblocks(
+                (kernel_filepath, 1)):
+            # divide by sum to normalize
+            kernel_block /= running_sum
+            kernel_band.WriteArray(
+                kernel_block, xoff=kernel_data['xoff'],
+                yoff=kernel_data['yoff'])
+
+        kernel_dataset.FlushCache()
+        kernel_band = None
+        kernel_dataset = None
+
 
 def build_file_registry(base_file_path_list, file_suffix):
     """Combine file suffixes with key names, base filenames, and directories.


Test patch:
diff --git a/tests/test_args_specs.py b/tests/test_args_specs.py
index 8152c8f48..3eff75870 100644
--- a/tests/test_args_specs.py
+++ b/tests/test_args_specs.py
@@ -2,9 +2,10 @@ import importlib
 import re
 import unittest
 
-from natcap.invest import cli, MODEL_METADATA
 import pint
 
+from natcap.invest import cli
+from natcap.invest import MODEL_METADATA
 
 valid_nested_types = {
     'all': {
@@ -99,19 +100,19 @@ class ValidateArgsSpecs(unittest.TestCase):
         """
         with self.subTest(nested_arg_name=name):
             for attr in required_attrs:
-                self.assertTrue(attr in arg)
+                self.assertIn(attr, arg)
 
             # arg['type'] can be either a string or a set of strings
             types = arg['type'] if isinstance(
                 arg['type'], set) else [arg['type']]
             attrs = set(arg.keys())
             for t in types:
-                self.assertTrue(t in valid_types)
+                self.assertIn(t, valid_types)
 
                 if t == 'option_string':
                     # option_string type should have an options property that
                     # describes the valid options
-                    self.assertTrue('options' in arg)
+                    self.assertIn('options', arg)
                     # May be a list or dict because some option sets are self
                     # explanatory and others need a description
                     self.assertTrue(isinstance(arg['options'], dict) or
@@ -139,7 +140,7 @@ class ValidateArgsSpecs(unittest.TestCase):
 
                 elif t == 'number':
                     # number type should have a units property
-                    self.assertTrue('units' in arg)
+                    self.assertIn('units', arg)
                     # Undefined units should use the custom u.none unit
                     self.assertTrue(isinstance(arg['units'], pint.Unit))
                     attrs.remove('units')
@@ -157,7 +158,7 @@ class ValidateArgsSpecs(unittest.TestCase):
                 elif t == 'raster':
                     # raster type should have a bands property that maps each band
                     # index to a nested type dictionary describing the band's data
-                    self.assertTrue('bands' in arg)
+                    self.assertIn('bands', arg)
                     self.assertTrue(isinstance(arg['bands'], dict))
                     for band in arg['bands']:
                         self.assertTrue(isinstance(band, int))
@@ -187,7 +188,7 @@ class ValidateArgsSpecs(unittest.TestCase):
                     # - a fields property that maps each field header to a nested
                     #   type dictionary describing the data in that field
                     # - a geometries property: the set of valid geometry types
-                    self.assertTrue('fields' in arg)
+                    self.assertIn('fields', arg)
                     self.assertTrue(isinstance(arg['fields'], dict))
                     for field in arg['fields']:
                         self.assertTrue(isinstance(field, str))
@@ -196,7 +197,7 @@ class ValidateArgsSpecs(unittest.TestCase):
                             f'{name}.fields.{field}',
                             valid_types=valid_nested_types['vector'])
 
-                    self.assertTrue('geometries' in arg)
+                    self.assertIn('geometries', arg)
                     self.assertTrue(isinstance(arg['geometries'], set))
 
                     attrs.remove('fields')
@@ -252,7 +253,7 @@ class ValidateArgsSpecs(unittest.TestCase):
                     # directory type should have a contents property that maps each
                     # expected path name/pattern within the directory to a nested
                     # type dictionary describing the data at that filepath
-                    self.assertTrue('contents' in arg)
+                    self.assertIn('contents', arg)
                     self.assertTrue(isinstance(arg['contents'], dict))
                     for path in arg['contents']:
                         self.assertTrue(isinstance(path, str))
diff --git a/tests/test_urban_nature_access.py b/tests/test_urban_nature_access.py
index 99dede28c..dc57a101e 100644
--- a/tests/test_urban_nature_access.py
+++ b/tests/test_urban_nature_access.py
@@ -1,9 +1,11 @@
 # coding=UTF-8
 """Tests for the Urban Nature Access Model."""
+import math
 import os
 import random
 import shutil
 import tempfile
+import textwrap
 import unittest
 
 import numpy
@@ -86,6 +88,81 @@ class UNATests(unittest.TestCase):
                     population_array.sum(), resampled_population_array.sum(),
                     rtol=1e-3)
 
+    def test_dichotomous_decay_simple(self):
+        """UNA: Test dichotomous decay on a simple case."""
+        from natcap.invest import urban_nature_access
+
+        expected_distance = 5
+        kernel_filepath = os.path.join(self.workspace_dir, 'kernel.tif')
+
+        urban_nature_access.dichotomous_decay_kernel_raster(
+            expected_distance, kernel_filepath)
+
+        expected_array = numpy.array([
+            [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
+            [0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0],
+            [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
+            [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
+            [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
+            [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
+            [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
+            [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
+            [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
+            [0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0],
+            [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0]], dtype=numpy.uint8)
+
+        extracted_kernel_array = pygeoprocessing.raster_to_numpy_array(
+            kernel_filepath)
+        numpy.testing.assert_array_equal(
+            expected_array, extracted_kernel_array)
+
+    def test_dichotomous_decay_large(self):
+        """UNA: Test dichotomous decay on a very large pixel radius."""
+        from natcap.invest import urban_nature_access
+
+        # kernel with > 268 million pixels.  This is big enough to force my
+        # laptop to noticeably hang while swapping memory on an all in-memory
+        # implementation.
+        expected_distance = 2**13
+        kernel_filepath = os.path.join(self.workspace_dir, 'kernel.tif')
+
+        urban_nature_access.dichotomous_decay_kernel_raster(
+            expected_distance, kernel_filepath)
+
+        expected_shape = (expected_distance*2+1, expected_distance*2+1)
+        expected_n_1_pixels = math.pi*expected_distance**2
+
+        kernel_info = pygeoprocessing.get_raster_info(kernel_filepath)
+        n_1_pixels = 0
+        for _, block in pygeoprocessing.iterblocks((kernel_filepath, 1)):
+            n_1_pixels += numpy.count_nonzero(block)
+
+        # 210828417 is only a slight overestimate from the area of the circle
+        # at this radius: math.pi*expected_distance**2 = 210828714.13315654
+        numpy.testing.assert_allclose(
+            n_1_pixels, expected_n_1_pixels, rtol=1e-5)
+        self.assertEqual(kernel_info['raster_size'], expected_shape)
+
+    def test_density_decay(self):
+        """UNA: Test density decay."""
+        from natcap.invest import urban_nature_access
+
+        expected_distance = 200
+        kernel_filepath = os.path.join(self.workspace_dir, 'kernel.tif')
+
+        urban_nature_access.density_decay_kernel_raster(
+            expected_distance, kernel_filepath)
+
+        expected_shape = (expected_distance*2+1,) * 2
+        kernel_info = pygeoprocessing.get_raster_info(kernel_filepath)
+        kernel_array = pygeoprocessing.raster_to_numpy_array(kernel_filepath)
+        self.assertEqual(kernel_info['raster_size'], expected_shape)
+        numpy.testing.assert_allclose(
+            47123.867,  # obtained from manual inspection
+            kernel_array.sum())
+        self.assertEqual(0.75, kernel_array.max())
+        self.assertEqual(0, kernel_array.min())
+
     def test_model(self):
         """UNA: Run through the model."""
         from natcap.invest import urban_nature_access
@@ -96,6 +173,10 @@ class UNATests(unittest.TestCase):
             'population_raster_path': os.path.join(
                 self.workspace_dir, 'population.tif'),
             'lulc_raster_path': os.path.join(self.workspace_dir, 'lulc.tif'),
+            'lulc_attribute_table': os.path.join(
+                self.workspace_dir, 'lulc_attributes.csv'),
+            'decay_function': 'gaussian',
+            'search_radius': 100.0,  # meters
         }
 
         random.seed(-1)  # for our random number generation
@@ -128,7 +209,40 @@ class UNATests(unittest.TestCase):
             projection_wkt=population_wkt,
             target_path=args['lulc_raster_path'])
 
+        with open(args['lulc_attribute_table'], 'w') as attr_table:
+            attr_table.write(textwrap.dedent(
+                """lucode,greenspace
+                0,0
+                1,1
+                2,0
+                3,1
+                4,0
+                5,1
+                6,0
+                7,1
+                8,0
+                9,1"""))
+
         urban_nature_access.execute(args)
 
-        # TODO: Assertions will be added here later when I have something to
-        # test.
+        # Since we're doing a semi-manual alignment step, assert that the
+        # aligned LULC and population rasters have the same pixel sizes, origin
+        # and raster dimensions.
+        aligned_lulc_raster_info = pygeoprocessing.get_raster_info(
+            os.path.join(args['workspace_dir'], 'intermediate',
+                         f"aligned_lulc_{args['results_suffix']}.tif"))
+        aligned_population_raster_info = pygeoprocessing.get_raster_info(
+            os.path.join(args['workspace_dir'], 'intermediate',
+                         f"aligned_population_{args['results_suffix']}.tif"))
+        numpy.testing.assert_allclose(
+            aligned_lulc_raster_info['pixel_size'],
+            aligned_population_raster_info['pixel_size'])
+        numpy.testing.assert_allclose(
+            aligned_lulc_raster_info['raster_size'],
+            aligned_population_raster_info['raster_size'])
+        numpy.testing.assert_allclose(
+            aligned_lulc_raster_info['geotransform'],
+            aligned_population_raster_info['geotransform'])
+        numpy.testing.assert_allclose(
+            aligned_lulc_raster_info['bounding_box'],
+            aligned_population_raster_info['bounding_box'])
diff --git a/tests/test_utils.py b/tests/test_utils.py
index 5409e72fd..d98d04f2f 100644
--- a/tests/test_utils.py
+++ b/tests/test_utils.py
@@ -1,26 +1,27 @@
 """Module for testing the natcap.invest.utils module.."""
 import codecs
-import unittest
-import os
-import tempfile
-import shutil
+import glob
 import logging
 import logging.handlers
-import threading
-import warnings
+import os
+import queue
 import re
-import glob
+import shutil
+import tempfile
 import textwrap
-import queue
+import threading
+import unittest
+import warnings
 
 import numpy
+import numpy.testing
+import pygeoprocessing
+import scipy.ndimage
 from osgeo import gdal
 from osgeo import ogr
 from osgeo import osr
-from shapely.geometry import Polygon
 from shapely.geometry import Point
-
-import pygeoprocessing
+from shapely.geometry import Polygon
 
 
 class SuffixUtilsTests(unittest.TestCase):
@@ -211,6 +212,61 @@ class ExponentialDecayUtilsTests(unittest.TestCase):
         numpy.testing.assert_allclose(model_array, reg_array, atol=1e-6)
 
 
+class GaussianDecayUtilsTests(unittest.TestCase):
+    """Tests for natcap.invest.utils.gaussian_decay_kernel_raster."""
+
+    def setUp(self):
+        """Setup workspace."""
+        self.workspace_dir = tempfile.mkdtemp()
+
+    def tearDown(self):
+        """Delete workspace."""
+        shutil.rmtree(self.workspace_dir)
+
+    def test_gaussian_kernel_raster(self):
+        """Utils: test gaussian_decay_kernel_raster."""
+        from natcap.invest import utils
+
+        sigma = 10
+        # The kernel should have a radius of 3 standard deviations.
+        # We then double that for the diameter and add 1 to ensure there's a
+        # single pixel at the center.
+        kernel_dimension = sigma * 3 * 2 + 1
+
+        def gkern():
+            """
+            Creates a gaussian kernel with side length ``kernel_dimension``
+            and a sigma of ``sigma``.
+
+            This function adapted from the following stackoverflow answer:
+            https://stackoverflow.com/a/43346070/299084
+            """
+            ax = numpy.linspace(
+                -(kernel_dimension - 1) / 2.,
+                (kernel_dimension - 1) / 2.,
+                kernel_dimension)
+            gauss = numpy.exp(-0.5 * numpy.square(ax) / numpy.square(sigma))
+            kernel = numpy.outer(gauss, gauss)
+
+            dist_from_center = numpy.hypot(
+                *numpy.mgrid[
+                    -kernel_dimension/2:kernel_dimension/2,
+                    -kernel_dimension/2:kernel_dimension/2])
+            # The sigma*3 is the maximum radius from the center
+            # Anything greater than that distance should be set to 0 by the
+            # gaussian kernel creation function.
+            kernel[dist_from_center > (sigma * 3)] = 0.0
+            return kernel / numpy.sum(kernel)
+
+        expected_matrix = gkern()
+
+        kernel_filepath = os.path.join(self.workspace_dir, 'kernel.tif')
+        utils.gaussian_decay_kernel_raster(sigma, kernel_filepath)
+
+        disk_kernel = pygeoprocessing.raster_to_numpy_array(kernel_filepath)
+        numpy.testing.assert_allclose(expected_matrix, disk_kernel, atol=1e-4)
+
+
 class SandboxTempdirTests(unittest.TestCase):
     """Test Sandbox Tempdir."""