diff --git a/HISTORY.rst b/HISTORY.rst
index 522ab8477a..dca729789d 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -49,6 +49,15 @@ Unreleased Changes
       (`#1374 <https://github.com/natcap/invest/issues/1374>`_)
     * Datastack archives will now be correctly extracted
       (`#1308 <https://github.com/natcap/invest/issues/1308>`_)
+    * Updated to ``pygeoprocessing`` 2.4.2. This includes an update to
+      ``pygeoprocessing.zonal_statistics``, which is now more correct on certain
+      edge cases. Aggregated model results may change slightly.
+    * Removed the ``utils`` functions ``array_equals_nodata``,
+      ``exponential_decay_kernel_raster``, and ``gaussian_decay_kernel_raster``,
+      which were obsoleted by new ``pygeoprocessing`` features.
+* Pollination
+    * Replaced custom kernel implementation with ``pygeoprocessing.kernels``.
+      Convolution results may be slightly different (more accurate).
 * NDR
     * Fixing an issue where minor geometric issues in the watersheds input
       (such as a ring self-intersection) would raise an error in the model.
diff --git a/requirements.txt b/requirements.txt
index d9dc564d7e..5b2583b0fb 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -17,7 +17,7 @@ numpy>=1.11.0,!=1.16.0
 Rtree>=0.8.2,!=0.9.1
 shapely>=2.0.0
 scipy>=1.9.0
-pygeoprocessing==2.4.0
+pygeoprocessing>=2.4.2  # pip-only
 taskgraph>=0.11.0
 psutil>=5.6.6
 chardet>=3.0.4
diff --git a/src/natcap/invest/annual_water_yield.py b/src/natcap/invest/annual_water_yield.py
index 2bc37ac486..314ddb8153 100644
--- a/src/natcap/invest/annual_water_yield.py
+++ b/src/natcap/invest/annual_water_yield.py
@@ -733,12 +733,12 @@ def execute(args):
 
     LOGGER.info('Calculate PET from Ref Evap times Kc')
     calculate_pet_task = graph.add_task(
-        func=pygeoprocessing.raster_calculator,
-        args=([(eto_path, 1), (tmp_Kc_raster_path, 1),
-               (nodata_dict['eto'], 'raw'),
-               (nodata_dict['out_nodata'], 'raw')],
-              pet_op, tmp_pet_path, gdal.GDT_Float32,
-              nodata_dict['out_nodata']),
+        func=pygeoprocessing.raster_map,
+        kwargs=dict(
+            op=numpy.multiply,  # PET = ET0 * KC
+            rasters=[eto_path, tmp_Kc_raster_path],
+            target_path=tmp_pet_path,
+            target_nodata=nodata_dict['out_nodata']),
         target_path_list=[tmp_pet_path],
         dependent_task_list=[create_Kc_raster_task],
         task_name='calculate_pet')
@@ -764,12 +764,12 @@ def execute(args):
 
     LOGGER.info('Performing wyield operation')
     calculate_wyield_task = graph.add_task(
-        func=pygeoprocessing.raster_calculator,
-        args=([(fractp_path, 1), (precip_path, 1),
-               (nodata_dict['precip'], 'raw'),
-               (nodata_dict['out_nodata'], 'raw')],
-              wyield_op, wyield_path, gdal.GDT_Float32,
-              nodata_dict['out_nodata']),
+        func=pygeoprocessing.raster_map,
+        kwargs=dict(
+            op=wyield_op,
+            rasters=[fractp_path, precip_path],
+            target_path=wyield_path,
+            target_nodata=nodata_dict['out_nodata']),
         target_path_list=[wyield_path],
         dependent_task_list=[calculate_fractp_task, align_raster_stack_task],
         task_name='calculate_wyield')
@@ -777,11 +777,12 @@ def execute(args):
 
     LOGGER.debug('Performing aet operation')
     calculate_aet_task = graph.add_task(
-        func=pygeoprocessing.raster_calculator,
-        args=([(fractp_path, 1), (precip_path, 1),
-               (nodata_dict['precip'], 'raw'),
-               (nodata_dict['out_nodata'], 'raw')],
-              aet_op, aet_path, gdal.GDT_Float32, nodata_dict['out_nodata']),
+        func=pygeoprocessing.raster_map,
+        kwargs=dict(
+            op=numpy.multiply,  # AET = fractp * precip
+            rasters=[fractp_path, precip_path],
+            target_path=aet_path,
+            target_nodata=nodata_dict['out_nodata']),
         target_path_list=[aet_path],
         dependent_task_list=[
             calculate_fractp_task, create_veg_raster_task,
@@ -866,6 +867,10 @@ def execute(args):
     graph.join()
 
 
+# wyield equation to pass to raster_map
+def wyield_op(fractp, precip): return (1 - fractp) * precip
+
+
 def copy_vector(base_vector_path, target_vector_path):
     """Wrapper around CreateCopy that handles opening & closing the dataset.
 
@@ -977,55 +982,6 @@ def zonal_stats_tofile(base_vector_path, raster_path, target_stats_pickle):
         picklefile.write(pickle.dumps(ws_stats_dict))
 
 
-def aet_op(fractp, precip, precip_nodata, output_nodata):
-    """Compute actual evapotranspiration values.
-
-    Args:
-        fractp (numpy.ndarray float): fractp raster values.
-        precip (numpy.ndarray): precipitation raster values (mm).
-        precip_nodata (float): nodata value from the precip raster.
-        output_nodata (float): nodata value assigned to output of
-            raster_calculator.
-
-    Returns:
-        numpy.ndarray of actual evapotranspiration values (mm).
-
-    """
-    result = numpy.empty_like(fractp)
-    result[:] = output_nodata
-    # checking if fractp >= 0 because it's a value that's between 0 and 1
-    # and the nodata value is a large negative number.
-    valid_mask = fractp >= 0
-    if precip_nodata is not None:
-        valid_mask &= ~utils.array_equals_nodata(precip, precip_nodata)
-    result[valid_mask] = fractp[valid_mask] * precip[valid_mask]
-    return result
-
-
-def wyield_op(fractp, precip, precip_nodata, output_nodata):
-    """Calculate water yield.
-
-    Args:
-        fractp (numpy.ndarray float): fractp raster values.
-        precip (numpy.ndarray): precipitation raster values (mm).
-        precip_nodata (float): nodata value from the precip raster.
-        output_nodata (float): nodata value assigned to output of
-            raster_calculator.
-
-    Returns:
-        numpy.ndarray of water yield value (mm).
-
-    """
-    result = numpy.empty_like(fractp)
-    result[:] = output_nodata
-    # output_nodata is defined above, should never be None
-    valid_mask = ~utils.array_equals_nodata(fractp, output_nodata)
-    if precip_nodata is not None:
-        valid_mask &= ~utils.array_equals_nodata(precip, precip_nodata)
-    result[valid_mask] = (1 - fractp[valid_mask]) * precip[valid_mask]
-    return result
-
-
 def fractp_op(
         Kc, eto, precip, root, soil, pawc, veg,
         nodata_dict, seasonality_constant):
@@ -1062,19 +1018,19 @@ def fractp_op(
     # and retain their original nodata values.
     # out_nodata is defined above and should never be None.
     valid_mask = (
-        ~utils.array_equals_nodata(Kc, nodata_dict['out_nodata']) &
-        ~utils.array_equals_nodata(root, nodata_dict['out_nodata']) &
-        ~utils.array_equals_nodata(veg, nodata_dict['out_nodata']) &
-        ~utils.array_equals_nodata(precip, 0))
+        ~pygeoprocessing.array_equals_nodata(Kc, nodata_dict['out_nodata']) &
+        ~pygeoprocessing.array_equals_nodata(root, nodata_dict['out_nodata']) &
+        ~pygeoprocessing.array_equals_nodata(veg, nodata_dict['out_nodata']) &
+        ~pygeoprocessing.array_equals_nodata(precip, 0))
     if nodata_dict['eto'] is not None:
-        valid_mask &= ~utils.array_equals_nodata(eto, nodata_dict['eto'])
+        valid_mask &= ~pygeoprocessing.array_equals_nodata(eto, nodata_dict['eto'])
     if nodata_dict['precip'] is not None:
-        valid_mask &= ~utils.array_equals_nodata(precip, nodata_dict['precip'])
+        valid_mask &= ~pygeoprocessing.array_equals_nodata(precip, nodata_dict['precip'])
     if nodata_dict['depth_root'] is not None:
-        valid_mask &= ~utils.array_equals_nodata(
+        valid_mask &= ~pygeoprocessing.array_equals_nodata(
             soil, nodata_dict['depth_root'])
     if nodata_dict['pawc'] is not None:
-        valid_mask &= ~utils.array_equals_nodata(pawc, nodata_dict['pawc'])
+        valid_mask &= ~pygeoprocessing.array_equals_nodata(pawc, nodata_dict['pawc'])
 
     # Compute Budyko Dryness index
     # Use the original AET equation if the land cover type is vegetation
@@ -1121,30 +1077,6 @@ def fractp_op(
     return fractp
 
 
-def pet_op(eto_pix, Kc_pix, eto_nodata, output_nodata):
-    """Calculate the plant potential evapotranspiration.
-
-    Args:
-        eto_pix (numpy.ndarray): a numpy array of ETo
-        Kc_pix (numpy.ndarray): a numpy array of  Kc coefficient
-        precip_nodata (float): nodata value from the precip raster
-        output_nodata (float): nodata value assigned to output of
-            raster_calculator
-
-    Returns:
-        numpy.ndarray of potential evapotranspiration (mm)
-
-    """
-    result = numpy.empty(eto_pix.shape, dtype=numpy.float32)
-    result[:] = output_nodata
-
-    valid_mask = ~utils.array_equals_nodata(Kc_pix, output_nodata)
-    if eto_nodata is not None:
-        valid_mask &= ~utils.array_equals_nodata(eto_pix, eto_nodata)
-    result[valid_mask] = eto_pix[valid_mask] * Kc_pix[valid_mask]
-    return result
-
-
 def compute_watershed_valuation(watershed_results_vector_path, val_df):
     """Compute net present value and energy for the watersheds.
 
diff --git a/src/natcap/invest/carbon.py b/src/natcap/invest/carbon.py
index 624644c0b4..0ea0025a7e 100644
--- a/src/natcap/invest/carbon.py
+++ b/src/natcap/invest/carbon.py
@@ -296,9 +296,6 @@
 
 # -1.0 since carbon stocks are 0 or greater
 _CARBON_NODATA = -1.0
-# use min float32 which is unlikely value to see in a NPV raster
-_VALUE_NODATA = float(numpy.finfo(numpy.float32).min)
-
 
 def execute(args):
     """Carbon.
@@ -435,8 +432,12 @@ def execute(args):
             "Calculate carbon storage for '%s'", output_key)
 
         sum_rasters_task = graph.add_task(
-            _sum_rasters,
-            args=(storage_path_list, file_registry[output_key]),
+            func=pygeoprocessing.raster_map,
+            kwargs=dict(
+                op=sum_op,
+                rasters=storage_path_list,
+                target_path=file_registry[output_key],
+                target_nodata=_CARBON_NODATA),
             target_path_list=[file_registry[output_key]],
             dependent_task_list=carbon_map_task_lookup[scenario_type],
             task_name='sum_rasters_for_total_c_%s' % output_key)
@@ -450,13 +451,15 @@ def execute(args):
             continue
         output_key = 'delta_cur_' + scenario_type
         LOGGER.info("Calculate sequestration scenario '%s'", output_key)
-        storage_path_list = [
-            file_registry['tot_c_cur'],
-            file_registry['tot_c_' + scenario_type]]
 
         diff_rasters_task = graph.add_task(
-            _diff_rasters,
-            args=(storage_path_list, file_registry[output_key]),
+            func=pygeoprocessing.raster_map,
+            kwargs=dict(
+                op=numpy.subtract,  # delta = scenario C - current C
+                rasters=[file_registry['tot_c_' + scenario_type],
+                         file_registry['tot_c_cur']],
+                target_path=file_registry[output_key],
+                target_nodata=_CARBON_NODATA),
             target_path_list=[file_registry[output_key]],
             dependent_task_list=[
                 sum_rasters_task_lookup['cur'],
@@ -513,6 +516,10 @@ def execute(args):
                 file_registry[tmp_filename_key], os_error)
 
 
+# element-wise sum function to pass to raster_map
+def sum_op(*xs): return numpy.sum(xs, axis=0)
+
+
 def _accumulate_totals(raster_path):
     """Sum all non-nodata pixels in `raster_path` and return result."""
     nodata = pygeoprocessing.get_raster_info(raster_path)['nodata'][0]
@@ -522,7 +529,7 @@ def _accumulate_totals(raster_path):
         # the sum.  Users calculated the sum with ArcGIS zonal statistics,
         # noticed a difference and wrote to us about it on the forum.
         raster_sum += numpy.sum(
-            block[~utils.array_equals_nodata(
+            block[~pygeoprocessing.array_equals_nodata(
                     block, nodata)], dtype=numpy.float64)
     return raster_sum
 
@@ -555,43 +562,6 @@ def _generate_carbon_map(
         gdal.GDT_Float32, _CARBON_NODATA, reclass_error_details)
 
 
-def _sum_rasters(storage_path_list, output_sum_path):
-    """Sum all the rasters in `storage_path_list` to `output_sum_path`."""
-    def _sum_op(*storage_arrays):
-        """Sum all the arrays or nodata a pixel stack if one exists."""
-        valid_mask = reduce(
-            lambda x, y: x & y, [
-                ~utils.array_equals_nodata(_, _CARBON_NODATA)
-                for _ in storage_arrays])
-        result = numpy.empty(storage_arrays[0].shape)
-        result[:] = _CARBON_NODATA
-        result[valid_mask] = numpy.sum([
-            _[valid_mask] for _ in storage_arrays], axis=0)
-        return result
-
-    pygeoprocessing.raster_calculator(
-        [(x, 1) for x in storage_path_list], _sum_op, output_sum_path,
-        gdal.GDT_Float32, _CARBON_NODATA)
-
-
-def _diff_rasters(storage_path_list, output_diff_path):
-    """Subtract rasters in `storage_path_list` to `output_sum_path`."""
-    def _diff_op(base_array, future_array):
-        """Subtract future_array from base_array and ignore nodata."""
-        result = numpy.empty(base_array.shape, dtype=numpy.float32)
-        result[:] = _CARBON_NODATA
-        valid_mask = (
-            ~utils.array_equals_nodata(base_array, _CARBON_NODATA) &
-            ~utils.array_equals_nodata(future_array, _CARBON_NODATA))
-        result[valid_mask] = (
-            future_array[valid_mask] - base_array[valid_mask])
-        return result
-
-    pygeoprocessing.raster_calculator(
-        [(x, 1) for x in storage_path_list], _diff_op, output_diff_path,
-        gdal.GDT_Float32, _CARBON_NODATA)
-
-
 def _calculate_valuation_constant(
         lulc_cur_year, lulc_fut_year, discount_rate, rate_change,
         price_per_metric_ton_of_c):
@@ -639,17 +609,10 @@ def _calculate_npv(delta_carbon_path, valuation_constant, npv_out_path):
     Returns:
         None.
     """
-    def _npv_value_op(carbon_array):
-        """Calculate the NPV given carbon storage or loss values."""
-        result = numpy.empty(carbon_array.shape, dtype=numpy.float32)
-        result[:] = _VALUE_NODATA
-        valid_mask = ~utils.array_equals_nodata(carbon_array,  _CARBON_NODATA)
-        result[valid_mask] = carbon_array[valid_mask] * valuation_constant
-        return result
-
-    pygeoprocessing.raster_calculator(
-        [(delta_carbon_path, 1)], _npv_value_op, npv_out_path,
-        gdal.GDT_Float32, _VALUE_NODATA)
+    pygeoprocessing.raster_map(
+        op=lambda carbon: carbon * valuation_constant,
+        rasters=[delta_carbon_path],
+        target_path=npv_out_path)
 
 
 def _generate_report(raster_file_set, model_args, file_registry):
diff --git a/src/natcap/invest/coastal_blue_carbon/coastal_blue_carbon.py b/src/natcap/invest/coastal_blue_carbon/coastal_blue_carbon.py
index 9ce6ef5536..e610aef976 100644
--- a/src/natcap/invest/coastal_blue_carbon/coastal_blue_carbon.py
+++ b/src/natcap/invest/coastal_blue_carbon/coastal_blue_carbon.py
@@ -1514,28 +1514,13 @@ def _calculate_npv(
                 prices_by_year[year] / (
                     (1 + discount_rate) ** years_since_baseline))
 
-        def _npv(*sequestration_matrices):
-            npv = numpy.empty(sequestration_matrices[0].shape,
-                              dtype=numpy.float32)
-            npv[:] = NODATA_FLOAT32_MIN
-
-            matrix_sum = numpy.zeros(npv.shape, dtype=numpy.float32)
-            valid_pixels = numpy.ones(npv.shape, dtype=bool)
-            for matrix in sequestration_matrices:
-                valid_pixels &= ~utils.array_equals_nodata(matrix, NODATA_FLOAT32_MIN)
-                matrix_sum[valid_pixels] += matrix[valid_pixels]
-
-            npv[valid_pixels] = (
-                matrix_sum[valid_pixels] * valuation_factor)
-            return npv
-
-        raster_path_band_tuples = [
-            (path, 1) for (year, path) in net_sequestration_rasters.items() if
-            year <= target_raster_year]
-
-        pygeoprocessing.raster_calculator(
-            raster_path_band_tuples, _npv, target_raster_path,
-            gdal.GDT_Float32, NODATA_FLOAT32_MIN)
+        pygeoprocessing.raster_map(
+            op=lambda *seq_arrays: numpy.sum(
+                seq_arrays, axis=0) * valuation_factor,
+            rasters=[
+                path for year, path in net_sequestration_rasters.items() if
+                year <= target_raster_year],
+            target_path=target_raster_path)
 
 
 def _calculate_stocks_after_baseline_period(
@@ -1562,32 +1547,10 @@ def _calculate_stocks_after_baseline_period(
         ``None``.
 
     """
-    # Both of these values are assumed to be defined from earlier in the
-    # model's execution.
-    baseline_nodata = pygeoprocessing.get_raster_info(
-        baseline_stock_raster_path)['nodata'][0]
-    accum_nodata = pygeoprocessing.get_raster_info(
-        yearly_accumulation_raster_path)['nodata'][0]
-
-    def _calculate_accumulation_over_years(baseline_matrix, accum_matrix):
-        target_matrix = numpy.empty(baseline_matrix.shape, dtype=numpy.float32)
-        target_matrix[:] = NODATA_FLOAT32_MIN
-
-        valid_pixels = (
-            ~utils.array_equals_nodata(baseline_matrix, baseline_nodata) &
-            ~utils.array_equals_nodata(accum_matrix, accum_nodata))
-
-        target_matrix[valid_pixels] = (
-            baseline_matrix[valid_pixels] + (
-                accum_matrix[valid_pixels] * n_years))
-
-        return target_matrix
-
-    pygeoprocessing.raster_calculator(
-        [(baseline_stock_raster_path, 1),
-         (yearly_accumulation_raster_path, 1)],
-        _calculate_accumulation_over_years, target_raster_path,
-        gdal.GDT_Float32, NODATA_FLOAT32_MIN)
+    pygeoprocessing.raster_map(
+        op=lambda baseline, accum: baseline + (accum * n_years),
+        rasters=[baseline_stock_raster_path, yearly_accumulation_raster_path],
+        target_path=target_raster_path)
 
 
 def _calculate_accumulation_over_time(
@@ -1615,9 +1578,9 @@ def _calculate_accumulation_over_time(
     target_matrix[:] = NODATA_FLOAT32_MIN
 
     valid_pixels = (
-        ~utils.array_equals_nodata(annual_biomass_matrix, NODATA_FLOAT32_MIN) &
-        ~utils.array_equals_nodata(annual_soil_matrix, NODATA_FLOAT32_MIN) &
-        ~utils.array_equals_nodata(annual_litter_matrix, NODATA_FLOAT32_MIN))
+        ~pygeoprocessing.array_equals_nodata(annual_biomass_matrix, NODATA_FLOAT32_MIN) &
+        ~pygeoprocessing.array_equals_nodata(annual_soil_matrix, NODATA_FLOAT32_MIN) &
+        ~pygeoprocessing.array_equals_nodata(annual_litter_matrix, NODATA_FLOAT32_MIN))
 
     target_matrix[valid_pixels] = (
         (annual_biomass_matrix[valid_pixels] +
@@ -1715,14 +1678,14 @@ def _track_disturbance(
             disturbance_magnitude_matrix.shape, dtype=numpy.float32)
         disturbed_carbon_volume[:] = NODATA_FLOAT32_MIN
         disturbed_carbon_volume[
-            ~utils.array_equals_nodata(disturbance_magnitude_matrix,
+            ~pygeoprocessing.array_equals_nodata(disturbance_magnitude_matrix,
                            NODATA_FLOAT32_MIN)] = 0.0
 
         if year_of_disturbance_band:
             known_transition_years_matrix = (
                 year_of_disturbance_band.ReadAsArray(**block_info))
             pixels_previously_disturbed = (
-                ~utils.array_equals_nodata(
+                ~pygeoprocessing.array_equals_nodata(
                     known_transition_years_matrix, NODATA_UINT16_MAX))
             year_last_disturbed[pixels_previously_disturbed] = (
                 known_transition_years_matrix[pixels_previously_disturbed])
@@ -1734,9 +1697,9 @@ def _track_disturbance(
 
         stock_matrix = stock_band.ReadAsArray(**block_info)
         pixels_changed_this_year = (
-            ~utils.array_equals_nodata(disturbance_magnitude_matrix, NODATA_FLOAT32_MIN) &
-            ~utils.array_equals_nodata(disturbance_magnitude_matrix, 0.0) &
-            ~utils.array_equals_nodata(stock_matrix, NODATA_FLOAT32_MIN)
+            ~pygeoprocessing.array_equals_nodata(disturbance_magnitude_matrix, NODATA_FLOAT32_MIN) &
+            ~pygeoprocessing.array_equals_nodata(disturbance_magnitude_matrix, 0.0) &
+            ~pygeoprocessing.array_equals_nodata(stock_matrix, NODATA_FLOAT32_MIN)
         )
 
         disturbed_carbon_volume[pixels_changed_this_year] = (
@@ -1801,7 +1764,7 @@ def _record_sequestration(accumulation_matrix, emissions_matrix):
                                                dtype=bool)
         if accumulation_nodata is not None:
             valid_accumulation_pixels &= (
-                ~utils.array_equals_nodata(
+                ~pygeoprocessing.array_equals_nodata(
                     accumulation_matrix, accumulation_nodata))
         target_matrix[valid_accumulation_pixels] += (
             accumulation_matrix[valid_accumulation_pixels])
@@ -1809,7 +1772,7 @@ def _record_sequestration(accumulation_matrix, emissions_matrix):
         valid_emissions_pixels = ~numpy.isclose(emissions_matrix, 0.0)
         if emissions_nodata is not None:
             valid_emissions_pixels &= (
-                ~utils.array_equals_nodata(emissions_matrix, emissions_nodata))
+                ~pygeoprocessing.array_equals_nodata(emissions_matrix, emissions_nodata))
 
         target_matrix[valid_emissions_pixels] = emissions_matrix[
             valid_emissions_pixels] * -1
@@ -1850,9 +1813,9 @@ def _calculate_emissions(
     zero_half_life = numpy.isclose(carbon_half_life_matrix, 0.0)
 
     valid_pixels = (
-        ~utils.array_equals_nodata(
+        ~pygeoprocessing.array_equals_nodata(
             carbon_disturbed_matrix, NODATA_FLOAT32_MIN) &
-        ~utils.array_equals_nodata(
+        ~pygeoprocessing.array_equals_nodata(
             year_of_last_disturbance_matrix, NODATA_UINT16_MAX) &
         ~zero_half_life)
 
@@ -1937,7 +1900,7 @@ def _sum_n_rasters(
             array = band.ReadAsArray(**block_info)
             valid_pixels = slice(None)
             if nodata is not None:
-                valid_pixels = ~utils.array_equals_nodata(array, nodata)
+                valid_pixels = ~pygeoprocessing.array_equals_nodata(array, nodata)
 
             sum_array[valid_pixels] += array[valid_pixels]
             pixels_touched[valid_pixels] = 1
@@ -2130,11 +2093,11 @@ def _reclassify_accumulation(
         valid_pixels = numpy.ones(landuse_transition_from_matrix.shape,
                                   dtype=bool)
         if from_nodata is not None:
-            valid_pixels &= ~utils.array_equals_nodata(
+            valid_pixels &= ~pygeoprocessing.array_equals_nodata(
                 landuse_transition_from_matrix, from_nodata)
 
         if to_nodata is not None:
-            valid_pixels &= ~utils.array_equals_nodata(
+            valid_pixels &= ~pygeoprocessing.array_equals_nodata(
                 landuse_transition_to_matrix, to_nodata)
 
         output_matrix[valid_pixels] = accumulation_rate_matrix[
@@ -2179,39 +2142,12 @@ def _reclassify_disturbance_magnitude(
         ``None``
 
     """
-    from_nodata = pygeoprocessing.get_raster_info(
-        landuse_transition_from_raster)['nodata'][0]
-    to_nodata = pygeoprocessing.get_raster_info(
-        landuse_transition_to_raster)['nodata'][0]
-
-    def _reclassify_disturbance(
-            landuse_transition_from_matrix, landuse_transition_to_matrix):
-        """Pygeoprocessing op to reclassify disturbances."""
-        output_matrix = numpy.empty(landuse_transition_from_matrix.shape,
-                                    dtype=numpy.float32)
-        output_matrix[:] = NODATA_FLOAT32_MIN
-
-        valid_pixels = numpy.ones(landuse_transition_from_matrix.shape,
-                                  dtype=bool)
-        if from_nodata is not None:
-            valid_pixels &= ~utils.array_equals_nodata(
-                landuse_transition_from_matrix, from_nodata)
-
-        if to_nodata is not None:
-            valid_pixels &= ~utils.array_equals_nodata(
-                landuse_transition_to_matrix, to_nodata)
-
-        disturbance_magnitude = disturbance_magnitude_matrix[
-            landuse_transition_from_matrix[valid_pixels],
-            landuse_transition_to_matrix[valid_pixels]].toarray().flatten()
-
-        output_matrix[valid_pixels] = disturbance_magnitude
-        return output_matrix
-
-    pygeoprocessing.raster_calculator(
-        [(landuse_transition_from_raster, 1),
-            (landuse_transition_to_raster, 1)], _reclassify_disturbance,
-        target_raster_path, gdal.GDT_Float32, NODATA_FLOAT32_MIN)
+    pygeoprocessing.raster_map(
+        op=lambda _from, _to: (
+            disturbance_magnitude_matrix[_from, _to].toarray().flatten()),
+        rasters=[landuse_transition_from_raster, landuse_transition_to_raster],
+        target_path=target_raster_path,
+        target_dtype=numpy.float32)
 
 
 @validation.invest_validator
diff --git a/src/natcap/invest/coastal_blue_carbon/preprocessor.py b/src/natcap/invest/coastal_blue_carbon/preprocessor.py
index 8c395e6fba..9cef76c8f0 100644
--- a/src/natcap/invest/coastal_blue_carbon/preprocessor.py
+++ b/src/natcap/invest/coastal_blue_carbon/preprocessor.py
@@ -293,8 +293,8 @@ def _create_transition_table(landcover_df, lulc_snapshot_list,
             # integer type.  When int matrices, we can compare directly to
             # None.
             valid_pixels = (
-                ~utils.array_equals_nodata(from_array, from_nodata) &
-                ~utils.array_equals_nodata(to_array, to_nodata))
+                ~pygeoprocessing.array_equals_nodata(from_array, from_nodata) &
+                ~pygeoprocessing.array_equals_nodata(to_array, to_nodata))
             transition_pairs = transition_pairs.union(
                 set(zip(from_array[valid_pixels].flatten(),
                         to_array[valid_pixels].flatten())))
diff --git a/src/natcap/invest/coastal_vulnerability.py b/src/natcap/invest/coastal_vulnerability.py
index dcfcfab1c4..ff247987df 100644
--- a/src/natcap/invest/coastal_vulnerability.py
+++ b/src/natcap/invest/coastal_vulnerability.py
@@ -3,11 +3,14 @@
 import math
 import os
 import pickle
+import shutil
+import tempfile
 import time
 
 import numpy
 import pandas
 import pygeoprocessing
+import pygeoprocessing.kernels
 import rtree
 import shapely
 import shapely.errors
@@ -1786,7 +1789,7 @@ def extract_bathymetry_along_ray(
         # if nodata value is undefined, all pixels are valid
         valid_mask = slice(None)
         if bathy_nodata is not None:
-            valid_mask = ~utils.array_equals_nodata(values, bathy_nodata)
+            valid_mask = ~pygeoprocessing.array_equals_nodata(values, bathy_nodata)
         if numpy.any(valid_mask):
             # take mean of valids and move on
             value = numpy.mean(values[valid_mask])
@@ -2130,7 +2133,7 @@ def mean_op(array, mask):
 
     _aggregate_raster_values_in_radius(
         base_shore_point_vector_path, positive_dem_path, dem_averaging_radius,
-        target_relief_pickle_path, mean_op)
+        workspace_dir, target_relief_pickle_path, mean_op)
 
     if missing_values:
         LOGGER.warning(
@@ -2145,7 +2148,7 @@ def zero_negative_values(depth_array, nodata):
     """Convert negative values to zero for relief."""
     result_array = numpy.empty_like(depth_array)
     if nodata is not None:
-        valid_mask = ~utils.array_equals_nodata(depth_array, nodata)
+        valid_mask = ~pygeoprocessing.array_equals_nodata(depth_array, nodata)
         result_array[:] = nodata
         result_array[valid_mask] = 0
     else:
@@ -2270,7 +2273,7 @@ def density_op(array, mask):
 
     _aggregate_raster_values_in_radius(
         base_shore_point_vector_path, clipped_projected_pop_path,
-        search_radius, target_pickle_path, density_op)
+        search_radius, workspace_dir, target_pickle_path, density_op)
 
     if missing_values:
         LOGGER.warning(
@@ -2339,6 +2342,7 @@ def _schedule_habitat_tasks(
                       habitat_row['path'],
                       target_habitat_pickle_path,
                       model_resolution,
+                      working_dir,
                       file_suffix),
                 target_path_list=[target_habitat_pickle_path],
                 task_name=f'searching for {_id}'))
@@ -2349,7 +2353,7 @@ def _schedule_habitat_tasks(
 def search_for_raster_habitat(
         base_shore_point_vector_path, search_radius, habitat_rank,
         habitat_id, habitat_raster_path, target_habitat_pickle_path,
-        model_resolution, file_suffix):
+        model_resolution, workspace_dir, file_suffix):
     """Search for habitat raster within a radius of each shore point.
 
     If habitat is present within the search radius, assign the habitat_rank
@@ -2405,7 +2409,8 @@ def habitat_op(array, mask):
 
     _aggregate_raster_values_in_radius(
         base_shore_point_vector_path, clipped_projected_hab_path,
-        search_radius, target_habitat_pickle_path, habitat_op)
+        search_radius, workspace_dir, target_habitat_pickle_path,
+        habitat_op)
     LOGGER.info(f'Finished searching for {habitat_id}')
 
 
@@ -3106,7 +3111,7 @@ def clip_and_project_vector(
 
 def _aggregate_raster_values_in_radius(
         base_point_vector_path, base_raster_path, sample_distance,
-        target_pickle_path, aggregation_op):
+        work_dir, target_pickle_path, aggregation_op):
     """Aggregate raster values in radius around a point.
 
     Do the radial search by constructing a rectangular kernel mask
@@ -3132,6 +3137,7 @@ def _aggregate_raster_values_in_radius(
             projection matching base_point_vector_path.
         sample_distance (float): radius around each point to search
             for valid pixels.
+        work_dir (string): path to directory for intermediate files.
         target_pickle_path (string): path to pickle file storing dict
             keyed by point id.
         aggregation_op (function): takes a signal array and a mask array
@@ -3150,21 +3156,17 @@ def _aggregate_raster_values_in_radius(
 
     # we can assume square pixels at this point because
     # we already warped input raster and defined square pixels
-    pixel_dist = int(abs(
-        sample_distance / (geotransform[1])))
-
-    # kernel dimensions will be 2 * pixel_dist + 1 so that
-    # point feature is always inside the center pixel of the kernel.
-    # create rectangular kernel and a mask so it looks circular.
-    X, Y = numpy.ogrid[:(1 + pixel_dist), :(1 + pixel_dist)]
-    lr_quadrant = numpy.hypot(X, Y)
-    ll_quadrant = numpy.flip(lr_quadrant[:, 1:], axis=1)
-    bottom_half = numpy.concatenate((ll_quadrant, lr_quadrant), axis=1)
-    top_half = numpy.flip(bottom_half[1:, :], axis=0)
-    kernel_index_distances = numpy.concatenate((top_half, bottom_half), axis=0)
-
-    radial_kernel_mask = numpy.where(
-        kernel_index_distances > pixel_dist, False, True)
+    pixel_dist = int(abs(sample_distance / (geotransform[1])))
+
+    temp_dir = tempfile.mkdtemp(dir=work_dir)
+    kernel_path = os.path.join(temp_dir, 'kernel.tif')
+    pygeoprocessing.kernels.dichotomous_kernel(
+        target_kernel_path=kernel_path,
+        max_distance=pixel_dist,
+        normalize=False)
+    radial_kernel_mask = pygeoprocessing.raster_to_numpy_array(
+    kernel_path).astype(bool)
+    shutil.rmtree(temp_dir, ignore_errors=True)
 
     result = {}
     vector = gdal.OpenEx(
@@ -3233,7 +3235,7 @@ def _aggregate_raster_values_in_radius(
             xoff=pixel_x, yoff=pixel_y, win_xsize=win_xsize,
             win_ysize=win_ysize)
         if nodata is not None:
-            kernel_mask &= ~utils.array_equals_nodata(array, nodata)
+            kernel_mask &= ~pygeoprocessing.array_equals_nodata(array, nodata)
 
         result[shore_id] = aggregation_op(array, kernel_mask)
     with open(target_pickle_path, 'wb') as pickle_file:
diff --git a/src/natcap/invest/crop_production_percentile.py b/src/natcap/invest/crop_production_percentile.py
index c23913ef88..c91fb5ae96 100644
--- a/src/natcap/invest/crop_production_percentile.py
+++ b/src/natcap/invest/crop_production_percentile.py
@@ -770,41 +770,12 @@ def calculate_crop_production(lulc_path, yield_path, crop_lucode,
     Returns:
         None
     """
-
-    lulc_nodata = pygeoprocessing.get_raster_info(lulc_path)['nodata'][0]
-    yield_nodata = pygeoprocessing.get_raster_info(yield_path)['nodata'][0]
-
-    def _crop_production_op(lulc_array, yield_array):
-        """Mask in yields that overlap with `crop_lucode`.
-
-        Args:
-            lulc_array (numpy.ndarray): landcover raster values
-            yield_array (numpy.ndarray): interpolated yield raster values
-
-        Returns:
-            numpy.ndarray with float values of yields for the current crop
-
-        """
-        result = numpy.full(lulc_array.shape, _NODATA_YIELD,
-                            dtype=numpy.float32)
-
-        valid_mask = numpy.full(lulc_array.shape, True)
-        if lulc_nodata is not None:
-            valid_mask &= ~utils.array_equals_nodata(lulc_array, lulc_nodata)
-        if yield_nodata is not None:
-            valid_mask &= ~utils.array_equals_nodata(yield_array, yield_nodata)
-        result[valid_mask] = 0
-
-        lulc_mask = lulc_array == crop_lucode
-        result[valid_mask & lulc_mask] = (
-            yield_array[valid_mask & lulc_mask] * pixel_area_ha)
-        return result
-
-    pygeoprocessing.raster_calculator(
-        [(lulc_path, 1), (yield_path, 1)],
-        _crop_production_op,
-        target_path,
-        gdal.GDT_Float32, _NODATA_YIELD),
+    pygeoprocessing.raster_map(
+        op=lambda lulc, _yield: numpy.where(
+            lulc == crop_lucode, _yield * pixel_area_ha, 0),
+        rasters=[lulc_path, yield_path],
+        target_path=target_path,
+        target_nodata=_NODATA_YIELD)
 
 
 def _zero_observed_yield_op(observed_yield_array, observed_yield_nodata):
@@ -823,7 +794,7 @@ def _zero_observed_yield_op(observed_yield_array, observed_yield_nodata):
     result[:] = 0
     valid_mask = slice(None)
     if observed_yield_nodata is not None:
-        valid_mask = ~utils.array_equals_nodata(
+        valid_mask = ~pygeoprocessing.array_equals_nodata(
             observed_yield_array, observed_yield_nodata)
     result[valid_mask] = observed_yield_array[valid_mask]
     return result
@@ -849,7 +820,7 @@ def _mask_observed_yield_op(
     result = numpy.empty(lulc_array.shape, dtype=numpy.float32)
     if landcover_nodata is not None:
         result[:] = observed_yield_nodata
-        valid_mask = ~utils.array_equals_nodata(lulc_array, landcover_nodata)
+        valid_mask = ~pygeoprocessing.array_equals_nodata(lulc_array, landcover_nodata)
         result[valid_mask] = 0
     else:
         result[:] = 0
@@ -923,7 +894,7 @@ def tabulate_results(
                 # if nodata value undefined, assume all pixels are valid
                 valid_mask = numpy.full(yield_block.shape, True)
                 if observed_yield_nodata is not None:
-                    valid_mask = ~utils.array_equals_nodata(
+                    valid_mask = ~pygeoprocessing.array_equals_nodata(
                         yield_block, observed_yield_nodata)
                 production_pixel_count += numpy.count_nonzero(
                     valid_mask & (yield_block > 0))
@@ -943,7 +914,7 @@ def tabulate_results(
                         (yield_percentile_raster_path, 1)):
                     # _NODATA_YIELD will always have a value (defined above)
                     yield_sum += numpy.sum(
-                        yield_block[~utils.array_equals_nodata(
+                        yield_block[~pygeoprocessing.array_equals_nodata(
                             yield_block, _NODATA_YIELD)])
                 production_lookup[yield_percentile_id] = yield_sum
                 result_table.write(",%f" % yield_sum)
@@ -970,7 +941,7 @@ def tabulate_results(
                 (landcover_raster_path, 1)):
             if landcover_nodata is not None:
                 total_area += numpy.count_nonzero(
-                    ~utils.array_equals_nodata(band_values, landcover_nodata))
+                    ~pygeoprocessing.array_equals_nodata(band_values, landcover_nodata))
             else:
                 total_area += band_values.size
         result_table.write(
diff --git a/src/natcap/invest/crop_production_regression.py b/src/natcap/invest/crop_production_regression.py
index 540ef3ed73..da29ec52f4 100644
--- a/src/natcap/invest/crop_production_regression.py
+++ b/src/natcap/invest/crop_production_regression.py
@@ -716,12 +716,14 @@ def execute(args):
                 crop_name, file_suffix))
 
         calc_min_NKP_task = task_graph.add_task(
-            func=pygeoprocessing.raster_calculator,
-            args=([(nitrogen_yield_raster_path, 1),
-                   (phosphorus_yield_raster_path, 1),
-                   (potassium_yield_raster_path, 1)],
-                  _min_op, crop_production_raster_path,
-                  gdal.GDT_Float32, _NODATA_YIELD),
+            func=pygeoprocessing.raster_map,
+            kwargs=dict(
+                op=_min_op,
+                rasters=[nitrogen_yield_raster_path,
+                         phosphorus_yield_raster_path,
+                         potassium_yield_raster_path],
+                target_path=crop_production_raster_path,
+                target_nodata=_NODATA_YIELD),
             target_path_list=[crop_production_raster_path],
             dependent_task_list=dependent_task_list,
             task_name='calc_min_of_NKP')
@@ -855,9 +857,9 @@ def _x_yield_op(
     result = numpy.empty(b_x.shape, dtype=numpy.float32)
     result[:] = _NODATA_YIELD
     valid_mask = (
-        ~utils.array_equals_nodata(y_max, _NODATA_YIELD) &
-        ~utils.array_equals_nodata(b_x, _NODATA_YIELD) &
-        ~utils.array_equals_nodata(c_x, _NODATA_YIELD) &
+        ~pygeoprocessing.array_equals_nodata(y_max, _NODATA_YIELD) &
+        ~pygeoprocessing.array_equals_nodata(b_x, _NODATA_YIELD) &
+        ~pygeoprocessing.array_equals_nodata(c_x, _NODATA_YIELD) &
         (lulc_array == crop_lucode))
     result[valid_mask] = pixel_area_ha * y_max[valid_mask] * (
         1 - b_x[valid_mask] * numpy.exp(
@@ -866,19 +868,8 @@ def _x_yield_op(
     return result
 
 
-def _min_op(y_n, y_p, y_k):
-    """Calculate the min of the three inputs and multiply by Ymax."""
-    result = numpy.empty(y_n.shape, dtype=numpy.float32)
-    result[:] = _NODATA_YIELD
-    valid_mask = (
-        ~utils.array_equals_nodata(y_n, _NODATA_YIELD) &
-        ~utils.array_equals_nodata(y_k, _NODATA_YIELD) &
-        ~utils.array_equals_nodata(y_p, _NODATA_YIELD))
-    result[valid_mask] = (
-        numpy.min(
-            [y_n[valid_mask], y_k[valid_mask], y_p[valid_mask]],
-            axis=0))
-    return result
+"""equation for raster_map: calculate min of inputs and multiply by Ymax."""
+def _min_op(y_n, y_p, y_k): return numpy.min([y_n, y_k, y_p], axis=0)
 
 
 def _zero_observed_yield_op(observed_yield_array, observed_yield_nodata):
@@ -897,7 +888,7 @@ def _zero_observed_yield_op(observed_yield_array, observed_yield_nodata):
     result[:] = 0
     valid_mask = slice(None)
     if observed_yield_nodata is not None:
-        valid_mask = ~utils.array_equals_nodata(
+        valid_mask = ~pygeoprocessing.array_equals_nodata(
             observed_yield_array, observed_yield_nodata)
     result[valid_mask] = observed_yield_array[valid_mask]
     return result
@@ -923,7 +914,7 @@ def _mask_observed_yield_op(
     result = numpy.empty(lulc_array.shape, dtype=numpy.float32)
     if landcover_nodata is not None:
         result[:] = observed_yield_nodata
-        valid_mask = ~utils.array_equals_nodata(lulc_array, landcover_nodata)
+        valid_mask = ~pygeoprocessing.array_equals_nodata(lulc_array, landcover_nodata)
         result[valid_mask] = 0
     else:
         result[:] = 0
@@ -985,7 +976,7 @@ def tabulate_regression_results(
                 # if nodata value undefined, assume all pixels are valid
                 valid_mask = numpy.full(yield_block.shape, True)
                 if observed_yield_nodata is not None:
-                    valid_mask = ~utils.array_equals_nodata(
+                    valid_mask = ~pygeoprocessing.array_equals_nodata(
                         yield_block, observed_yield_nodata)
                 production_pixel_count += numpy.count_nonzero(
                     valid_mask & (yield_block > 0.0))
@@ -1003,7 +994,7 @@ def tabulate_regression_results(
                     (crop_production_raster_path, 1)):
                 yield_sum += numpy.sum(
                     # _NODATA_YIELD will always have a value (defined above)
-                    yield_block[~utils.array_equals_nodata(
+                    yield_block[~pygeoprocessing.array_equals_nodata(
                         yield_block, _NODATA_YIELD)])
             production_lookup['modeled'] = yield_sum
             result_table.write(",%f" % yield_sum)
@@ -1029,7 +1020,7 @@ def tabulate_regression_results(
                 (landcover_raster_path, 1)):
             if landcover_nodata is not None:
                 total_area += numpy.count_nonzero(
-                    ~utils.array_equals_nodata(band_values, landcover_nodata))
+                    ~pygeoprocessing.array_equals_nodata(band_values, landcover_nodata))
             else:
                 total_area += band_values.size
         result_table.write(
diff --git a/src/natcap/invest/delineateit/delineateit.py b/src/natcap/invest/delineateit/delineateit.py
index dae31c3ec4..9644744007 100644
--- a/src/natcap/invest/delineateit/delineateit.py
+++ b/src/natcap/invest/delineateit/delineateit.py
@@ -357,7 +357,7 @@ def _threshold_streams(flow_accum, src_nodata, out_nodata, threshold):
 
     valid_pixels = slice(None)
     if src_nodata is not None:
-        valid_pixels = ~utils.array_equals_nodata(flow_accum, src_nodata)
+        valid_pixels = ~pygeoprocessing.array_equals_nodata(flow_accum, src_nodata)
 
     over_threshold = flow_accum > threshold
     out_matrix[valid_pixels & over_threshold] = 1
diff --git a/src/natcap/invest/forest_carbon_edge_effect.py b/src/natcap/invest/forest_carbon_edge_effect.py
index 34b07964cd..ee73f3ac45 100644
--- a/src/natcap/invest/forest_carbon_edge_effect.py
+++ b/src/natcap/invest/forest_carbon_edge_effect.py
@@ -523,7 +523,7 @@ def combine_carbon_maps(*carbon_maps):
     nodata_mask = numpy.empty(carbon_maps[0].shape, dtype=bool)
     nodata_mask[:] = True
     for carbon_map in carbon_maps:
-        valid_mask = ~utils.array_equals_nodata(carbon_map, NODATA_VALUE)
+        valid_mask = ~pygeoprocessing.array_equals_nodata(carbon_map, NODATA_VALUE)
         nodata_mask &= ~valid_mask
         result[valid_mask] += carbon_map[valid_mask]
     result[nodata_mask] = NODATA_VALUE
@@ -696,26 +696,12 @@ def _map_distance_from_tropical_forest_edge(
     lulc_nodata = pygeoprocessing.get_raster_info(
         base_lulc_raster_path)['nodata']
 
-    forest_mask_nodata = 255
-
-    def mask_non_forest_op(lulc_array):
-        """Convert forest lulc codes to 0.
-        Args:
-            lulc_array (numpy.ndarray): array representing a LULC raster where
-                each forest LULC code is in `forest_codes`.
-        Returns:
-            numpy.ndarray with the same shape as lulc_array. All pixels are
-                0 (forest), 1 (non-forest), or 255 (nodata).
-        """
-        non_forest_mask = ~numpy.isin(lulc_array, forest_codes)
-        nodata_mask = lulc_array == lulc_nodata
-        # where LULC has nodata, set value to nodata value (255)
-        # where LULC has data, set to 0 if LULC is a forest type, 1 if it's not
-        return numpy.where(nodata_mask, forest_mask_nodata, non_forest_mask)
-
-    pygeoprocessing.raster_calculator(
-        [(base_lulc_raster_path, 1)], mask_non_forest_op,
-        target_non_forest_mask_path, gdal.GDT_Byte, forest_mask_nodata)
+    pygeoprocessing.raster_map(
+        op=lambda lulc_array: ~numpy.isin(lulc_array, forest_codes),
+        rasters=[base_lulc_raster_path],
+        target_path=target_non_forest_mask_path,
+        target_dtype=numpy.uint8,
+        target_nodata=255)
 
     # Do the distance transform on non-forest pixels
     # This is the distance from each pixel to the nearest pixel with value 1.
@@ -735,7 +721,7 @@ def mask_non_forest_op(lulc_array):
         # where LULC has nodata, overwrite edge distance with nodata value
         lulc_block = lulc_band.ReadAsArray(**offset_dict)
         distance_block = edge_distance_band.ReadAsArray(**offset_dict)
-        nodata_mask = utils.array_equals_nodata(lulc_block, lulc_nodata)
+        nodata_mask = pygeoprocessing.array_equals_nodata(lulc_block, lulc_nodata)
         distance_block[nodata_mask] = lulc_nodata
         edge_distance_band.WriteArray(
             distance_block,
diff --git a/src/natcap/invest/habitat_quality.py b/src/natcap/invest/habitat_quality.py
index 7c0b15e54b..3f8de3f7e6 100644
--- a/src/natcap/invest/habitat_quality.py
+++ b/src/natcap/invest/habitat_quality.py
@@ -298,17 +298,6 @@
                 "filtered_[THREAT]_aligned.tif": {
                     "about": "Filtered threat raster",
                     "bands": {1: {"type": "ratio"}},
-                },
-                "kernels": {
-                    "type": "directory",
-                    "contents": {
-                        "kernel_[HABITAT]_[SCENARIO].tif": {
-                            "about": (
-                                "Convolution kernel for the given habitat and "
-                                "scenario"),
-                            "bands": {1: {"type": "integer"}}
-                        }
-                    }
                 }
             }
         },
@@ -781,30 +770,12 @@ def _calculate_habitat_quality(deg_hab_raster_list, quality_out_path, ksq):
     Returns:
         None
     """
-    def quality_op(degradation, habitat):
-        """Computes habitat quality given degradation and habitat values."""
-        out_array = numpy.empty_like(degradation)
-        out_array[:] = _OUT_NODATA
-        # Both these rasters are Float32, so the actual pixel values written
-        # might be *slightly* off of _OUT_NODATA but should still be
-        # interpreted as nodata.
-        # _OUT_NODATA (defined above) should never be None, so this is okay
-        valid_pixels = ~(
-            utils.array_equals_nodata(degradation, _OUT_NODATA) |
-            utils.array_equals_nodata(habitat, _OUT_NODATA))
-
-        out_array[valid_pixels] = (
-            habitat[valid_pixels] *
-            (1.0 - (degradation[valid_pixels]**_SCALING_PARAM) /
-                (degradation[valid_pixels]**_SCALING_PARAM + ksq)))
-        return out_array
-
-    deg_hab_raster_band_list = [
-        (path, 1) for path in deg_hab_raster_list]
-
-    pygeoprocessing.raster_calculator(
-        deg_hab_raster_band_list, quality_op, quality_out_path,
-        gdal.GDT_Float32, _OUT_NODATA)
+    pygeoprocessing.raster_map(
+        op=lambda degradation, habitat: (
+            habitat * (1 - (degradation**_SCALING_PARAM) /
+            (degradation**_SCALING_PARAM + ksq))),
+        rasters=deg_hab_raster_list,
+        target_path=quality_out_path)
 
 
 def _calculate_total_degradation(
@@ -822,7 +793,7 @@ def _calculate_total_degradation(
     Returns:
         None
     """
-    def total_degradation(*raster):
+    def total_degradation(*arrays):
         """Computes the total degradation value.
 
         Args:
@@ -842,27 +813,19 @@ def total_degradation(*raster):
         # we can not be certain how many threats the user will enter,
         # so we handle each filtered threat and sensitivity raster
         # in pairs
-        sum_degradation = numpy.zeros(raster[0].shape)
-        for index in range(len(raster) // 2):
+        sum_degradation = numpy.zeros(arrays[0].shape)
+        for index in range(len(arrays) // 2):
             step = index * 2
             sum_degradation += (
-                raster[step] * raster[step + 1] * weight_list[index])
-
-        nodata_mask = numpy.empty(raster[0].shape, dtype=numpy.int8)
-        nodata_mask[:] = 0
-        for array in raster:
-            nodata_mask = nodata_mask | utils.array_equals_nodata(
-                array, _OUT_NODATA)
-
-        # the last element in raster is access
-        return numpy.where(
-            nodata_mask, _OUT_NODATA, sum_degradation * raster[-1])
+                arrays[step] * arrays[step + 1] * weight_list[index])
 
-    deg_raster_band_list = [(path, 1) for path in deg_raster_list]
+        # the last element in arrays is access
+        return sum_degradation * arrays[-1]
 
-    pygeoprocessing.raster_calculator(
-        deg_raster_band_list, total_degradation, deg_sum_raster_path,
-        gdal.GDT_Float32, _OUT_NODATA)
+    pygeoprocessing.raster_map(
+        op=total_degradation,
+        rasters=deg_raster_list,
+        target_path=deg_sum_raster_path)
 
 
 def _compute_rarity_operation(
@@ -904,27 +867,11 @@ def _compute_rarity_operation(
     lulc_area = float(abs(lulc_pixel_size[0]) * abs(lulc_pixel_size[1]))
     lulc_nodata = lulc_raster_info['nodata'][0]
 
-    def trim_op(base, cover_x):
-        """Trim cover_x to the mask of base.
-
-        Args:
-            base (numpy.ndarray): base raster from 'lulc_base'
-            cover_x (numpy.ndarray): either future or current land
-                cover raster from 'lulc_path_band' above
-
-        Returns:
-            _OUT_NODATA where either array has nodata, otherwise cover_x.
-        """
-        result_array = numpy.full(cover_x.shape, _OUT_NODATA)
-        valid_mask = (
-            ~utils.array_equals_nodata(base, base_nodata) &
-            ~utils.array_equals_nodata(cover_x, lulc_nodata))
-        result_array[valid_mask] = cover_x[valid_mask]
-        return result_array
-
-    pygeoprocessing.raster_calculator(
-        [base_lulc_path_band, lulc_path_band], trim_op, new_cover_path[0],
-        gdal.GDT_Float32, _OUT_NODATA)
+    # Trim cover_x to the mask of base.
+    pygeoprocessing.raster_map(
+        op=lambda base, cover_x: cover_x,
+        rasters=[base_lulc_path_band[0], lulc_path_band[0]],
+        target_path=new_cover_path[0])
 
     LOGGER.info('Starting rarity computation on'
                 f' {os.path.basename(lulc_path_band[0])} land cover.')
@@ -1001,7 +948,7 @@ def _raster_values_in_bounds(raster_path_band, lower_bound, upper_bound):
     values_valid = True
 
     for _, raster_block in pygeoprocessing.iterblocks(raster_path_band):
-        nodata_mask = ~utils.array_equals_nodata(raster_block, raster_nodata)
+        nodata_mask = ~pygeoprocessing.array_equals_nodata(raster_block, raster_nodata)
         if ((raster_block[nodata_mask] < lower_bound) |
                 (raster_block[nodata_mask] > upper_bound)).any():
             values_valid = False
@@ -1032,7 +979,7 @@ def _decay_distance(dist_raster_path, max_dist, decay_type, target_path):
         dist_raster_path)['pixel_size']
 
     # convert max distance (given in KM) to meters
-    max_dist_m = max_dist * 1000.0
+    max_dist_m = max_dist * 1000
 
     # convert max distance from meters to the number of pixels that
     # represents on the raster
@@ -1041,31 +988,21 @@ def _decay_distance(dist_raster_path, max_dist, decay_type, target_path):
 
     def linear_op(dist):
         """Linear decay operation."""
-        valid_mask = ~utils.array_equals_nodata(dist, _OUT_NODATA)
-        result = numpy.empty(dist.shape, dtype=numpy.float32)
-        result[:] = _OUT_NODATA
-
-        result[valid_mask] = numpy.where(
-            dist[valid_mask] > max_dist_pixel, 0.0,
-            (max_dist_pixel - dist[valid_mask]) / max_dist_pixel)
-        return result
+        return numpy.where(
+            dist > max_dist_pixel, 0,
+            (max_dist_pixel - dist) / max_dist_pixel)
 
     def exp_op(dist):
         """Exponential decay operation."""
-        valid_mask = ~utils.array_equals_nodata(dist, _OUT_NODATA)
-        result = numpy.empty(dist.shape, dtype=numpy.float32)
-        result[:] = _OUT_NODATA
-
         # Some background on where the 2.99 constant comes from:
         # With the constant of 2.99, the impact of the threat is reduced by
         # 95% (to 5%) at the specified max threat distance. So I suspect it's
         # based on the traditional 95% cutoff that is used in statistics. We
         # could tweak this cutoff (e.g., 99% decay at max distance), if we
         # wanted. - Lisa Mandle
-        result[valid_mask] = numpy.where(
-            dist[valid_mask] > max_dist_pixel, 0.0,
-            numpy.exp((-dist[valid_mask] * 2.99) / max_dist_pixel))
-        return result
+        return numpy.where(
+            dist > max_dist_pixel, 0,
+            numpy.exp((-dist * 2.99) / max_dist_pixel))
 
     if decay_type == 'linear':
         decay_op = linear_op
@@ -1077,9 +1014,10 @@ def exp_op(dist):
             f" either 'linear' or 'exponential'. Input was '{decay_type}' for"
             f" output raster path : '{target_path}'")
 
-    pygeoprocessing.raster_calculator(
-        [(dist_raster_path, 1)], decay_op, target_path, gdal.GDT_Float32,
-        _OUT_NODATA)
+    pygeoprocessing.raster_map(
+        op=decay_op,
+        rasters=[dist_raster_path],
+        target_path=target_path)
 
 
 def _validate_threat_path(threat_path, lulc_key):
diff --git a/src/natcap/invest/hra.py b/src/natcap/invest/hra.py
index 7c21534ce2..99a8f035c9 100644
--- a/src/natcap/invest/hra.py
+++ b/src/natcap/invest/hra.py
@@ -1720,7 +1720,7 @@ def _translate_op(source_rating_array):
             source_rating_array.shape, _TARGET_NODATA_FLOAT32,
             dtype=numpy.float32)
         valid_mask = (
-             (~utils.array_equals_nodata(
+             (~pygeoprocessing.array_equals_nodata(
                  source_rating_array, source_nodata)) &
              (source_rating_array >= 0.0))
         target_rating_array[valid_mask] = source_rating_array[valid_mask]
@@ -2047,7 +2047,7 @@ def _no_buffer(stressor_presence_array):
         # distance and also valid in the source edt.
         pixels_within_buffer = (
             source_edt_block < buffer_distance_in_pixels)
-        nodata_pixels = utils.array_equals_nodata(
+        nodata_pixels = pygeoprocessing.array_equals_nodata(
             source_edt_block, edt_nodata)
         valid_pixels = ((~nodata_pixels) & pixels_within_buffer)
 
@@ -2154,7 +2154,7 @@ def _calc_criteria(attributes_list, habitat_mask_raster_path,
                     # Any habitat pixels with a nodata rating (no rating
                     # specified by the user) should be
                     # interpreted as having a rating of 0.
-                    rating[utils.array_equals_nodata(
+                    rating[pygeoprocessing.array_equals_nodata(
                         rating, rating_nodata)] = 0
                 finally:
                     rating_band = None
@@ -2339,7 +2339,7 @@ def _sum_op(*array_list):
         pixels_have_valid_values = numpy.zeros(result.shape, dtype=bool)
         valid_pixel_count = numpy.zeros(result.shape, dtype=numpy.uint16)
         for array, nodata in zip(array_list, nodata_list):
-            non_nodata_pixels = ~utils.array_equals_nodata(array, nodata)
+            non_nodata_pixels = ~pygeoprocessing.array_equals_nodata(array, nodata)
             pixels_have_valid_values |= non_nodata_pixels
             valid_pixel_count += non_nodata_pixels
 
diff --git a/src/natcap/invest/ndr/ndr.py b/src/natcap/invest/ndr/ndr.py
index 37cfe44cae..9d7370631e 100644
--- a/src/natcap/invest/ndr/ndr.py
+++ b/src/natcap/invest/ndr/ndr.py
@@ -704,8 +704,11 @@ def execute(args):
         task_name='calculate slope')
 
     threshold_slope_task = task_graph.add_task(
-        func=_slope_proportion_and_threshold,
-        args=(f_reg['slope_path'], f_reg['thresholded_slope_path']),
+        func=pygeoprocessing.raster_map,
+        kwargs=dict(
+            op=_slope_proportion_and_threshold_op,
+            rasters=[f_reg['slope_path']],
+            target_path=f_reg['thresholded_slope_path']),
         target_path_list=[f_reg['thresholded_slope_path']],
         dependent_task_list=[calculate_slope_task],
         task_name='threshold slope')
@@ -728,9 +731,13 @@ def execute(args):
         task_name='route s')
 
     s_bar_task = task_graph.add_task(
-        func=s_bar_calculate,
-        args=(f_reg['s_accumulation_path'], f_reg['flow_accumulation_path'],
-              f_reg['s_bar_path']),
+        func=pygeoprocessing.raster_map,
+        kwargs=dict(
+            op=numpy.divide,  # s_bar = s_accum / flow_accum
+            rasters=[f_reg['s_accumulation_path'], f_reg['flow_accumulation_path']],
+            target_path=f_reg['s_bar_path'],
+            target_dtype=numpy.float32,
+            target_nodata=_TARGET_NODATA),
         target_path_list=[f_reg['s_bar_path']],
         dependent_task_list=[s_task, flow_accum_task],
         task_name='calculate s bar')
@@ -744,8 +751,12 @@ def execute(args):
         task_name='d up')
 
     s_inv_task = task_graph.add_task(
-        func=invert_raster_values,
-        args=(f_reg['thresholded_slope_path'], f_reg['s_factor_inverse_path']),
+        func=pygeoprocessing.raster_map,
+        kwargs=dict(
+            op=_inverse_op,
+            rasters=[f_reg['thresholded_slope_path']],
+            target_path=f_reg['s_factor_inverse_path'],
+            target_nodata=_TARGET_NODATA),
         target_path_list=[f_reg['s_factor_inverse_path']],
         dependent_task_list=[threshold_slope_task],
         task_name='s inv')
@@ -810,9 +821,12 @@ def execute(args):
             task_name=f'{nutrient} load')
 
         modified_load_task = task_graph.add_task(
-            func=_multiply_rasters,
-            args=([load_path, f_reg['runoff_proxy_index_path']],
-                  _TARGET_NODATA, modified_load_path),
+            func=pygeoprocessing.raster_map,
+            kwargs=dict(
+                op=_mult_op,
+                rasters=[load_path, f_reg['runoff_proxy_index_path']],
+                target_path=modified_load_path,
+                target_nodata=_TARGET_NODATA),
             target_path_list=[modified_load_path],
             dependent_task_list=[load_task, runoff_proxy_index_task],
             task_name=f'modified load {nutrient}')
@@ -872,8 +886,12 @@ def execute(args):
 
         surface_export_path = f_reg[f'{nutrient}_surface_export_path']
         surface_export_task = task_graph.add_task(
-            func=_calculate_export,
-            args=(surface_load_path, ndr_path, surface_export_path),
+            func=pygeoprocessing.raster_map,
+            kwargs=dict(
+                op=numpy.multiply,  # export = load * ndr
+                rasters=[surface_load_path, ndr_path],
+                target_path=surface_export_path,
+                target_nodata=_TARGET_NODATA),
             target_path_list=[surface_export_path],
             dependent_task_list=[
                 load_task, ndr_task, surface_load_task],
@@ -907,9 +925,12 @@ def execute(args):
                 task_name='sub ndr n')
 
             subsurface_export_task = task_graph.add_task(
-                func=_calculate_export,
-                args=(f_reg['sub_load_n_path'], f_reg['sub_ndr_n_path'],
-                      f_reg['n_subsurface_export_path']),
+                func=pygeoprocessing.raster_map,
+                kwargs=dict(
+                    op=numpy.multiply,  # export = load * ndr
+                    rasters=[f_reg['sub_load_n_path'], f_reg['sub_ndr_n_path']],
+                    target_path=f_reg['n_subsurface_export_path'],
+                    target_nodata=_TARGET_NODATA),
                 target_path_list=[f_reg['n_subsurface_export_path']],
                 dependent_task_list=[
                     subsurface_load_task, subsurface_ndr_task],
@@ -918,9 +939,12 @@ def execute(args):
             # only need to calculate total for nitrogen because
             # phosphorus only has surface export
             total_export_task = task_graph.add_task(
-                func=_sum_rasters,
-                args=([surface_export_path, f_reg['n_subsurface_export_path']],
-                      _TARGET_NODATA, f_reg['n_total_export_path']),
+                func=pygeoprocessing.raster_map,
+                kwargs=dict(
+                    op=_sum_op,
+                    rasters=[surface_export_path, f_reg['n_subsurface_export_path']],
+                    target_path=f_reg['n_total_export_path'],
+                    target_nodata=_TARGET_NODATA),
                 target_path_list=[f_reg['n_total_export_path']],
                 dependent_task_list=[
                     surface_export_task, subsurface_export_task],
@@ -1001,6 +1025,23 @@ def execute(args):
     LOGGER.info(r' (_")  (_/(__)_) (__)  (__) ')
 
 
+# raster_map equation: Multiply a series of arrays element-wise
+def _mult_op(*array_list): return numpy.prod(numpy.stack(array_list), axis=0)
+
+# raster_map equation: Sum a list of arrays element-wise
+def _sum_op(*array_list): return numpy.sum(array_list, axis=0)
+
+# raster_map equation: calculate inverse of S factor
+def _inverse_op(base_val): return numpy.where(base_val == 0, 0, 1 / base_val)
+
+# raster_map equation: rescale and threshold slope between 0.005 and 1
+def _slope_proportion_and_threshold_op(slope):
+    slope_fraction = slope / 100
+    slope_fraction[slope_fraction < 0.005] = 0.005
+    slope_fraction[slope_fraction > 1] = 1
+    return slope_fraction
+
+
 def _create_mask_raster(source_raster_path, source_vector_path,
                         target_raster_path):
     """Create a mask raster from a vector.
@@ -1056,7 +1097,7 @@ def _mask_op(mask, raster):
         result = numpy.full(mask.shape, nodata,
                             dtype=source_raster_info['numpy_type'])
         valid_pixels = (
-            ~utils.array_equals_nodata(raster, nodata) &
+            ~pygeoprocessing.array_equals_nodata(raster, nodata) &
             (mask == 1))
         result[valid_pixels] = raster[valid_pixels]
         return result
@@ -1066,37 +1107,6 @@ def _mask_op(mask, raster):
         target_masked_raster_path, target_dtype, nodata)
 
 
-def _slope_proportion_and_threshold(slope_path, target_threshold_slope_path):
-    """Rescale slope to proportion and threshold to between 0.005 and 1.0.
-
-    Args:
-        slope_path (string): a raster with slope values in percent.
-        target_threshold_slope_path (string): generated raster with slope
-            values as a proportion (100% is 1.0) and thresholded to values
-            between 0.005 and 1.0.
-
-    Returns:
-        None.
-
-    """
-    slope_nodata = pygeoprocessing.get_raster_info(slope_path)['nodata'][0]
-
-    def _slope_proportion_and_threshold_op(slope):
-        """Rescale and threshold slope between 0.005 and 1.0."""
-        valid_mask = ~utils.array_equals_nodata(slope, slope_nodata)
-        result = numpy.empty(valid_mask.shape, dtype=numpy.float32)
-        result[:] = slope_nodata
-        slope_fraction = slope[valid_mask] / 100
-        slope_fraction[slope_fraction < 0.005] = 0.005
-        slope_fraction[slope_fraction > 1.0] = 1.0
-        result[valid_mask] = slope_fraction
-        return result
-
-    pygeoprocessing.raster_calculator(
-        [(slope_path, 1)], _slope_proportion_and_threshold_op,
-        target_threshold_slope_path, gdal.GDT_Float32, slope_nodata)
-
-
 def _add_fields_to_shapefile(field_pickle_map, target_vector_path):
     """Add fields and values to an OGR layer open for writing.
 
@@ -1201,42 +1211,21 @@ def _normalize_raster(base_raster_path_band, target_normalized_raster_path):
         None.
 
     """
-    value_sum = 0.0
-    value_count = 0.0
-    base_nodata = pygeoprocessing.get_raster_info(
-        base_raster_path_band[0])['nodata'][base_raster_path_band[1]-1]
-    for _, raster_block in pygeoprocessing.iterblocks(
-            base_raster_path_band):
-        valid_mask = slice(None)
-        if base_nodata is not None:
-            valid_mask = ~utils.array_equals_nodata(raster_block, base_nodata)
-
-        valid_block = raster_block[valid_mask]
-        value_sum += numpy.sum(valid_block)
-        value_count += valid_block.size
+    value_sum, value_count = pygeoprocessing.raster_reduce(
+        function=lambda sum_count, block:  # calculate both in one pass
+            (sum_count[0] + numpy.sum(block), sum_count[1] + block.size),
+        raster_path_band=base_raster_path_band,
+        initializer=(0, 0))
 
     value_mean = value_sum
-    if value_count > 0.0:
+    if value_count > 0:
         value_mean /= value_count
-    target_nodata = float(numpy.finfo(numpy.float32).min)
-
-    def _normalize_raster_op(array):
-        """Divide values by mean."""
-        result = numpy.empty(array.shape, dtype=numpy.float32)
-        result[:] = target_nodata
-
-        valid_mask = slice(None)
-        if base_nodata is not None:
-            valid_mask = ~utils.array_equals_nodata(array, base_nodata)
-        result[valid_mask] = array[valid_mask]
-        if value_mean != 0:
-            result[valid_mask] /= value_mean
-        return result
 
-    pygeoprocessing.raster_calculator(
-        [base_raster_path_band], _normalize_raster_op,
-        target_normalized_raster_path, gdal.GDT_Float32,
-        target_nodata)
+    pygeoprocessing.raster_map(
+        op=lambda array: array if value_mean == 0 else array / value_mean,
+        rasters=[base_raster_path_band[0]],
+        target_path=target_normalized_raster_path,
+        target_dtype=numpy.float32)
 
 
 def _calculate_load(lulc_raster_path, lucode_to_load, target_load_raster):
@@ -1253,99 +1242,28 @@ def _calculate_load(lulc_raster_path, lucode_to_load, target_load_raster):
         None.
 
     """
-    lulc_raster_info = pygeoprocessing.get_raster_info(lulc_raster_path)
-    nodata_landuse = lulc_raster_info['nodata'][0]
-    cell_area_ha = abs(numpy.prod(lulc_raster_info['pixel_size'])) * 0.0001
+    cell_area_ha = abs(numpy.prod(pygeoprocessing.get_raster_info(
+        lulc_raster_path)['pixel_size'])) * 0.0001
 
     def _map_load_op(lucode_array):
         """Convert unit load to total load & handle nodata."""
         result = numpy.empty(lucode_array.shape)
-        result[:] = _TARGET_NODATA
         for lucode in numpy.unique(lucode_array):
-            if lucode != nodata_landuse:
-                try:
-                    result[lucode_array == lucode] = (
-                        lucode_to_load[lucode] * cell_area_ha)
-                except KeyError:
-                    raise KeyError(
-                        'lucode: %d is present in the landuse raster but '
-                        'missing from the biophysical table' % lucode)
-        return result
-
-    pygeoprocessing.raster_calculator(
-        [(lulc_raster_path, 1)], _map_load_op, target_load_raster,
-        gdal.GDT_Float32, _TARGET_NODATA)
-
-
-def _multiply_rasters(raster_path_list, target_nodata, target_result_path):
-    """Multiply the rasters in `raster_path_list`.
-
-    Args:
-        raster_path_list (list): list of single band raster paths.
-        target_nodata (float): desired target nodata value.
-        target_result_path (string): path to float 32 target raster
-            multiplied where all rasters are not nodata.
-
-    Returns:
-        None.
-
-    """
-    def _mult_op(*array_nodata_list):
-        """Multiply non-nodata stacks."""
-        result = numpy.empty(array_nodata_list[0].shape)
-        result[:] = target_nodata
-        valid_mask = numpy.full(result.shape, True)
-        for array, nodata in zip(*[iter(array_nodata_list)]*2):
-            if nodata is not None:
-                valid_mask &= ~utils.array_equals_nodata(array, nodata)
-        result[valid_mask] = array_nodata_list[0][valid_mask]
-        for array in array_nodata_list[2::2]:
-            result[valid_mask] *= array[valid_mask]
-        return result
-
-    # make a list of (raster_path_band, nodata) tuples, then flatten it
-    path_nodata_list = list(itertools.chain(*[
-        ((path, 1),
-         (pygeoprocessing.get_raster_info(path)['nodata'][0], 'raw'))
-        for path in raster_path_list]))
-    pygeoprocessing.raster_calculator(
-        path_nodata_list, _mult_op, target_result_path,
-        gdal.GDT_Float32, target_nodata)
-
-
-def _sum_rasters(raster_path_list, target_nodata, target_result_path):
-    """Sum two or more rasters pixelwise.
-
-    The result has nodata where any input raster has nodata.
-
-    Args:
-        raster_path_list (list): list of raster paths to sum
-        target_nodata (float): desired target nodata value
-        target_result_path (string): path to write out the sum raster
-
-    Returns:
-        None.
-
-    """
-    nodata_list = [pygeoprocessing.get_raster_info(
-        path)['nodata'][0] for path in raster_path_list]
-
-    def _sum_op(*array_list):
-        """Sum arrays where all are valid."""
-        result = numpy.full(array_list[0].shape, target_nodata,
-                            dtype=numpy.float32)
-        valid_mask = numpy.full(result.shape, True)
-        for array, nodata in zip(array_list, nodata_list):
-            if nodata is not None:
-                valid_mask &= ~numpy.isclose(array, nodata)
-        result[valid_mask] = 0
-        for array in array_list:
-            result[valid_mask] += array[valid_mask]
+            try:
+                result[lucode_array == lucode] = (
+                    lucode_to_load[lucode] * cell_area_ha)
+            except KeyError:
+                raise KeyError(
+                    'lucode: %d is present in the landuse raster but '
+                    'missing from the biophysical table' % lucode)
         return result
 
-    pygeoprocessing.raster_calculator(
-        [(path, 1) for path in raster_path_list],
-        _sum_op, target_result_path, gdal.GDT_Float32, target_nodata)
+    pygeoprocessing.raster_map(
+        op=_map_load_op,
+        rasters=[lulc_raster_path],
+        target_path=target_load_raster,
+        target_dtype=numpy.float32,
+        target_nodata=_TARGET_NODATA)
 
 
 def _map_surface_load(
@@ -1366,9 +1284,6 @@ def _map_surface_load(
         None.
 
     """
-    lulc_raster_info = pygeoprocessing.get_raster_info(lulc_raster_path)
-    nodata_landuse = lulc_raster_info['nodata'][0]
-
     if lucode_to_subsurface_proportion is not None:
         keys = sorted(lucode_to_subsurface_proportion.keys())
         subsurface_values = numpy.array(
@@ -1378,23 +1293,16 @@ def _map_surface_load_op(lucode_array, modified_load_array):
         """Convert unit load to total load & handle nodata."""
         # If we don't have subsurface, just return 0.0.
         if lucode_to_subsurface_proportion is None:
-            return numpy.where(
-                ~utils.array_equals_nodata(lucode_array, nodata_landuse),
-                modified_load_array, _TARGET_NODATA)
-
-        valid_mask = ~utils.array_equals_nodata(lucode_array, nodata_landuse)
-        result = numpy.empty(valid_mask.shape, dtype=numpy.float32)
-        result[:] = _TARGET_NODATA
-        index = numpy.digitize(
-            lucode_array[valid_mask].ravel(), keys, right=True)
-        result[valid_mask] = (
-            modified_load_array[valid_mask] * (1 - subsurface_values[index]))
-        return result
+            return modified_load_array
+        index = numpy.digitize(lucode_array.ravel(), keys, right=True)
+        return modified_load_array * (1 - subsurface_values[index])
 
-    pygeoprocessing.raster_calculator(
-        [(lulc_raster_path, 1), (modified_load_path, 1)],
-        _map_surface_load_op, target_surface_load_path, gdal.GDT_Float32,
-        _TARGET_NODATA)
+    pygeoprocessing.raster_map(
+        op=_map_surface_load_op,
+        rasters=[lulc_raster_path, modified_load_path],
+        target_path=target_surface_load_path,
+        target_dtype=numpy.float32,
+        target_nodata=_TARGET_NODATA)
 
 
 def _map_subsurface_load(
@@ -1413,29 +1321,17 @@ def _map_subsurface_load(
         None.
 
     """
-    lulc_raster_info = pygeoprocessing.get_raster_info(lulc_raster_path)
-    nodata_landuse = lulc_raster_info['nodata'][0]
-
     keys = sorted(numpy.array(list(proportion_subsurface_map)))
     subsurface_permeance_values = numpy.array(
         [proportion_subsurface_map[x] for x in keys])
-
-    def _map_subsurface_load_op(lucode_array, modified_load_array):
-        """Convert unit load to total load & handle nodata."""
-        valid_mask = ~utils.array_equals_nodata(lucode_array, nodata_landuse)
-        result = numpy.empty(valid_mask.shape, dtype=numpy.float32)
-        result[:] = _TARGET_NODATA
-        index = numpy.digitize(
-            lucode_array[valid_mask].ravel(), keys, right=True)
-        result[valid_mask] = (
-            modified_load_array[valid_mask] *
-            subsurface_permeance_values[index])
-        return result
-
-    pygeoprocessing.raster_calculator(
-        [(lulc_raster_path, 1), (modified_load_path, 1)],
-        _map_subsurface_load_op, target_sub_load_path, gdal.GDT_Float32,
-        _TARGET_NODATA)
+    pygeoprocessing.raster_map(
+        op=lambda lulc, modified_load: (
+            modified_load * subsurface_permeance_values[
+                numpy.digitize(lulc.ravel(), keys, right=True)]),
+        rasters=[lulc_raster_path, modified_load_path],
+        target_path=target_sub_load_path,
+        target_dtype=numpy.float32,
+        target_nodata=_TARGET_NODATA)
 
 
 def _map_lulc_to_val_mask_stream(
@@ -1456,110 +1352,27 @@ def _map_lulc_to_val_mask_stream(
     """
     lucodes = sorted(lucodes_to_vals.keys())
     values = numpy.array([lucodes_to_vals[x] for x in lucodes])
-
-    nodata_landuse = pygeoprocessing.get_raster_info(
-        lulc_raster_path)['nodata'][0]
-    nodata_stream = pygeoprocessing.get_raster_info(stream_path)['nodata'][0]
-
-    def _map_eff_op(lucode_array, stream_array):
-        """Map efficiency from LULC and handle nodata/streams."""
-        valid_mask = (
-            ~utils.array_equals_nodata(lucode_array, nodata_landuse) &
-            ~utils.array_equals_nodata(stream_array, nodata_stream))
-        result = numpy.empty(valid_mask.shape, dtype=numpy.float32)
-        result[:] = _TARGET_NODATA
-        index = numpy.digitize(
-            lucode_array[valid_mask].ravel(), lucodes, right=True)
-        result[valid_mask] = (
-            values[index] * (1 - stream_array[valid_mask]))
-        return result
-
-    pygeoprocessing.raster_calculator(
-        ((lulc_raster_path, 1), (stream_path, 1)), _map_eff_op,
-        target_eff_path, gdal.GDT_Float32, _TARGET_NODATA)
-
-
-def s_bar_calculate(
-        s_accumulation_path, flow_accumulation_path, target_s_bar_path):
-    """Calculate bar op which is s/flow."""
-    s_nodata = pygeoprocessing.get_raster_info(
-        s_accumulation_path)['nodata'][0]
-    flow_nodata = pygeoprocessing.get_raster_info(
-        flow_accumulation_path)['nodata'][0]
-
-    def _bar_op(s_accumulation, flow_accumulation):
-        """Calculate bar operation of s_accum / flow_accum."""
-        valid_mask = (
-            ~utils.array_equals_nodata(s_accumulation, s_nodata) &
-            ~utils.array_equals_nodata(flow_accumulation, flow_nodata))
-        result = numpy.empty(valid_mask.shape, dtype=numpy.float32)
-        result[:] = _TARGET_NODATA
-        result[valid_mask] = (
-            s_accumulation[valid_mask] / flow_accumulation[valid_mask])
-        return result
-
-    pygeoprocessing.raster_calculator(
-        ((s_accumulation_path, 1), (flow_accumulation_path, 1)), _bar_op,
-        target_s_bar_path, gdal.GDT_Float32, _TARGET_NODATA)
+    pygeoprocessing.raster_map(
+        op=lambda lulc, stream: (
+            values[numpy.digitize(lulc.ravel(), lucodes, right=True)] *
+            (1 - stream)),
+        rasters=[lulc_raster_path, stream_path],
+        target_path=target_eff_path,
+        target_dtype=numpy.float32,
+        target_nodata=_TARGET_NODATA)
 
 
 def d_up_calculation(s_bar_path, flow_accum_path, target_d_up_path):
     """Calculate d_up = s_bar * sqrt(upslope area)."""
-    s_bar_info = pygeoprocessing.get_raster_info(s_bar_path)
-    s_bar_nodata = s_bar_info['nodata'][0]
-    flow_accum_nodata = pygeoprocessing.get_raster_info(
-        flow_accum_path)['nodata'][0]
-    cell_area_m2 = abs(numpy.prod(s_bar_info['pixel_size']))
-
-    def _d_up_op(s_bar, flow_accumulation):
-        """Calculate d_up index."""
-        valid_mask = (
-            ~utils.array_equals_nodata(s_bar, s_bar_nodata) &
-            ~utils.array_equals_nodata(flow_accumulation, flow_accum_nodata))
-        result = numpy.empty(valid_mask.shape, dtype=numpy.float32)
-        result[:] = _TARGET_NODATA
-        result[valid_mask] = (
-            s_bar[valid_mask] * numpy.sqrt(
-                flow_accumulation[valid_mask] * cell_area_m2))
-        return result
-
-    pygeoprocessing.raster_calculator(
-        [(s_bar_path, 1), (flow_accum_path, 1)], _d_up_op,
-        target_d_up_path, gdal.GDT_Float32, _TARGET_NODATA)
-
-
-def invert_raster_values(base_raster_path, target_raster_path):
-    """Invert (1/x) the values in `base`.
-
-    Args:
-        base_raster_path (string): path to floating point raster.
-        target_raster_path (string): path to created output raster whose
-            values are 1/x of base.
-
-    Returns:
-        None.
-
-    """
-    base_nodata = pygeoprocessing.get_raster_info(
-        base_raster_path)['nodata'][0]
-
-    def _inverse_op(base_val):
-        """Calculate inverse of S factor."""
-        result = numpy.empty(base_val.shape, dtype=numpy.float32)
-        result[:] = _TARGET_NODATA
-        valid_mask = slice(None)
-        if base_nodata is not None:
-            valid_mask = ~utils.array_equals_nodata(base_val, base_nodata)
-
-        zero_mask = base_val == 0.0
-        result[valid_mask & ~zero_mask] = (
-            1.0 / base_val[valid_mask & ~zero_mask])
-        result[zero_mask] = 0.0
-        return result
-
-    pygeoprocessing.raster_calculator(
-        ((base_raster_path, 1),), _inverse_op,
-        target_raster_path, gdal.GDT_Float32, _TARGET_NODATA)
+    cell_area_m2 = abs(numpy.prod(pygeoprocessing.get_raster_info(
+        s_bar_path)['pixel_size']))
+    pygeoprocessing.raster_map(
+        op=lambda s_bar, flow_accum: (
+            s_bar * numpy.sqrt(flow_accum * cell_area_m2)),
+        rasters=[s_bar_path, flow_accum_path],
+        target_path=target_d_up_path,
+        target_dtype=numpy.float32,
+        target_nodata=_TARGET_NODATA)
 
 
 def calculate_ic(d_up_path, d_dn_path, target_ic_path):
@@ -1571,8 +1384,8 @@ def calculate_ic(d_up_path, d_dn_path, target_ic_path):
     def _ic_op(d_up, d_dn):
         """Calculate IC0."""
         valid_mask = (
-            ~utils.array_equals_nodata(d_up, d_up_nodata) &
-            ~utils.array_equals_nodata(d_dn, d_dn_nodata) & (d_up != 0) &
+            ~pygeoprocessing.array_equals_nodata(d_up, d_up_nodata) &
+            ~pygeoprocessing.array_equals_nodata(d_dn, d_dn_nodata) & (d_up != 0) &
             (d_dn != 0))
         result = numpy.empty(valid_mask.shape, dtype=numpy.float32)
         result[:] = ic_nodata
@@ -1592,83 +1405,26 @@ def _calculate_ndr(
     ic_min, ic_max, _, _ = ic_factor_band.GetStatistics(0, 1)
     ic_factor_band = None
     ic_factor_raster = None
-    ic_0_param = (ic_min + ic_max) / 2.0
-    effective_retention_nodata = pygeoprocessing.get_raster_info(
-        effective_retention_path)['nodata'][0]
-    ic_nodata = pygeoprocessing.get_raster_info(ic_factor_path)['nodata'][0]
+    ic_0_param = (ic_min + ic_max) / 2
 
-    def _calculate_ndr_op(effective_retention_array, ic_array):
-        """Calculate NDR."""
-        valid_mask = (
-            ~utils.array_equals_nodata(
-                effective_retention_array, effective_retention_nodata) &
-            ~utils.array_equals_nodata(ic_array, ic_nodata))
-        result = numpy.empty(valid_mask.shape, dtype=numpy.float32)
-        result[:] = _TARGET_NODATA
-        result[valid_mask] = (
-            (1.0 - effective_retention_array[valid_mask]) /
-            (1.0 + numpy.exp(
-                (ic_0_param - ic_array[valid_mask]) / k_param)))
-        return result
-
-    pygeoprocessing.raster_calculator(
-        [(effective_retention_path, 1), (ic_factor_path, 1)],
-        _calculate_ndr_op, target_ndr_path, gdal.GDT_Float32, _TARGET_NODATA)
+    pygeoprocessing.raster_map(
+        op=lambda eff, ic: ((1 - eff) /
+                            (1 + numpy.exp((ic_0_param - ic) / k_param))),
+        rasters=[effective_retention_path, ic_factor_path],
+        target_path=target_ndr_path,
+        target_nodata=_TARGET_NODATA)
 
 
 def _calculate_sub_ndr(
         eff_sub, crit_len_sub, dist_to_channel_path, target_sub_ndr_path):
     """Calculate subsurface: subndr = eff_sub(1-e^(-5*l/crit_len)."""
-    dist_to_channel_nodata = pygeoprocessing.get_raster_info(
-        dist_to_channel_path)['nodata'][0]
-
-    def _sub_ndr_op(dist_to_channel_array):
-        """Calculate subsurface NDR."""
-        # nodata value from this intermediate output should always be
-        # defined by pygeoprocessing, not None
-        valid_mask = ~utils.array_equals_nodata(
-            dist_to_channel_array, dist_to_channel_nodata)
-        result = numpy.empty(valid_mask.shape, dtype=numpy.float32)
-        result[:] = _TARGET_NODATA
-        result[valid_mask] = 1.0 - eff_sub * (
-            1-numpy.exp(-5*dist_to_channel_array[valid_mask]/crit_len_sub))
-        return result
-
-    pygeoprocessing.raster_calculator(
-        [(dist_to_channel_path, 1)], _sub_ndr_op, target_sub_ndr_path,
-        gdal.GDT_Float32, _TARGET_NODATA)
-
-
-def _calculate_export(load_path, ndr_path, target_export_path):
-    """Calculate export.
-
-    Args:
-        load_path (str): path to nutrient load raster
-        ndr_path (str): path to corresponding ndr raster
-        target_export_path (str): path to write out export raster.
-
-    Returns:
-        None
-    """
-    load_nodata = pygeoprocessing.get_raster_info(load_path)['nodata'][0]
-    ndr_nodata = pygeoprocessing.get_raster_info(ndr_path)['nodata'][0]
-
-    def _calculate_export_op(load_array, ndr_array):
-        """Multiply load by NDR."""
-        # these intermediate outputs should always have defined nodata
-        # values assigned by pygeoprocessing
-        valid_mask = ~(
-            utils.array_equals_nodata(load_array, load_nodata) |
-            utils.array_equals_nodata(ndr_array, ndr_nodata))
-        result = numpy.full(valid_mask.shape, _TARGET_NODATA,
-                            dtype=numpy.float32)
-        result[valid_mask] = load_array[valid_mask] * ndr_array[valid_mask]
-        return result
-
-    pygeoprocessing.raster_calculator(
-        [(load_path, 1), (ndr_path, 1)],
-        _calculate_export_op, target_export_path, gdal.GDT_Float32,
-        _TARGET_NODATA)
+    pygeoprocessing.raster_map(
+        op=lambda dist_to_channel: (
+            1 - eff_sub *
+            (1 - numpy.exp(-5 * dist_to_channel / crit_len_sub))),
+        rasters=[dist_to_channel_path],
+        target_path=target_sub_ndr_path,
+        target_nodata=_TARGET_NODATA)
 
 
 def _aggregate_and_pickle_total(
diff --git a/src/natcap/invest/pollination.py b/src/natcap/invest/pollination.py
index adb5d851eb..c5e5fc3692 100644
--- a/src/natcap/invest/pollination.py
+++ b/src/natcap/invest/pollination.py
@@ -9,6 +9,7 @@
 
 import numpy
 import pygeoprocessing
+import pygeoprocessing.kernels
 import taskgraph
 from osgeo import gdal
 from osgeo import ogr
@@ -607,7 +608,7 @@ def execute(args):
             substrate_path_map = scenario_variables[
                 'farm_nesting_substrate_index_path']
         else:
-            dependent_task_list = landcover_substrate_index_tasks.values()
+            dependent_task_list = list(landcover_substrate_index_tasks.values())
             substrate_path_map = scenario_variables[
                 'nesting_substrate_index_path']
 
@@ -616,14 +617,13 @@ def execute(args):
                 intermediate_output_dir,
                 _HABITAT_NESTING_INDEX_FILE_PATTERN % (species, file_suffix)))
 
-        calculate_habitat_nesting_index_op = _CalculateHabitatNestingIndex(
-            substrate_path_map,
-            scenario_variables['species_substrate_index'][species],
-            scenario_variables['habitat_nesting_index_path'][species])
-
         habitat_nesting_tasks[species] = task_graph.add_task(
             task_name=f'calculate_habitat_nesting_{species}',
-            func=calculate_habitat_nesting_index_op,
+            func=_calculate_habitat_nesting_index,
+            args=(
+                substrate_path_map,
+                scenario_variables['species_substrate_index'][species],
+                scenario_variables['habitat_nesting_index_path'][species]),
             dependent_task_list=dependent_task_list,
             target_path_list=[
                 scenario_variables['habitat_nesting_index_path'][species]])
@@ -698,17 +698,14 @@ def execute(args):
             relative_abundance_path = (
                 scenario_variables['relative_floral_abundance_index_path'][
                     season])
-            mult_by_scalar_op = _MultByScalar(
-                scenario_variables['species_foraging_activity'][
-                    (species, season)])
             foraged_flowers_index_task_map[(species, season)] = (
                 task_graph.add_task(
                     task_name=f'calculate_foraged_flowers_{species}_{season}',
-                    func=pygeoprocessing.raster_calculator,
+                    func=_multiply_by_scalar,
                     args=(
-                        [(relative_abundance_path, 1)],
-                        mult_by_scalar_op, foraged_flowers_index_path,
-                        gdal.GDT_Float32, _INDEX_NODATA),
+                        relative_abundance_path,
+                        scenario_variables['species_foraging_activity'][(species, season)],
+                        foraged_flowers_index_path),
                     dependent_task_list=[
                         relative_floral_abudance_task_map[season]],
                     target_path_list=[foraged_flowers_index_path]))
@@ -737,7 +734,7 @@ def execute(args):
             func=pygeoprocessing.raster_calculator,
             args=(
                 foraged_flowers_path_band_list,
-                _SumRasters(), local_foraging_effectiveness_path,
+                _sum_arrays, local_foraging_effectiveness_path,
                 gdal.GDT_Float32, _INDEX_NODATA),
             target_path_list=[
                 local_foraging_effectiveness_path],
@@ -771,8 +768,11 @@ def execute(args):
         except:
             alpha_kernel_raster_task = task_graph.add_task(
                 task_name=f'decay_kernel_raster_{alpha}',
-                func=utils.exponential_decay_kernel_raster,
-                args=(alpha, kernel_path),
+                func=pygeoprocessing.kernels.exponential_decay_kernel,
+                kwargs=dict(
+                    target_kernel_path=kernel_path,
+                    max_distance=alpha * 5,
+                    expected_distance=alpha),
                 target_path_list=[kernel_path])
             alpha_kernel_map[kernel_path] = alpha_kernel_raster_task
 
@@ -803,17 +803,14 @@ def execute(args):
         pollinator_supply_index_path = os.path.join(
             output_dir, _POLLINATOR_SUPPLY_FILE_PATTERN % (
                 species, file_suffix))
-        ps_index_op = _PollinatorSupplyIndexOp(
-            scenario_variables['species_abundance'][species])
         pollinator_supply_task = task_graph.add_task(
             task_name=f'calculate_pollinator_supply_{species}',
-            func=pygeoprocessing.raster_calculator,
+            func=_calculate_pollinator_supply_index,
             args=(
-                [(scenario_variables['habitat_nesting_index_path'][species],
-                  1),
-                 (floral_resources_index_path, 1)], ps_index_op,
-                pollinator_supply_index_path, gdal.GDT_Float32,
-                _INDEX_NODATA),
+                scenario_variables['habitat_nesting_index_path'][species],
+                floral_resources_index_path,
+                scenario_variables['species_abundance'][species],
+                pollinator_supply_index_path),
             dependent_task_list=[
                 floral_resources_task, habitat_nesting_tasks[species]],
             target_path_list=[pollinator_supply_index_path])
@@ -850,13 +847,14 @@ def execute(args):
             pollinator_abundance_task_map[(species, season)] = (
                 task_graph.add_task(
                     task_name=f'calculate_poll_abudance_{species}',
-                    func=pygeoprocessing.raster_calculator,
-                    args=(
-                        [(foraged_flowers_index_path, 1),
-                         (floral_resources_index_path_map[species], 1),
-                         (convolve_ps_path, 1)],
-                        _PollinatorSupplyOp(), pollinator_abundance_path,
-                        gdal.GDT_Float32, _INDEX_NODATA),
+                    func=pygeoprocessing.raster_map,
+                    kwargs=dict(
+                        op=pollinator_supply_op,
+                        rasters=[
+                            foraged_flowers_index_path,
+                            floral_resources_index_path_map[species],
+                            convolve_ps_path],
+                        target_path=pollinator_abundance_path),
                     dependent_task_list=[
                         foraged_flowers_index_task_map[(species, season)],
                         floral_resources_index_task_map[species],
@@ -881,7 +879,7 @@ def execute(args):
             task_name=f'calculate_poll_abudnce_{species}_{season}',
             func=pygeoprocessing.raster_calculator,
             args=(
-                pollinator_abundance_season_path_band_list, _SumRasters(),
+                pollinator_abundance_season_path_band_list, _sum_arrays,
                 total_pollinator_abundance_index_path, gdal.GDT_Float32,
                 _INDEX_NODATA),
             dependent_task_list=[
@@ -930,12 +928,13 @@ def execute(args):
                 season, file_suffix))
         farm_pollinator_season_task_list.append(task_graph.add_task(
             task_name=f'farm_pollinator_{season}',
-            func=pygeoprocessing.raster_calculator,
-            args=(
-                [(half_saturation_raster_path, 1),
-                 (total_pollinator_abundance_index_path, 1)],
-                _OnFarmPollinatorAbundance(), farm_pollinator_season_path,
-                gdal.GDT_Float32, _INDEX_NODATA),
+            func=pygeoprocessing.raster_map,
+            kwargs=dict(
+                op=on_farm_pollinator_abundance_op,
+                rasters=[
+                    half_saturation_raster_path,
+                    total_pollinator_abundance_index_path],
+                target_path=farm_pollinator_season_path),
             dependent_task_list=[
                 half_saturation_task, total_pollinator_abundance_task[season]],
             target_path_list=[farm_pollinator_season_path]))
@@ -949,7 +948,7 @@ def execute(args):
         func=pygeoprocessing.raster_calculator,
         args=(
             [(path, 1) for path in farm_pollinator_season_path_list],
-            _SumRasters(), farm_pollinator_path, gdal.GDT_Float32,
+            _sum_arrays, farm_pollinator_path, gdal.GDT_Float32,
             _INDEX_NODATA),
         dependent_task_list=farm_pollinator_season_task_list,
         target_path_list=[farm_pollinator_path])
@@ -972,11 +971,11 @@ def execute(args):
         output_dir, _TOTAL_POLLINATOR_YIELD_FILE_PATTERN % file_suffix)
     pyt_task = task_graph.add_task(
         task_name='calculate_total_pollinators',
-        func=pygeoprocessing.raster_calculator,
-        args=(
-            [(managed_pollinator_path, 1), (farm_pollinator_path, 1)],
-            _PYTOp(), total_pollinator_yield_path, gdal.GDT_Float32,
-            _INDEX_NODATA),
+        func=pygeoprocessing.raster_map,
+        kwargs=dict(
+            op=pyt_op,
+            rasters=[managed_pollinator_path, farm_pollinator_path],
+            target_path=total_pollinator_yield_path),
         dependent_task_list=[farm_pollinator_task, managed_pollinator_task],
         target_path_list=[total_pollinator_yield_path])
 
@@ -985,11 +984,11 @@ def execute(args):
         output_dir, _WILD_POLLINATOR_YIELD_FILE_PATTERN % file_suffix)
     wild_pollinator_task = task_graph.add_task(
         task_name='calcualte_wild_pollinators',
-        func=pygeoprocessing.raster_calculator,
-        args=(
-            [(managed_pollinator_path, 1), (total_pollinator_yield_path, 1)],
-            _PYWOp(), wild_pollinator_yield_path, gdal.GDT_Float32,
-            _INDEX_NODATA),
+        func=pygeoprocessing.raster_map,
+        kwargs=dict(
+            op=pyw_op,
+            rasters=[managed_pollinator_path, total_pollinator_yield_path],
+            target_path=wild_pollinator_yield_path),
         dependent_task_list=[pyt_task, managed_pollinator_task],
         target_path_list=[wild_pollinator_yield_path])
 
@@ -1068,6 +1067,25 @@ def execute(args):
     task_graph.join()
 
 
+def pollinator_supply_op(
+        foraged_flowers_array, floral_resources_array,
+        convolve_ps_array):
+    """raster_map equation: calculate (RA*fa)/FR * convolve(PS)."""
+    return numpy.where(
+        floral_resources_array == 0, 0,
+        foraged_flowers_array / floral_resources_array * convolve_ps_array)
+
+def on_farm_pollinator_abundance_op(h, pat):
+    """raster_map equation: return (pat * (1 - h)) / (h * (1 - 2*pat)+pat))"""
+    return (pat * (1 - h)) / (h * (1 - 2 * pat) + pat)
+
+# raster_map equation: return min(mp_array+FP_array, 1)
+def pyt_op(mp_array, FP_array): return numpy.minimum(mp_array + FP_array, 1)
+
+# raster_map equation: return max(0,PYT_array-mp_array)
+def pyw_op(mp_array, PYT_array): return numpy.maximum(PYT_array - mp_array, 0)
+
+
 def _rasterize_vector_onto_base(
         base_raster_path, base_vector_path, attribute_id,
         target_raster_path, filter_string=None):
@@ -1367,288 +1385,73 @@ def _parse_scenario_variables(args):
     return result
 
 
-class _CalculateHabitatNestingIndex(object):
-    """Closure for HN(x, s) = max_n(N(x, n) ns(s,n)) calculation."""
-
-    def __init__(
-            self, substrate_path_map, species_substrate_index_map,
-            target_habitat_nesting_index_path):
-        """Define parameters necessary for HN(x,s) calculation.
-
-        Args:
-            substrate_path_map (dict): map substrate name to substrate index
-                raster path. (N(x, n))
-            species_substrate_index_map (dict): map substrate name to
-                scalar value of species substrate suitability. (ns(s,n))
-            target_habitat_nesting_index_path (string): path to target
-                raster
-        """
-        # try to get the source code of __call__ so task graph will recompute
-        # if the function has changed
-        try:
-            self.__name__ = hashlib.sha1(inspect.getsource(
-                    _CalculateHabitatNestingIndex.__call__
-                ).encode('utf-8')).hexdigest()
-        except IOError:
-            # default to the classname if it doesn't work
-            self.__name__ = _CalculateHabitatNestingIndex.__name__
-        self.__name__ += str([
-            substrate_path_map, species_substrate_index_map,
-            target_habitat_nesting_index_path])
-        self.substrate_path_list = [
-            substrate_path_map[substrate_id]
-            for substrate_id in sorted(substrate_path_map)]
-
-        self.species_substrate_suitability_index_array = numpy.array([
-            species_substrate_index_map[substrate_id]
-            for substrate_id in sorted(substrate_path_map)]).reshape(
-                (len(species_substrate_index_map), 1))
-
-        self.target_habitat_nesting_index_path = (
-            target_habitat_nesting_index_path)
-
-    def __call__(self):
-        """Calculate HN(x, s) = max_n(N(x, n) ns(s,n))."""
-        def max_op(*substrate_index_arrays):
-            """Return the max of index_array[n] * ns[n]."""
-            result = numpy.max(
-                numpy.stack([x.flatten() for x in substrate_index_arrays]) *
-                self.species_substrate_suitability_index_array, axis=0)
-            result = result.reshape(substrate_index_arrays[0].shape)
-            nodata_mask = utils.array_equals_nodata(
-                substrate_index_arrays[0], _INDEX_NODATA)
-            result[nodata_mask] = _INDEX_NODATA
-            return result
-
-        pygeoprocessing.raster_calculator(
-            [(x, 1) for x in self.substrate_path_list], max_op,
-            self.target_habitat_nesting_index_path,
-            gdal.GDT_Float32, _INDEX_NODATA)
-
-
-class _SumRasters(object):
-    """Sum all rasters where nodata is 0 unless the entire stack is nodata."""
-
-    def __init__(self):
-        # try to get the source code of __call__ so task graph will recompute
-        # if the function has changed
-        try:
-            self.__name__ = hashlib.sha1(
-                inspect.getsource(
-                    _SumRasters.__call__
-                ).encode('utf-8')).hexdigest()
-        except IOError:
-            # default to the classname if it doesn't work
-            self.__name__ = (
-                _SumRasters.__name__)
-
-    def __call__(self, *array_list):
-        """Calculate sum of array_list and account for nodata."""
-        valid_mask = numpy.zeros(array_list[0].shape, dtype=bool)
-        result = numpy.empty_like(array_list[0])
-        result[:] = 0
-        for array in array_list:
-            local_valid_mask = ~utils.array_equals_nodata(array, _INDEX_NODATA)
-            result[local_valid_mask] += array[local_valid_mask]
-            valid_mask |= local_valid_mask
-        result[~valid_mask] = _INDEX_NODATA
-        return result
-
-
-class _PollinatorSupplyOp(object):
-    """Calc PA=RA*fa/FR * convolve(PS)."""
-
-    def __init__(self):
-        # try to get the source code of __call__ so task graph will recompute
-        # if the function has changed
-        try:
-            self.__name__ = hashlib.sha1(
-                inspect.getsource(
-                    _PollinatorSupplyOp.__call__
-                ).encode('utf-8')).hexdigest()
-        except IOError:
-            # default to the classname if it doesn't work
-            self.__name__ = (
-                _PollinatorSupplyOp.__name__)
-
-    def __call__(
-            self, foraged_flowers_array, floral_resources_array,
-            convolve_ps_array):
-        """Calculating (RA*fa)/FR * convolve(PS)."""
-        valid_mask = ~utils.array_equals_nodata(
-            foraged_flowers_array, _INDEX_NODATA)
-        result = numpy.empty_like(foraged_flowers_array)
-        result[:] = _INDEX_NODATA
-        zero_mask = floral_resources_array == 0
-        result[zero_mask & valid_mask] = 0
-        result_mask = valid_mask & ~zero_mask
-        result[result_mask] = (
-            foraged_flowers_array[result_mask] /
-            floral_resources_array[result_mask] *
-            convolve_ps_array[result_mask])
-        return result
-
-
-class _PollinatorSupplyIndexOp(object):
-    """Calculate PS(x,s) = FR(x,s) * HN(x,s) * sa(s)."""
-
-    def __init__(self, species_abundance):
-        """Create a closure for species abundance to multiply later.
-
-        Args:
-            species_abundance (float): value to use in `__call__` when
-                calculating pollinator abundance.
-
-        Returns:
-            None.
-        """
-        self.species_abundance = species_abundance
-        # try to get the source code of __call__ so task graph will recompute
-        # if the function has changed
-        try:
-            self.__name__ = hashlib.sha1(
-                inspect.getsource(
-                    _PollinatorSupplyIndexOp.__call__
-                ).encode('utf-8')).hexdigest()
-        except IOError:
-            # default to the classname if it doesn't work
-            self.__name__ = (
-                _PollinatorSupplyIndexOp.__name__)
-        self.__name__ += str(species_abundance)
-
-    def __call__(
-            self, floral_resources_array, habitat_nesting_suitability_array):
-        """Calculate f_r * h_n * self.species_abundance."""
-        result = numpy.empty_like(floral_resources_array)
-        result[:] = _INDEX_NODATA
-        valid_mask = ~utils.array_equals_nodata(
-            floral_resources_array, _INDEX_NODATA)
-        result[valid_mask] = (
-            self.species_abundance * floral_resources_array[valid_mask] *
-            habitat_nesting_suitability_array[valid_mask])
-        return result
-
-
-class _MultByScalar(object):
-    """Calculate a raster * scalar.  Mask through nodata."""
-
-    def __init__(self, scalar):
-        """Create a closure for multiplying an array by a scalar.
-
-        Args:
-            scalar (float): value to use in `__call__` when multiplying by
-                its parameter.
-
-        Returns:
-            None.
-        """
-        self.scalar = scalar
-        # try to get the source code of __call__ so task graph will recompute
-        # if the function has changed
-        try:
-            self.__name__ = hashlib.sha1(
-                inspect.getsource(
-                    _MultByScalar.__call__
-                ).encode('utf-8')).hexdigest()
-        except IOError:
-            # default to the classname if it doesn't work
-            self.__name__ = (
-                _MultByScalar.__name__)
-        self.__name__ += str(scalar)
-
-    def __call__(self, array):
-        """Return array * self.scalar accounting for nodata."""
-        result = numpy.empty_like(array)
-        result[:] = _INDEX_NODATA
-        valid_mask = ~utils.array_equals_nodata(array, _INDEX_NODATA)
-        result[valid_mask] = array[valid_mask] * self.scalar
-        return result
-
-
-class _OnFarmPollinatorAbundance(object):
-    """Calculate FP(x) = (PAT * (1 - h)) / (h * (1 - 2*pat)+pat))."""
-
-    def __init__(self):
-        # try to get the source code of __call__ so task graph will recompute
-        # if the function has changed
-        try:
-            self.__name__ = hashlib.sha1(
-                inspect.getsource(
-                    _OnFarmPollinatorAbundance.__call__
-                ).encode('utf-8')).hexdigest()
-        except IOError:
-            # default to the classname if it doesn't work
-            self.__name__ = (
-                _OnFarmPollinatorAbundance.__name__)
-
-    def __call__(self, h_array, pat_array):
-        """Return (pad * (1 - h)) / (h * (1 - 2*pat)+pat)) tolerate nodata."""
-        result = numpy.empty_like(h_array)
-        result[:] = _INDEX_NODATA
-
-        valid_mask = (
-            ~utils.array_equals_nodata(h_array, _INDEX_NODATA) &
-            ~utils.array_equals_nodata(pat_array, _INDEX_NODATA))
-
-        result[valid_mask] = (
-            (pat_array[valid_mask]*(1-h_array[valid_mask])) /
-            (h_array[valid_mask]*(1-2*pat_array[valid_mask]) +
-             pat_array[valid_mask]))
-        return result
-
-
-class _PYTOp(object):
-    """Calculate PYT=min((mp+FP), 1)."""
-
-    def __init__(self):
-        # try to get the source code of __call__ so task graph will recompute
-        # if the function has changed
-        try:
-            self.__name__ = hashlib.sha1(
-                inspect.getsource(
-                    _PYTOp.__call__
-                ).encode('utf-8')).hexdigest()
-        except IOError:
-            # default to the classname if it doesn't work
-            self.__name__ = (
-                _PYTOp.__name__)
-
-    def __call__(self, mp_array, FP_array):
-        """Return min(mp_array+FP_array, 1) accounting for nodata."""
-        valid_mask = ~utils.array_equals_nodata(mp_array, _INDEX_NODATA)
-        result = numpy.empty_like(mp_array)
-        result[:] = _INDEX_NODATA
-        result[valid_mask] = mp_array[valid_mask]+FP_array[valid_mask]
-        min_mask = valid_mask & (result > 1)
-        result[min_mask] = 1
-        return result
-
-
-class _PYWOp(object):
-    """Calculate PYW=max(0,PYT-mp)."""
-
-    def __init__(self):
-        # try to get the source code of __call__ so task graph will recompute
-        # if the function has changed
-        try:
-            self.__name__ = hashlib.sha1(
-                inspect.getsource(
-                    _PYWOp.__call__
-                ).encode('utf-8')).hexdigest()
-        except IOError:
-            # default to the classname if it doesn't work
-            self.__name__ = (
-                _PYWOp.__name__)
-
-    def __call__(self, mp_array, PYT_array):
-        """Return max(0,PYT_array-mp_array) accounting for nodata."""
-        valid_mask = ~utils.array_equals_nodata(mp_array, _INDEX_NODATA)
-        result = numpy.empty_like(mp_array)
-        result[:] = _INDEX_NODATA
-        result[valid_mask] = PYT_array[valid_mask]-mp_array[valid_mask]
-        max_mask = valid_mask & (result < 0)
-        result[max_mask] = 0
-        return result
+def _sum_arrays(*array_list):
+    """Calculate sum of array_list and account for nodata."""
+    valid_mask = numpy.zeros(array_list[0].shape, dtype=bool)
+    result = numpy.empty_like(array_list[0])
+    result[:] = 0
+    for array in array_list:
+        local_valid_mask = ~pygeoprocessing.array_equals_nodata(array, _INDEX_NODATA)
+        result[local_valid_mask] += array[local_valid_mask]
+        valid_mask |= local_valid_mask
+    result[~valid_mask] = _INDEX_NODATA
+    return result
+
+
+def _calculate_habitat_nesting_index(
+        substrate_path_map, species_substrate_index_map,
+        target_habitat_nesting_index_path):
+    """Calculate HN(x, s) = max_n(N(x, n) ns(s,n))."""
+
+    substrate_path_list = [
+        substrate_path_map[substrate_id]
+        for substrate_id in sorted(substrate_path_map)]
+
+    species_substrate_suitability_index_array = numpy.array([
+        species_substrate_index_map[substrate_id]
+        for substrate_id in sorted(substrate_path_map)]).reshape(
+            (len(species_substrate_index_map), 1))
+
+    def max_op(*substrate_index_arrays):
+        """Return the max of index_array[n] * ns[n]."""
+        return numpy.max(
+            numpy.stack([x.flatten() for x in substrate_index_arrays]) *
+            species_substrate_suitability_index_array, axis=0
+        ).reshape(substrate_index_arrays[0].shape)
+
+    pygeoprocessing.raster_map(
+        op=max_op,
+        rasters=substrate_path_list,
+        target_path=target_habitat_nesting_index_path)
+
+
+def _multiply_by_scalar(raster_path, scalar, target_path):
+    """Multiply a raster by a scalar and write out result."""
+    pygeoprocessing.raster_map(
+        op=lambda array: array * scalar,
+        rasters=[raster_path],
+        target_path=target_path)
+
+
+def _calculate_pollinator_supply_index(
+        habitat_nesting_suitability_path, floral_resources_path,
+        species_abundance, target_path):
+    """Calculate pollinator supply index..
+
+    Args:
+        habitat_nesting_suitability_path (str): path to habitat nesting
+            suitability raster
+        floral_resources_path (str): path to floral resources raster
+        species_abundance (float): species abundance value
+        target_path (str): path to write out result raster
+
+    Returns:
+        None.
+    """
+    pygeoprocessing.raster_map(
+        op=lambda f_r, h_n: species_abundance * f_r * h_n,
+        rasters=[habitat_nesting_suitability_path, floral_resources_path],
+        target_path=target_path)
 
 
 @validation.invest_validator
diff --git a/src/natcap/invest/scenario_gen_proximity.py b/src/natcap/invest/scenario_gen_proximity.py
index 62c4acd8b0..bb62ce8fae 100644
--- a/src/natcap/invest/scenario_gen_proximity.py
+++ b/src/natcap/invest/scenario_gen_proximity.py
@@ -11,6 +11,7 @@
 
 import numpy
 import pygeoprocessing
+import pygeoprocessing.kernels
 import scipy
 import taskgraph
 from osgeo import gdal
@@ -445,17 +446,18 @@ 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
-    utils.gaussian_decay_kernel_raster(
-        1.0, tmp_file_registry['gaussian_kernel'])
+    pygeoprocessing.kernels.normal_distribution_kernel(
+        tmp_file_registry['gaussian_kernel'], sigma=1)
 
     # create the output raster first as a copy of the base landcover so it can
     # be looped on for each step
     lulc_raster_info = pygeoprocessing.get_raster_info(base_lulc_path)
     lulc_nodata = lulc_raster_info['nodata'][0]
     mask_nodata = 2
-    pygeoprocessing.raster_calculator(
-        [(base_lulc_path, 1)], lambda x: x, output_landscape_raster_path,
-        gdal.GDT_Int32, lulc_nodata)
+    pygeoprocessing.raster_map(
+        op=lambda x: x,
+        rasters=[base_lulc_path],
+        target_path=output_landscape_raster_path)
 
     # convert everything furthest from edge for each of n_steps
     pixel_area_ha = (
@@ -493,13 +495,14 @@ def _mask_base_op(lulc_array):
                         lulc_array.shape)
                 if invert_mask:
                     base_mask = ~base_mask
-                return numpy.where(
-                    utils.array_equals_nodata(lulc_array, lulc_nodata),
-                    mask_nodata, base_mask)
-            pygeoprocessing.raster_calculator(
-                [(output_landscape_raster_path, 1)], _mask_base_op,
-                tmp_file_registry[mask_id], gdal.GDT_Byte,
-                mask_nodata)
+                return base_mask
+
+            pygeoprocessing.raster_map(
+                op=_mask_base_op,
+                rasters=[output_landscape_raster_path],
+                target_path=tmp_file_registry[mask_id],
+                target_dtype=numpy.uint8,
+                target_nodata=mask_nodata)
 
             # create distance transform for the current mask
             pygeoprocessing.distance_transform_edt(
diff --git a/src/natcap/invest/scenic_quality/scenic_quality.py b/src/natcap/invest/scenic_quality/scenic_quality.py
index 9ff1679c19..3c7bbf799f 100644
--- a/src/natcap/invest/scenic_quality/scenic_quality.py
+++ b/src/natcap/invest/scenic_quality/scenic_quality.py
@@ -617,7 +617,7 @@ def _determine_valid_viewpoints(dem_path, structures_path):
                 (viewpoint[0] - dem_gt[0]) // dem_gt[1]) - block_data['xoff']
             iy_viewpoint = int(
                 (viewpoint[1] - dem_gt[3]) // dem_gt[5]) - block_data['yoff']
-            if utils.array_equals_nodata(
+            if pygeoprocessing.array_equals_nodata(
                     numpy.array(dem_block[iy_viewpoint][ix_viewpoint]),
                     dem_nodata).any():
                 LOGGER.info(
@@ -712,14 +712,14 @@ def _sum_valuation_rasters(dem_path, valuation_filepaths, target_path):
     dem_nodata = pygeoprocessing.get_raster_info(dem_path)['nodata'][0]
 
     def _sum_rasters(dem, *valuation_rasters):
-        valid_dem_pixels = ~utils.array_equals_nodata(dem, dem_nodata)
+        valid_dem_pixels = ~pygeoprocessing.array_equals_nodata(dem, dem_nodata)
         raster_sum = numpy.empty(dem.shape, dtype=numpy.float64)
         raster_sum[:] = _VALUATION_NODATA
         raster_sum[valid_dem_pixels] = 0
 
         for valuation_matrix in valuation_rasters:
             valid_pixels = (
-                ~utils.array_equals_nodata(valuation_matrix, _VALUATION_NODATA)
+                ~pygeoprocessing.array_equals_nodata(valuation_matrix, _VALUATION_NODATA)
                 & valid_dem_pixels)
             raster_sum[valid_pixels] += valuation_matrix[valid_pixels]
         return raster_sum
@@ -848,7 +848,7 @@ def _valuation(distance, visibility):
                                 dtype=numpy.float64) * pixel_size_in_m
 
         valid_distances = (dist_in_m <= max_valuation_radius)
-        nodata = utils.array_equals_nodata(vis_block, vis_nodata)
+        nodata = pygeoprocessing.array_equals_nodata(vis_block, vis_nodata)
         valid_indexes = (valid_distances & (~nodata))
 
         visibility_value[valid_indexes] = _valuation(dist_in_m[valid_indexes],
@@ -920,7 +920,7 @@ def _clip_and_mask_dem(dem_path, aoi_path, target_path, working_dir):
     dem_nodata = dem_raster_info['nodata'][0]
 
     def _mask_op(dem, aoi_mask):
-        valid_pixels = (~utils.array_equals_nodata(dem, dem_nodata) &
+        valid_pixels = (~pygeoprocessing.array_equals_nodata(dem, dem_nodata) &
                         (aoi_mask == 1))
         masked_dem = numpy.empty(dem.shape)
         masked_dem[:] = dem_nodata
@@ -977,7 +977,7 @@ def _count_and_weight_visible_structures(visibility_raster_path_list, weights,
     for block_data in pygeoprocessing.iterblocks((clipped_dem_path, 1),
                                                  offset_only=True):
         dem_block = dem_band.ReadAsArray(**block_data)
-        valid_mask = ~utils.array_equals_nodata(dem_block, dem_nodata)
+        valid_mask = ~pygeoprocessing.array_equals_nodata(dem_block, dem_nodata)
 
         visibility_sum = numpy.empty(dem_block.shape, dtype=numpy.float32)
         visibility_sum[:] = target_nodata
@@ -1055,7 +1055,7 @@ def _mask_zeros(valuation_matrix):
         """Assign zeros to nodata, excluding them from percentile calc."""
         valid_mask = ~numpy.isclose(valuation_matrix, 0.0)
         if raster_nodata is not None:
-            valid_mask &= ~utils.array_equals_nodata(
+            valid_mask &= ~pygeoprocessing.array_equals_nodata(
                 valuation_matrix, raster_nodata)
         visual_quality = numpy.empty(valuation_matrix.shape,
                                      dtype=numpy.float64)
@@ -1081,22 +1081,12 @@ def _mask_zeros(valuation_matrix):
     percentile_bins = numpy.array(percentile_values)
     LOGGER.info('Mapping percentile breaks %s', percentile_bins)
 
-    def _map_percentiles(valuation_matrix):
-        nonzero = (valuation_matrix != 0)
-        nodata = utils.array_equals_nodata(valuation_matrix, raster_nodata)
-        valid_indexes = (~nodata & nonzero)
-        visual_quality = numpy.empty(valuation_matrix.shape,
-                                     dtype=numpy.uint8)
-        visual_quality[:] = _BYTE_NODATA
-        visual_quality[~nonzero & ~nodata] = 0
-        visual_quality[valid_indexes] = numpy.digitize(
-            valuation_matrix[valid_indexes], percentile_bins)
-        return visual_quality
-
-    pygeoprocessing.raster_calculator(
-        [(source_raster_path, 1)], _map_percentiles, target_path,
-        gdal.GDT_Byte, _BYTE_NODATA,
-        raster_driver_creation_tuple=BYTE_GTIFF_CREATION_OPTIONS)
+    pygeoprocessing.raster_map(
+        op=lambda val_array: numpy.where(
+            val_array == 0, 0, numpy.digitize(val_array, percentile_bins)),
+        rasters=[source_raster_path],
+        target_path=target_path,
+        target_dtype=numpy.uint8)
 
 
 @validation.invest_validator
diff --git a/src/natcap/invest/sdr/sdr.py b/src/natcap/invest/sdr/sdr.py
index 6eb954bd5a..5833ac174d 100644
--- a/src/natcap/invest/sdr/sdr.py
+++ b/src/natcap/invest/sdr/sdr.py
@@ -581,8 +581,11 @@ def execute(args):
         task_name='calculate slope')
 
     threshold_slope_task = task_graph.add_task(
-        func=_threshold_slope,
-        args=(f_reg['slope_path'], f_reg['thresholded_slope_path']),
+        func=pygeoprocessing.raster_map,
+        kwargs=dict(
+            op=threshold_slope_op,
+            rasters=[f_reg['slope_path']],
+            target_path=f_reg['thresholded_slope_path']),
         target_path_list=[f_reg['thresholded_slope_path']],
         dependent_task_list=[slope_task],
         task_name='threshold slope')
@@ -631,10 +634,12 @@ def execute(args):
 
     if drainage_present:
         drainage_task = task_graph.add_task(
-            func=_add_drainage(
-                f_reg['stream_path'],
-                f_reg['aligned_drainage_path'],
-                f_reg['stream_and_drainage_path']),
+            func=pygeoprocessing.raster_map,
+            kwargs=dict(
+                op=add_drainage_op,
+                rasters=[f_reg['stream_path'], f_reg['aligned_drainage_path']],
+                target_path=f_reg['stream_and_drainage_path'],
+                target_dtype=numpy.uint8),
             target_path_list=[f_reg['stream_and_drainage_path']],
             dependent_task_list=[stream_task, align_task],
             task_name='add drainage')
@@ -678,11 +683,11 @@ def execute(args):
         task_name='calculate RKLS')
 
     usle_task = task_graph.add_task(
-        func=_calculate_usle,
-        args=(
-            f_reg['rkls_path'],
-            f_reg['cp_factor_path'],
-            f_reg['usle_path']),
+        func=pygeoprocessing.raster_map,
+        kwargs=dict(
+            op=usle_op,
+            rasters=[f_reg['rkls_path'], f_reg['cp_factor_path']],
+            target_path=f_reg['usle_path']),
         target_path_list=[f_reg['usle_path']],
         dependent_task_list=[rkls_task, cp_task],
         task_name='calculate USLE')
@@ -722,10 +727,12 @@ def execute(args):
         task_name='calculate Dup')
 
     inverse_ws_factor_task = task_graph.add_task(
-        func=_calculate_inverse_ws_factor,
-        args=(
-            f_reg['thresholded_slope_path'], f_reg['thresholded_w_path'],
-            f_reg['ws_inverse_path']),
+        func=pygeoprocessing.raster_map,
+        kwargs=dict(
+            op=inverse_ws_op,
+            rasters=[f_reg['thresholded_w_path'],
+                     f_reg['thresholded_slope_path']],
+            target_path=f_reg['ws_inverse_path']),
         target_path_list=[f_reg['ws_inverse_path']],
         dependent_task_list=[threshold_slope_task, threshold_w_task],
         task_name='calculate inverse ws factor')
@@ -762,9 +769,11 @@ def execute(args):
         task_name='calculate sdr')
 
     sed_export_task = task_graph.add_task(
-        func=_calculate_sed_export,
-        args=(
-            f_reg['usle_path'], f_reg['sdr_path'], f_reg['sed_export_path']),
+        func=pygeoprocessing.raster_map,
+        kwargs=dict(
+            op=numpy.multiply,  # export = USLE * SDR
+            rasters=[f_reg['usle_path'], f_reg['sdr_path']],
+            target_path=f_reg['sed_export_path']),
         target_path_list=[f_reg['sed_export_path']],
         dependent_task_list=[usle_task, sdr_task],
         task_name='calculate sed export')
@@ -789,19 +798,23 @@ def execute(args):
         task_name='sediment deposition')
 
     avoided_erosion_task = task_graph.add_task(
-        func=_calculate_avoided_erosion,
-        args=(
-            f_reg['rkls_path'], f_reg['usle_path'],
-            f_reg['avoided_erosion_path']),
+        func=pygeoprocessing.raster_map,
+        kwargs=dict(
+            op=numpy.subtract,  # avoided erosion = rkls - usle
+            rasters=[f_reg['rkls_path'], f_reg['usle_path']],
+            target_path=f_reg['avoided_erosion_path']),
         dependent_task_list=[rkls_task, usle_task],
         target_path_list=[f_reg['avoided_erosion_path']],
         task_name='calculate avoided erosion')
 
     avoided_export_task = task_graph.add_task(
-        func=_calculate_avoided_export,
-        args=(
-            f_reg['avoided_erosion_path'], f_reg['sdr_path'],
-            f_reg['sed_deposition_path'], f_reg['avoided_export_path']),
+        func=pygeoprocessing.raster_map,
+        kwargs=dict(
+            op=_avoided_export_op,
+            rasters=[f_reg['avoided_erosion_path'],
+                     f_reg['sdr_path'],
+                     f_reg['sed_deposition_path']],
+            target_path=f_reg['avoided_export_path']),
         dependent_task_list=[avoided_erosion_task, sdr_task,
                              sed_deposition_task],
         target_path_list=[f_reg['avoided_export_path']],
@@ -832,103 +845,41 @@ def execute(args):
     task_graph.join()
 
 
-def _calculate_avoided_export(
-        avoided_erosion_path, sdr_path, sed_deposition_path,
-        target_avoided_export_path):
-    """Calculate total retention.
+def _avoided_export_op(avoided_erosion, sdr, sed_deposition):
+    """raster_map equation: calculate total retention.
 
     Args:
-        avoided_erosion_path (string): The string path to an avoided erosion
-            raster.
-        sdr_path (string): The path to a raster containing computed SDR values.
-        sed_deposition_path (string): The path to a raster containing computed
-            sediment deposition values.
-        target_avoided_export_path (string): The path to the calculated total
-            retention raster, produced by this function.
+        avoided_erosion (numpy.array): Avoided erosion values.
+        sdr (numpy.array): SDR values.
+        sed_deposition (numpy.array): Sediment deposition values.
 
     Returns:
-        ``None``
+        A ``numpy.array`` of computed total retention matching the shape of
+        the input numpy arrays.
     """
-    avoided_erosion_nodata = pygeoprocessing.get_raster_info(
-        avoided_erosion_path)['nodata'][0]
-    sdr_nodata = pygeoprocessing.get_raster_info(sdr_path)['nodata'][0]
-    sed_deposition_nodata = pygeoprocessing.get_raster_info(
-        sed_deposition_path)['nodata'][0]
-
-    def _avoided_export_function(avoided_erosion, sdr, sed_deposition):
-        """Calculate total retention.
-
-        Args:
-            avoided_erosion (numpy.array): Avoided erosion values.
-            sdr (numpy.array): SDR values.
-            sed_deposition (numpy.array): Sediment deposition values.
-
-        Returns:
-            A ``numpy.array`` of computed total retention matching the shape of
-            the input numpy arrays.
-        """
-        result = numpy.full(avoided_erosion.shape, _TARGET_NODATA,
-                            dtype=numpy.float32)
-        valid_mask = (
-            (~utils.array_equals_nodata(avoided_erosion,
-                                        avoided_erosion_nodata)) &
-            (~utils.array_equals_nodata(sdr, sdr_nodata)) &
-            (~utils.array_equals_nodata(sed_deposition,
-                                        sed_deposition_nodata)))
-
-        # avoided_erosion represents RLKS - RKLSCP (where RKLSCP is also
-        # known as a modified USLE)
-        result[valid_mask] = (
-            avoided_erosion[valid_mask] * sdr[valid_mask] +
-            sed_deposition[valid_mask])
-        return result
-
-    pygeoprocessing.raster_calculator(
-        [(avoided_erosion_path, 1), (sdr_path, 1),
-         (sed_deposition_path, 1)],
-        _avoided_export_function, target_avoided_export_path,
-        gdal.GDT_Float32, _TARGET_NODATA)
-
-
-def _calculate_avoided_erosion(
-        rkls_path, usle_path, target_avoided_erosion_path):
-    """Calculate avoided erosion.
+    # avoided_erosion represents RLKS - RKLSCP (where RKLSCP is also
+    # known as a modified USLE)
+    return avoided_erosion * sdr + sed_deposition
 
+def add_drainage_op(stream, drainage):
+    """raster_map equation: add drainage mask to stream layer.
 
     Args:
-        rkls_path (string): The path to the RKLS raster on disk.
-        usle_path (string): The path to the USLE raster on disk.
-        target_avoided_erosion_path (string): The path to the target
-            local retention raster, created by this function.
+        stream (numpy.array): binary array where 1 indicates
+            a stream, and 0 is a valid landscape pixel but not a stream.
+        drainage (numpy.array): binary array where 1 indicates any water
+            reaching that pixel drains to a stream.
 
     Returns:
-        ``None``
+        numpy.array combination of stream and drainage
     """
-    rkls_nodata = pygeoprocessing.get_raster_info(rkls_path)['nodata'][0]
-    usle_nodata = pygeoprocessing.get_raster_info(usle_path)['nodata'][0]
+    return numpy.where(drainage == 1, 1, stream)
 
-    def _avoided_erosion_function(rkls, usle):
-        """Calculate avoided erosion.
+# raster_map equation: calculate USLE
+def usle_op(rkls, cp_factor): return rkls * cp_factor
 
-        Args:
-            rkls (numpy.array): Computed RKLS values.
-            usle (numpy.array): Computed USLE values.
-
-        Returns:
-            A ``numpy.array`` of the same size and shape of the input numpy
-            arrays containing computed avoided erosion values.
-        """
-        result = numpy.full(rkls.shape, _TARGET_NODATA, dtype=numpy.float32)
-        valid_mask = (
-            (~utils.array_equals_nodata(rkls, rkls_nodata)) &
-            (~utils.array_equals_nodata(usle, usle_nodata)))
-
-        result[valid_mask] = rkls[valid_mask] - usle[valid_mask]
-        return result
-
-    pygeoprocessing.raster_calculator(
-        [(rkls_path, 1), (usle_path, 1)], _avoided_erosion_function,
-        target_avoided_erosion_path, gdal.GDT_Float32, _TARGET_NODATA)
+# raster_map equation: calculate the inverse ws factor
+def inverse_ws_op(w_factor, s_factor): return 1 / (w_factor * s_factor)
 
 
 def _calculate_what_drains_to_stream(
@@ -979,10 +930,10 @@ def _what_drains_to_stream(flow_dir_mfd, dist_to_channel):
         """
         drains_to_stream = numpy.full(
             flow_dir_mfd.shape, _BYTE_NODATA, dtype=numpy.uint8)
-        valid_flow_dir = ~utils.array_equals_nodata(
+        valid_flow_dir = ~pygeoprocessing.array_equals_nodata(
             flow_dir_mfd, flow_dir_mfd_nodata)
         valid_dist_to_channel = (
-            ~utils.array_equals_nodata(
+            ~pygeoprocessing.array_equals_nodata(
                 dist_to_channel, dist_to_channel_nodata) &
             valid_flow_dir)
 
@@ -1060,15 +1011,11 @@ def _calculate_ls_factor(
         None
 
     """
-    slope_nodata = pygeoprocessing.get_raster_info(slope_path)['nodata'][0]
-
-    flow_accumulation_info = pygeoprocessing.get_raster_info(
-        flow_accumulation_path)
-    flow_accumulation_nodata = flow_accumulation_info['nodata'][0]
-    cell_size = abs(flow_accumulation_info['pixel_size'][0])
+    cell_size = abs(pygeoprocessing.get_raster_info(
+        flow_accumulation_path)['pixel_size'][0])
     cell_area = cell_size ** 2
 
-    def ls_factor_function(percent_slope, flow_accumulation, l_max):
+    def ls_factor_function(percent_slope, flow_accumulation):
         """Calculate the LS factor.
 
         Args:
@@ -1080,15 +1027,6 @@ def ls_factor_function(percent_slope, flow_accumulation, l_max):
             ls_factor
 
         """
-        # avg aspect intermediate output should always have a defined
-        # nodata value from pygeoprocessing
-        valid_mask = (
-            ~utils.array_equals_nodata(percent_slope, slope_nodata) &
-            ~utils.array_equals_nodata(
-                flow_accumulation, flow_accumulation_nodata))
-        result = numpy.empty(valid_mask.shape, dtype=numpy.float32)
-        result[:] = _TARGET_NODATA
-
         # Although Desmet & Govers (1996) discusses "upstream contributing
         # area", this is not strictly defined. We decided to use the square
         # root of the upstream contributing area here as an estimate, which
@@ -1097,9 +1035,8 @@ def ls_factor_function(percent_slope, flow_accumulation, l_max):
         # We subtract 1 from the flow accumulation because FA includes itself
         # in its count of pixels upstream and our LS factor equation wants only
         # those pixels that are strictly upstream.
-        contributing_area = numpy.sqrt(
-            (flow_accumulation[valid_mask]-1) * cell_area)
-        slope_in_radians = numpy.arctan(percent_slope[valid_mask] / 100.0)
+        contributing_area = numpy.sqrt((flow_accumulation - 1) * cell_area)
+        slope_in_radians = numpy.arctan(percent_slope / 100)
 
         aspect_length = (numpy.fabs(numpy.sin(slope_in_radians)) +
                          numpy.fabs(numpy.cos(slope_in_radians)))
@@ -1107,7 +1044,7 @@ def ls_factor_function(percent_slope, flow_accumulation, l_max):
         # From Equation 4 in "Extension and validation of a geographic
         # information system ..."
         slope_factor = numpy.where(
-            percent_slope[valid_mask] < 9.0,
+            percent_slope < 9,
             10.8 * numpy.sin(slope_in_radians) + 0.03,
             16.8 * numpy.sin(slope_in_radians) - 0.5)
 
@@ -1121,10 +1058,9 @@ def ls_factor_function(percent_slope, flow_accumulation, l_max):
         slope_table = numpy.array([1., 3.5, 5., 9.])
         m_table = numpy.array([0.2, 0.3, 0.4, 0.5])
         # mask where slopes are larger than lookup table
-        big_slope_mask = percent_slope[valid_mask] > slope_table[-1]
+        big_slope_mask = percent_slope > slope_table[-1]
         m_indexes = numpy.digitize(
-            percent_slope[valid_mask][~big_slope_mask], slope_table,
-            right=True)
+            percent_slope[~big_slope_mask], slope_table, right=True)
         m_exp = numpy.empty(big_slope_mask.shape, dtype=numpy.float32)
         m_exp[big_slope_mask] = (
             beta[big_slope_mask] / (1 + beta[big_slope_mask]))
@@ -1139,14 +1075,12 @@ def ls_factor_function(percent_slope, flow_accumulation, l_max):
         # threshold L factor to l_max
         l_factor[l_factor > l_max] = l_max
 
-        result[valid_mask] = l_factor * slope_factor
-        return result
+        return l_factor * slope_factor
 
-    pygeoprocessing.raster_calculator(
-        [(path, 1) for path in [slope_path, flow_accumulation_path]] + [
-            (l_max, 'raw')],
-        ls_factor_function, target_ls_factor_path, gdal.GDT_Float32,
-        _TARGET_NODATA)
+    pygeoprocessing.raster_map(
+        op=ls_factor_function,
+        rasters=[slope_path, flow_accumulation_path],
+        target_path=target_ls_factor_path)
 
 
 def _calculate_rkls(
@@ -1198,13 +1132,13 @@ def rkls_function(ls_factor, erosivity, erodibility, stream):
         """
         rkls = numpy.empty(ls_factor.shape, dtype=numpy.float32)
         nodata_mask = (
-            ~utils.array_equals_nodata(ls_factor, _TARGET_NODATA) &
-            ~utils.array_equals_nodata(stream, stream_nodata))
+            ~pygeoprocessing.array_equals_nodata(ls_factor, _TARGET_NODATA) &
+            ~pygeoprocessing.array_equals_nodata(stream, stream_nodata))
         if erosivity_nodata is not None:
-            nodata_mask &= ~utils.array_equals_nodata(
+            nodata_mask &= ~pygeoprocessing.array_equals_nodata(
                 erosivity, erosivity_nodata)
         if erodibility_nodata is not None:
-            nodata_mask &= ~utils.array_equals_nodata(
+            nodata_mask &= ~pygeoprocessing.array_equals_nodata(
                 erodibility, erodibility_nodata)
 
         valid_mask = nodata_mask & (stream == 0)
@@ -1225,62 +1159,15 @@ def rkls_function(ls_factor, erosivity, erodibility, stream):
         rkls_function, rkls_path, gdal.GDT_Float32, _TARGET_NODATA)
 
 
-def _threshold_slope(slope_path, out_thresholded_slope_path):
-    """Threshold the slope between 0.005 and 1.0.
-
-    Args:
-        slope_path (string): path to a raster of slope in percent
-        out_thresholded_slope_path (string): path to output raster of
-            thresholded slope between 0.005 and 1.0
-
-    Returns:
-        None
-
-    """
-    slope_nodata = pygeoprocessing.get_raster_info(slope_path)['nodata'][0]
-
-    def threshold_slope(slope):
-        """Convert slope to m/m and clamp at 0.005 and 1.0.
-
-        As desribed in Cavalli et al., 2013.
-        """
-        valid_slope = ~utils.array_equals_nodata(slope, slope_nodata)
-        slope_m = slope[valid_slope] / 100.0
-        slope_m[slope_m < 0.005] = 0.005
-        slope_m[slope_m > 1.0] = 1.0
-        result = numpy.empty(valid_slope.shape, dtype=numpy.float32)
-        result[:] = slope_nodata
-        result[valid_slope] = slope_m
-        return result
-
-    pygeoprocessing.raster_calculator(
-        [(slope_path, 1)], threshold_slope, out_thresholded_slope_path,
-        gdal.GDT_Float32, slope_nodata)
-
-
-def _add_drainage(stream_path, drainage_path, out_stream_and_drainage_path):
-    """Combine stream and drainage masks into one raster mask.
-
-    Args:
-        stream_path (string): path to stream raster mask where 1 indicates
-            a stream, and 0 is a valid landscape pixel but not a stream.
-        drainage_raster_path (string): path to 1/0 mask of drainage areas.
-            1 indicates any water reaching that pixel drains to a stream.
-        out_stream_and_drainage_path (string): output raster of a logical
-            OR of stream and drainage inputs
-
-    Returns:
-        None
+def threshold_slope_op(slope):
+    """raster_map equation: convert slope to m/m and clamp at 0.005 and 1.0.
 
+    As desribed in Cavalli et al., 2013.
     """
-    def add_drainage_op(stream, drainage):
-        """Add drainage mask to stream layer."""
-        return numpy.where(drainage == 1, 1, stream)
-
-    stream_nodata = pygeoprocessing.get_raster_info(stream_path)['nodata'][0]
-    pygeoprocessing.raster_calculator(
-        [(path, 1) for path in [stream_path, drainage_path]], add_drainage_op,
-        out_stream_and_drainage_path, gdal.GDT_Byte, stream_nodata)
+    slope_m = slope / 100
+    slope_m[slope_m < 0.005] = 0.005
+    slope_m[slope_m > 1] = 1
+    return slope_m
 
 
 def _calculate_w(
@@ -1316,17 +1203,10 @@ def _calculate_w(
         (lulc_path, 1), lulc_to_c, w_factor_path, gdal.GDT_Float32,
         _TARGET_NODATA, reclass_error_details)
 
-    def threshold_w(w_val):
-        """Threshold w to 0.001."""
-        w_val_copy = w_val.copy()
-        nodata_mask = utils.array_equals_nodata(w_val, _TARGET_NODATA)
-        w_val_copy[w_val < 0.001] = 0.001
-        w_val_copy[nodata_mask] = _TARGET_NODATA
-        return w_val_copy
-
-    pygeoprocessing.raster_calculator(
-        [(w_factor_path, 1)], threshold_w, out_thresholded_w_factor_path,
-        gdal.GDT_Float32, _TARGET_NODATA)
+    pygeoprocessing.raster_map(
+        op=lambda w_val: numpy.where(w_val < 0.001, 0.001, w_val),
+        rasters=[w_factor_path],
+        target_path=out_thresholded_w_factor_path)
 
 
 def _calculate_cp(lulc_to_cp, lulc_path, cp_factor_path):
@@ -1357,24 +1237,6 @@ def _calculate_cp(lulc_to_cp, lulc_path, cp_factor_path):
         _TARGET_NODATA, reclass_error_details)
 
 
-def _calculate_usle(
-        rkls_path, cp_factor_path, out_usle_path):
-    """Calculate USLE, multiply RKLS by CP."""
-    def usle_op(rkls, cp_factor):
-        """Calculate USLE."""
-        result = numpy.empty(rkls.shape, dtype=numpy.float32)
-        result[:] = _TARGET_NODATA
-        valid_mask = (
-            ~utils.array_equals_nodata(rkls, _TARGET_NODATA) &
-            ~utils.array_equals_nodata(cp_factor, _TARGET_NODATA))
-        result[valid_mask] = rkls[valid_mask] * cp_factor[valid_mask]
-        return result
-
-    pygeoprocessing.raster_calculator(
-        [(rkls_path, 1), (cp_factor_path, 1)], usle_op,
-        out_usle_path, gdal.GDT_Float32, _TARGET_NODATA)
-
-
 def _calculate_bar_factor(
         flow_direction_path, factor_path, flow_accumulation_path,
         accumulation_path, out_bar_path):
@@ -1397,9 +1259,6 @@ def _calculate_bar_factor(
         None.
 
     """
-    flow_accumulation_nodata = pygeoprocessing.get_raster_info(
-        flow_accumulation_path)['nodata'][0]
-
     LOGGER.debug("doing flow accumulation mfd on %s", factor_path)
     # manually setting compression to DEFLATE because we got some LZW
     # errors when testing with large data.
@@ -1410,22 +1269,10 @@ def _calculate_bar_factor(
             'TILED=YES', 'BIGTIFF=YES', 'COMPRESS=DEFLATE',
             'PREDICTOR=3']))
 
-    def bar_op(base_accumulation, flow_accumulation):
-        """Aggregate accumulation from base divided by the flow accum."""
-        result = numpy.empty(base_accumulation.shape, dtype=numpy.float32)
-        # flow accumulation intermediate output should always have a defined
-        # nodata value from pygeoprocessing
-        valid_mask = ~(
-            utils.array_equals_nodata(base_accumulation, _TARGET_NODATA) |
-            utils.array_equals_nodata(
-                flow_accumulation, flow_accumulation_nodata))
-        result[:] = _TARGET_NODATA
-        result[valid_mask] = (
-            base_accumulation[valid_mask] / flow_accumulation[valid_mask])
-        return result
-    pygeoprocessing.raster_calculator(
-        [(accumulation_path, 1), (flow_accumulation_path, 1)], bar_op,
-        out_bar_path, gdal.GDT_Float32, _TARGET_NODATA)
+    pygeoprocessing.raster_map(
+        op=lambda base_accum, flow_accum: base_accum / flow_accum,
+        rasters=[accumulation_path, flow_accumulation_path],
+        target_path=out_bar_path)
 
 
 def _calculate_d_up(
@@ -1433,31 +1280,11 @@ def _calculate_d_up(
     """Calculate w_bar * s_bar * sqrt(flow accumulation * cell area)."""
     cell_area = abs(
         pygeoprocessing.get_raster_info(w_bar_path)['pixel_size'][0])**2
-    flow_accumulation_nodata = pygeoprocessing.get_raster_info(
-        flow_accumulation_path)['nodata'][0]
-
-    def d_up_op(w_bar, s_bar, flow_accumulation):
-        """Calculate the d_up index.
-
-        w_bar * s_bar * sqrt(upslope area)
-
-        """
-        valid_mask = (
-            ~utils.array_equals_nodata(w_bar, _TARGET_NODATA) &
-            ~utils.array_equals_nodata(s_bar, _TARGET_NODATA) &
-            ~utils.array_equals_nodata(
-                flow_accumulation, flow_accumulation_nodata))
-        d_up_array = numpy.empty(valid_mask.shape, dtype=numpy.float32)
-        d_up_array[:] = _TARGET_NODATA
-        d_up_array[valid_mask] = (
-            w_bar[valid_mask] * s_bar[valid_mask] * numpy.sqrt(
-                flow_accumulation[valid_mask] * cell_area))
-        return d_up_array
-
-    pygeoprocessing.raster_calculator(
-        [(path, 1) for path in [
-            w_bar_path, s_bar_path, flow_accumulation_path]], d_up_op,
-        out_d_up_path, gdal.GDT_Float32, _TARGET_NODATA)
+    pygeoprocessing.raster_map(
+        op=lambda w_bar, s_bar, flow_accum: (
+            w_bar * s_bar * numpy.sqrt(flow_accum * cell_area)),
+        rasters=[w_bar_path, s_bar_path, flow_accumulation_path],
+        target_path=out_d_up_path)
 
 
 def _calculate_d_up_bare(
@@ -1465,71 +1292,11 @@ def _calculate_d_up_bare(
     """Calculate s_bar * sqrt(flow accumulation * cell area)."""
     cell_area = abs(
         pygeoprocessing.get_raster_info(s_bar_path)['pixel_size'][0])**2
-    flow_accumulation_nodata = pygeoprocessing.get_raster_info(
-        flow_accumulation_path)['nodata'][0]
-
-    def d_up_op(s_bar, flow_accumulation):
-        """Calculate the bare d_up index.
-
-        s_bar * sqrt(upslope area)
-
-        """
-        valid_mask = (
-            ~utils.array_equals_nodata(
-                flow_accumulation, flow_accumulation_nodata) &
-            ~utils.array_equals_nodata(s_bar, _TARGET_NODATA))
-        d_up_array = numpy.empty(valid_mask.shape, dtype=numpy.float32)
-        d_up_array[:] = _TARGET_NODATA
-        d_up_array[valid_mask] = (
-            numpy.sqrt(flow_accumulation[valid_mask] * cell_area) *
-            s_bar[valid_mask])
-        return d_up_array
-
-    pygeoprocessing.raster_calculator(
-        [(s_bar_path, 1), (flow_accumulation_path, 1)], d_up_op,
-        out_d_up_bare_path, gdal.GDT_Float32, _TARGET_NODATA)
-
-
-def _calculate_inverse_ws_factor(
-        thresholded_slope_path, thresholded_w_factor_path,
-        out_ws_factor_inverse_path):
-    """Calculate 1/(w*s)."""
-    slope_nodata = pygeoprocessing.get_raster_info(
-        thresholded_slope_path)['nodata'][0]
-
-    def ws_op(w_factor, s_factor):
-        """Calculate the inverse ws factor."""
-        valid_mask = (
-            ~utils.array_equals_nodata(w_factor, _TARGET_NODATA) &
-            ~utils.array_equals_nodata(s_factor, slope_nodata))
-        result = numpy.empty(valid_mask.shape, dtype=numpy.float32)
-        result[:] = _TARGET_NODATA
-        result[valid_mask] = (
-            1.0 / (w_factor[valid_mask] * s_factor[valid_mask]))
-        return result
-
-    pygeoprocessing.raster_calculator(
-        [(thresholded_w_factor_path, 1), (thresholded_slope_path, 1)], ws_op,
-        out_ws_factor_inverse_path, gdal.GDT_Float32, _TARGET_NODATA)
-
-
-def _calculate_inverse_s_factor(
-        thresholded_slope_path, out_s_factor_inverse_path):
-    """Calculate 1/s."""
-    slope_nodata = pygeoprocessing.get_raster_info(
-        thresholded_slope_path)['nodata'][0]
-
-    def s_op(s_factor):
-        """Calculate the inverse s factor."""
-        valid_mask = ~utils.array_equals_nodata(s_factor, slope_nodata)
-        result = numpy.empty(valid_mask.shape, dtype=numpy.float32)
-        result[:] = _TARGET_NODATA
-        result[valid_mask] = 1.0 / s_factor[valid_mask]
-        return result
-
-    pygeoprocessing.raster_calculator(
-        [(thresholded_slope_path, 1)], s_op,
-        out_s_factor_inverse_path, gdal.GDT_Float32, _TARGET_NODATA)
+    pygeoprocessing.raster_map(
+        op=lambda s_bar, flow_accum: (
+            numpy.sqrt(flow_accum * cell_area) * s_bar),
+        rasters=[s_bar_path, flow_accumulation_path],
+        target_path=out_d_up_bare_path)
 
 
 def _calculate_ic(d_up_path, d_dn_path, out_ic_factor_path):
@@ -1540,8 +1307,8 @@ def _calculate_ic(d_up_path, d_dn_path, out_ic_factor_path):
     def ic_op(d_up, d_dn):
         """Calculate IC factor."""
         valid_mask = (
-            ~utils.array_equals_nodata(d_up, _TARGET_NODATA) &
-            ~utils.array_equals_nodata(d_dn, d_dn_nodata) & (d_dn != 0) &
+            ~pygeoprocessing.array_equals_nodata(d_up, _TARGET_NODATA) &
+            ~pygeoprocessing.array_equals_nodata(d_dn, d_dn_nodata) & (d_dn != 0) &
             (d_up != 0))
         ic_array = numpy.empty(valid_mask.shape, dtype=numpy.float32)
         ic_array[:] = _IC_NODATA
@@ -1560,7 +1327,7 @@ def _calculate_sdr(
     def sdr_op(ic_factor, stream):
         """Calculate SDR factor."""
         valid_mask = (
-            ~utils.array_equals_nodata(ic_factor, _IC_NODATA) & (stream != 1))
+            ~pygeoprocessing.array_equals_nodata(ic_factor, _IC_NODATA) & (stream != 1))
         result = numpy.empty(valid_mask.shape, dtype=numpy.float32)
         result[:] = _TARGET_NODATA
         result[valid_mask] = (
@@ -1573,30 +1340,13 @@ def sdr_op(ic_factor, stream):
         gdal.GDT_Float32, _TARGET_NODATA)
 
 
-def _calculate_sed_export(usle_path, sdr_path, target_sed_export_path):
-    """Calculate USLE * SDR."""
-    def sed_export_op(usle, sdr):
-        """Sediment export."""
-        valid_mask = (
-            ~utils.array_equals_nodata(usle, _TARGET_NODATA) &
-            ~utils.array_equals_nodata(sdr, _TARGET_NODATA))
-        result = numpy.empty(valid_mask.shape, dtype=numpy.float32)
-        result[:] = _TARGET_NODATA
-        result[valid_mask] = usle[valid_mask] * sdr[valid_mask]
-        return result
-
-    pygeoprocessing.raster_calculator(
-        [(usle_path, 1), (sdr_path, 1)], sed_export_op,
-        target_sed_export_path, gdal.GDT_Float32, _TARGET_NODATA)
-
-
 def _calculate_e_prime(usle_path, sdr_path, stream_path, target_e_prime):
     """Calculate USLE * (1-SDR)."""
     def e_prime_op(usle, sdr, streams):
         """Wash that does not reach stream."""
         valid_mask = (
-            ~utils.array_equals_nodata(usle, _TARGET_NODATA) &
-            ~utils.array_equals_nodata(sdr, _TARGET_NODATA))
+            ~pygeoprocessing.array_equals_nodata(usle, _TARGET_NODATA) &
+            ~pygeoprocessing.array_equals_nodata(sdr, _TARGET_NODATA))
         result = numpy.empty(valid_mask.shape, dtype=numpy.float32)
         result[:] = _TARGET_NODATA
         result[valid_mask] = usle[valid_mask] * (1-sdr[valid_mask])
diff --git a/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py b/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py
index 1d260110b2..e969be9e04 100644
--- a/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py
+++ b/src/natcap/invest/seasonal_water_yield/seasonal_water_yield.py
@@ -841,8 +841,11 @@ def execute(args):
             quick_flow_task_list.append(monthly_quick_flow_task)
 
         qf_task = task_graph.add_task(
-            func=_calculate_annual_qfi,
-            args=(file_registry['qfm_path_list'], file_registry['qf_path']),
+            func=pygeoprocessing.raster_map,
+            kwargs=dict(
+                op=qfi_sum_op,
+                rasters=file_registry['qfm_path_list'],
+                target_path=file_registry['qf_path']),
             target_path_list=[file_registry['qf_path']],
             dependent_task_list=quick_flow_task_list,
             task_name='calculate QFi')
@@ -962,6 +965,18 @@ def execute(args):
     LOGGER.info('      //_| //_|')
 
 
+# raster_map equation: sum the monthly qfis
+def qfi_sum_op(*qf_values): return numpy.sum(qf_values)
+
+
+def _calculate_l_avail(l_path, gamma, target_l_avail_path):
+    """l avail = l * gamma."""
+    pygeoprocessing.raster_map(
+        op=lambda l: numpy.min(numpy.stack((gamma * l, l)), axis=0),
+        rasters=[l_path],
+        target_path=target_l_avail_path)
+
+
 def _calculate_vri(l_path, target_vri_path):
     """Calculate VRI as li_array / qb_sum.
 
@@ -979,7 +994,7 @@ def _calculate_vri(l_path, target_vri_path):
 
     for _, block in pygeoprocessing.iterblocks((l_path, 1)):
         valid_mask = (
-            ~utils.array_equals_nodata(block, l_nodata) &
+            ~pygeoprocessing.array_equals_nodata(block, l_nodata) &
             (~numpy.isinf(block)))
         qb_sum += numpy.sum(block[valid_mask])
         qb_valid_count += numpy.count_nonzero(valid_mask)
@@ -990,50 +1005,19 @@ def vri_op(li_array):
         result = numpy.empty_like(li_array)
         result[:] = li_nodata
         if qb_sum > 0:
-            valid_mask = ~utils.array_equals_nodata(li_array, li_nodata)
+            valid_mask = ~pygeoprocessing.array_equals_nodata(li_array, li_nodata)
             try:
                 result[valid_mask] = li_array[valid_mask] / qb_sum
             except RuntimeWarning:
                 LOGGER.exception(qb_sum)
                 raise
         return result
+
     pygeoprocessing.raster_calculator(
         [(l_path, 1)], vri_op, target_vri_path, gdal.GDT_Float32,
         li_nodata)
 
 
-def _calculate_annual_qfi(qfm_path_list, target_qf_path):
-    """Calculate annual quickflow.
-
-    Args:
-        qfm_path_list (list): list of monthly quickflow raster paths.
-        target_qf_path (str): path to target annual quickflow raster.
-
-    Returns:
-        None.
-
-    """
-    qf_nodata = -1
-
-    def qfi_sum_op(*qf_values):
-        """Sum the monthly qfis."""
-
-        # only calculate the sum where data is available for all 12 months
-        valid_mask = numpy.full(qf_values[0].shape, True)
-        for qf_array in qf_values:
-            valid_mask &= ~utils.array_equals_nodata(qf_array, qf_nodata)
-
-        qf_sum = numpy.full(qf_values[0].shape, qf_nodata, dtype=numpy.float32)
-        qf_sum[valid_mask] = 0
-        for qf_array in qf_values:
-            qf_sum[valid_mask] += qf_array[valid_mask]
-        return qf_sum
-
-    pygeoprocessing.raster_calculator(
-        [(path, 1) for path in qfm_path_list],
-        qfi_sum_op, target_qf_path, gdal.GDT_Float32, qf_nodata)
-
-
 def _calculate_monthly_quick_flow(precip_path, n_events_path, stream_path,
         si_path, qf_monthly_path):
     """Calculate quick flow for a month.
@@ -1068,8 +1052,8 @@ def qf_op(p_im, s_i, n_m, stream):
         Returns:
             quick flow (numpy.array)
         """
-        valid_p_mask = ~utils.array_equals_nodata(p_im, p_nodata)
-        valid_n_mask = ~utils.array_equals_nodata(n_m, n_nodata)
+        valid_p_mask = ~pygeoprocessing.array_equals_nodata(p_im, p_nodata)
+        valid_n_mask = ~pygeoprocessing.array_equals_nodata(n_m, n_nodata)
         # precip mask: both p_im and n_m are defined and greater than 0
         precip_mask = valid_p_mask & valid_n_mask & (p_im > 0) & (n_m > 0)
         stream_mask = stream == 1
@@ -1078,8 +1062,8 @@ def qf_op(p_im, s_i, n_m, stream):
         valid_mask = (
           valid_p_mask &
           valid_n_mask &
-          ~utils.array_equals_nodata(stream, stream_nodata) &
-          ~utils.array_equals_nodata(s_i, si_nodata))
+          ~pygeoprocessing.array_equals_nodata(stream, stream_nodata) &
+          ~pygeoprocessing.array_equals_nodata(s_i, si_nodata))
 
         # QF is defined in terms of three cases:
         #
@@ -1203,23 +1187,14 @@ def _calculate_curve_number_raster(
     Returns:
         None
     """
-    soil_nodata = pygeoprocessing.get_raster_info(
-        soil_group_path)['nodata'][0]
     map_soil_type_to_header = {
         1: 'cn_a',
         2: 'cn_b',
         3: 'cn_c',
         4: 'cn_d',
     }
-    # curve numbers are always positive so -1 a good nodata choice
     lulc_to_soil = {}
-    lulc_nodata = pygeoprocessing.get_raster_info(
-        lulc_raster_path)['nodata'][0]
-
     lucodes = biophysical_df.index.to_list()
-    if lulc_nodata is not None:
-        lucodes.append(lulc_nodata)
-
     for soil_id, soil_column in map_soil_type_to_header.items():
         lulc_to_soil[soil_id] = {
             'lulc_values': [],
@@ -1227,14 +1202,9 @@ def _calculate_curve_number_raster(
         }
 
         for lucode in sorted(lucodes):
-            if lucode != lulc_nodata:
-                lulc_to_soil[soil_id]['cn_values'].append(
-                    biophysical_df[soil_column][lucode])
-                lulc_to_soil[soil_id]['lulc_values'].append(lucode)
-            else:
-                # handle the lulc nodata with cn nodata
-                lulc_to_soil[soil_id]['lulc_values'].append(lulc_nodata)
-                lulc_to_soil[soil_id]['cn_values'].append(TARGET_NODATA)
+            lulc_to_soil[soil_id]['cn_values'].append(
+                biophysical_df[soil_column][lucode])
+            lulc_to_soil[soil_id]['lulc_values'].append(lucode)
 
         # Making the landcover array a float32 in case the user provides a
         # float landcover map like Kate did.
@@ -1247,7 +1217,7 @@ def _calculate_curve_number_raster(
 
     # Use set of table lucodes in cn_op
     lucodes_set = set(lucodes)
-    valid_soil_groups = set([soil_nodata, *map_soil_type_to_header.keys()])
+    valid_soil_groups = set(map_soil_type_to_header.keys())
 
     def cn_op(lulc_array, soil_group_array):
         """Map lulc code and soil to a curve number."""
@@ -1276,11 +1246,10 @@ def cn_op(lulc_array, soil_group_array):
                 "The soil group raster must only have groups 1, 2, 3 or 4. "
                 f"Invalid group(s) {', '.join(invalid_soil_groups)} were "
                 f"found in soil group raster {soil_group_path} "
-                f"(nodata value: {soil_nodata})")
+                "(nodata value: "
+                f"{pygeoprocessing.get_raster_info(soil_group_path)['nodata'][0]})")
 
         for soil_group_id in unique_soil_groups:
-            if soil_group_id == soil_nodata:
-                continue
             current_soil_mask = (soil_group_array == soil_group_id)
             index = numpy.digitize(
                 lulc_array.ravel(),
@@ -1291,9 +1260,10 @@ def cn_op(lulc_array, soil_group_array):
             cn_result[current_soil_mask] = cn_values[current_soil_mask]
         return cn_result
 
-    pygeoprocessing.raster_calculator(
-        [(lulc_raster_path, 1), (soil_group_path, 1)], cn_op, cn_path,
-        gdal.GDT_Float32, TARGET_NODATA)
+    pygeoprocessing.raster_map(
+        op=cn_op,
+        rasters=[lulc_raster_path, soil_group_path],
+        target_path=cn_path)
 
 
 def _calculate_si_raster(cn_path, stream_path, si_path):
@@ -1312,7 +1282,7 @@ def _calculate_si_raster(cn_path, stream_path, si_path):
     def si_op(ci_factor, stream_mask):
         """Calculate si factor."""
         valid_mask = (
-            ~utils.array_equals_nodata(ci_factor, cn_nodata) &
+            ~pygeoprocessing.array_equals_nodata(ci_factor, cn_nodata) &
             (ci_factor > 0))
         si_array = numpy.empty(ci_factor.shape)
         si_array[:] = TARGET_NODATA
@@ -1399,24 +1369,6 @@ def _aggregate_recharge(
     aggregate_vector = None
 
 
-def _calculate_l_avail(l_path, gamma, target_l_avail_path):
-    """l avail = l * gamma."""
-    li_nodata = pygeoprocessing.get_raster_info(l_path)['nodata'][0]
-
-    def l_avail_op(l_array):
-        """Calculate equation [8] L_avail = min(gamma*L, L)."""
-        result = numpy.empty(l_array.shape)
-        result[:] = li_nodata
-        valid_mask = ~utils.array_equals_nodata(l_array, li_nodata)
-        result[valid_mask] = numpy.min(numpy.stack(
-            (gamma*l_array[valid_mask], l_array[valid_mask])), axis=0)
-        return result
-
-    pygeoprocessing.raster_calculator(
-        [(l_path, 1)], l_avail_op, target_l_avail_path, gdal.GDT_Float32,
-        li_nodata)
-
-
 @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/stormwater.py b/src/natcap/invest/stormwater.py
index a765b327dd..5bd7f53664 100644
--- a/src/natcap/invest/stormwater.py
+++ b/src/natcap/invest/stormwater.py
@@ -5,6 +5,7 @@
 
 import numpy
 import pygeoprocessing
+import pygeoprocessing.kernels
 import taskgraph
 from osgeo import gdal
 from osgeo import ogr
@@ -625,16 +626,15 @@ def execute(args):
         # Using the averaged retention ratio raster and boolean
         # "within radius" rasters, adjust the retention ratios
         adjust_retention_ratio_task = task_graph.add_task(
-            func=pygeoprocessing.raster_calculator,
-            args=([
-                (files['retention_ratio_path'], 1),
-                (files['ratio_average_path'], 1),
-                (files['near_connected_lulc_path'], 1),
-                (files['near_road_path'], 1)],
-                adjust_op,
-                files['adjusted_retention_ratio_path'],
-                gdal.GDT_Float32,
-                FLOAT_NODATA),
+            func=pygeoprocessing.raster_map,
+            kwargs=dict(
+                op=adjust_op,
+                rasters=[
+                    files['retention_ratio_path'],
+                    files['ratio_average_path'],
+                    files['near_connected_lulc_path'],
+                    files['near_road_path']],
+                target_path=files['adjusted_retention_ratio_path']),
             target_path_list=[files['adjusted_retention_ratio_path']],
             task_name='adjust stormwater retention ratio',
             dependent_task_list=[retention_ratio_task, average_ratios_task,
@@ -665,13 +665,11 @@ def execute(args):
 
     # Calculate stormwater runoff ratios and volume
     runoff_ratio_task = task_graph.add_task(
-        func=pygeoprocessing.raster_calculator,
-        args=(
-            [(final_retention_ratio_path, 1)],
-            retention_to_runoff_op,
-            files['runoff_ratio_path'],
-            gdal.GDT_Float32,
-            FLOAT_NODATA),
+        func=pygeoprocessing.raster_map,
+        kwargs=dict(
+            op=retention_to_runoff_op,
+            rasters=[final_retention_ratio_path],
+            target_path=files['runoff_ratio_path']),
         target_path_list=[files['runoff_ratio_path']],
         dependent_task_list=[final_retention_ratio_task],
         task_name='calculate stormwater runoff ratio'
@@ -862,37 +860,19 @@ def lookup_ratios(lulc_path, soil_group_path, ratio_lookup, sorted_lucodes,
     Returns:
         None
     """
-    lulc_nodata = pygeoprocessing.get_raster_info(lulc_path)['nodata'][0]
-    soil_group_nodata = pygeoprocessing.get_raster_info(
-        soil_group_path)['nodata'][0]
     # insert a column on the left side of the array so that the soil
     # group codes 1-4 line up with their indexes. this is faster than
     # decrementing every value in a large raster.
     ratio_lookup = numpy.insert(ratio_lookup, 0,
                                 numpy.zeros(ratio_lookup.shape[0]), axis=1)
-
-    def ratio_op(lulc_array, soil_group_array):
-        output_ratio_array = numpy.full(lulc_array.shape, FLOAT_NODATA,
-                                        dtype=numpy.float32)
-        valid_mask = numpy.full(lulc_array.shape, True)
-        if lulc_nodata is not None:
-            valid_mask &= ~utils.array_equals_nodata(lulc_array, lulc_nodata)
-        if soil_group_nodata is not None:
-            valid_mask &= ~utils.array_equals_nodata(
-                soil_group_array, soil_group_nodata)
-        # the index of each lucode in the sorted lucodes array
-        lulc_index = numpy.digitize(lulc_array[valid_mask], sorted_lucodes,
-                                    right=True)
-        output_ratio_array[valid_mask] = (
-            ratio_lookup[lulc_index, soil_group_array[valid_mask]])
-        return output_ratio_array
-
-    pygeoprocessing.raster_calculator(
-        [(lulc_path, 1), (soil_group_path, 1)],
-        ratio_op,
-        output_path,
-        gdal.GDT_Float32,
-        FLOAT_NODATA)
+    pygeoprocessing.raster_map(
+        op=lambda lulc_array, soil_group_array: ratio_lookup[
+            # the index of each lucode in the sorted lucodes array
+            numpy.digitize(lulc_array, sorted_lucodes, right=True),
+            soil_group_array],
+        rasters=[lulc_path, soil_group_path],
+        target_path=output_path,
+        target_dtype=numpy.float32)
 
 
 def volume_op(ratio_array, precip_array, precip_nodata, pixel_area):
@@ -911,9 +891,9 @@ def volume_op(ratio_array, precip_array, precip_nodata, pixel_area):
     """
     volume_array = numpy.full(ratio_array.shape, FLOAT_NODATA,
                               dtype=numpy.float32)
-    valid_mask = ~utils.array_equals_nodata(ratio_array, FLOAT_NODATA)
+    valid_mask = ~pygeoprocessing.array_equals_nodata(ratio_array, FLOAT_NODATA)
     if precip_nodata is not None:
-        valid_mask &= ~utils.array_equals_nodata(precip_array, precip_nodata)
+        valid_mask &= ~pygeoprocessing.array_equals_nodata(precip_array, precip_nodata)
 
     # precipitation (mm/yr) * pixel area (m^2) *
     # 0.001 (m/mm) * ratio = volume (m^3/yr)
@@ -934,11 +914,7 @@ def retention_to_runoff_op(retention_array):
     Returns:
         numpy.ndarray of stormwater runoff ratios
     """
-    runoff_array = numpy.full(retention_array.shape, FLOAT_NODATA,
-                              dtype=numpy.float32)
-    valid_mask = ~utils.array_equals_nodata(retention_array, FLOAT_NODATA)
-    runoff_array[valid_mask] = 1 - retention_array[valid_mask]
-    return runoff_array
+    return 1 - retention_array
 
 
 def pollutant_load_op(lulc_array, lulc_nodata, volume_array, sorted_lucodes,
@@ -968,9 +944,9 @@ def pollutant_load_op(lulc_array, lulc_nodata, volume_array, sorted_lucodes,
     """
     load_array = numpy.full(
         lulc_array.shape, FLOAT_NODATA, dtype=numpy.float32)
-    valid_mask = ~utils.array_equals_nodata(volume_array, FLOAT_NODATA)
+    valid_mask = ~pygeoprocessing.array_equals_nodata(volume_array, FLOAT_NODATA)
     if lulc_nodata is not None:
-        valid_mask &= ~utils.array_equals_nodata(lulc_array, lulc_nodata)
+        valid_mask &= ~pygeoprocessing.array_equals_nodata(lulc_array, lulc_nodata)
 
     # bin each value in the LULC array such that
     # lulc_array[i,j] == sorted_lucodes[lulc_index[i,j]]. thus,
@@ -1000,7 +976,7 @@ def retention_value_op(retention_volume_array, replacement_cost):
     """
     value_array = numpy.full(retention_volume_array.shape, FLOAT_NODATA,
                              dtype=numpy.float32)
-    valid_mask = ~utils.array_equals_nodata(
+    valid_mask = ~pygeoprocessing.array_equals_nodata(
         retention_volume_array, FLOAT_NODATA)
 
     # retention (m^3/yr) * replacement cost ($/m^3) = retention value ($/yr)
@@ -1013,8 +989,8 @@ def adjust_op(ratio_array, avg_ratio_array, near_connected_lulc_array,
               near_road_array):
     """Apply the retention ratio adjustment algorithm to an array of ratios.
 
-    This is meant to be used with raster_calculator. Assumes that the nodata
-    value for all four input arrays is the global FLOAT_NODATA.
+    This is meant to be used with raster_map. Assumes that nodata is already
+    filtered out.
 
     Args:
         ratio_array (numpy.ndarray): 2D array of stormwater retention ratios
@@ -1028,36 +1004,14 @@ def adjust_op(ratio_array, avg_ratio_array, near_connected_lulc_array,
         2D numpy array of adjusted retention ratios. Has the same shape as
         ``retention_ratio_array``.
     """
-    adjusted_ratio_array = numpy.full(ratio_array.shape, FLOAT_NODATA,
-                                      dtype=numpy.float32)
-    adjustment_factor_array = numpy.full(ratio_array.shape, FLOAT_NODATA,
-                                         dtype=numpy.float32)
-    valid_mask = (
-        ~utils.array_equals_nodata(ratio_array, FLOAT_NODATA) &
-        ~utils.array_equals_nodata(avg_ratio_array, FLOAT_NODATA) &
-        (near_connected_lulc_array != UINT8_NODATA) &
-        (near_road_array != UINT8_NODATA))
-
     # adjustment factor:
     # - 0 if any of the nearby pixels are impervious/connected;
     # - average of nearby pixels, otherwise
-    is_not_connected = ~(
-        near_connected_lulc_array[valid_mask] |
-        near_road_array[valid_mask]).astype(bool)
-    adjustment_factor_array[valid_mask] = (avg_ratio_array[valid_mask] *
-                                           is_not_connected)
-
-    adjustment_factor_array[valid_mask] = (
-        avg_ratio_array[valid_mask] * ~(
-            near_connected_lulc_array[valid_mask] |
-            near_road_array[valid_mask]
-        ).astype(bool))
-
     # equation 2-4: Radj_ij = R_ij + (1 - R_ij) * C_ij
-    adjusted_ratio_array[valid_mask] = (
-        ratio_array[valid_mask] +
-        (1 - ratio_array[valid_mask]) * adjustment_factor_array[valid_mask])
-    return adjusted_ratio_array
+    return ratio_array + (1 - ratio_array) * (
+        avg_ratio_array * ~(
+            near_connected_lulc_array | near_road_array
+        ).astype(bool))
 
 
 def aggregate_results(base_aggregate_areas_path, target_vector_path, srs_wkt,
@@ -1140,60 +1094,15 @@ def is_near(input_path, radius, distance_path, out_path):
         None
     """
     # Calculate the distance from each pixel to the nearest '1' pixel
-    pygeoprocessing.distance_transform_edt(
-        (input_path, 1),
-        distance_path)
-
-    def lte_threshold_op(array, threshold):
-        """Binary array of elements less than or equal to the threshold."""
-        # no need to mask nodata because distance_transform_edt doesn't
-        # output any nodata pixels
-        return array <= threshold
+    pygeoprocessing.distance_transform_edt((input_path, 1), distance_path)
 
     # Threshold that to a binary array so '1' means it's within the radius
-    pygeoprocessing.raster_calculator(
-        [(distance_path, 1), (radius, 'raw')],
-        lte_threshold_op,
-        out_path,
-        gdal.GDT_Byte,
-        UINT8_NODATA)
-
-
-def make_search_kernel(raster_path, radius):
-    """Make a search kernel for a raster that marks pixels within a radius.
-
-    Args:
-        raster_path (str): path to a raster to make kernel for. It is assumed
-            that the raster has square pixels.
-        radius (float): distance around each pixel's centerpoint to search
-            in raster coordinate system units
-
-    Returns:
-        2D boolean numpy.ndarray. '1' pixels are within ``radius`` of the
-        center pixel, measured centerpoint-to-centerpoint. '0' pixels are
-        outside the radius. The array dimensions are as small as possible
-        while still including the entire radius.
-    """
-    raster_info = pygeoprocessing.get_raster_info(raster_path)
-    pixel_radius = radius / abs(raster_info['pixel_size'][0])
-    pixel_margin = math.floor(pixel_radius)
-    # the search kernel is just large enough to contain all pixels that
-    # *could* be within the radius of the center pixel
-    search_kernel_shape = tuple([pixel_margin * 2 + 1] * 2)
-    # arrays of the column index and row index of each pixel
-    col_indices, row_indices = numpy.indices(search_kernel_shape)
-    # adjust them so that (0, 0) is the center pixel
-    col_indices -= pixel_margin
-    row_indices -= pixel_margin
-    # hypotenuse_i = sqrt(col_indices_i**2 + row_indices_i**2) for each pixel i
-    hypotenuse = numpy.hypot(col_indices, row_indices)
-    # boolean kernel where 1=pixel centerpoint is within the radius of the
-    # center pixel's centerpoint
-    search_kernel = numpy.array(hypotenuse <= pixel_radius, dtype=numpy.uint8)
-    LOGGER.debug(
-        f'Search kernel for {raster_path} with radius {radius}:'
-        f'\n{search_kernel}')
-    return search_kernel
+    pygeoprocessing.raster_map(
+        op=lambda dist: dist <= radius,
+        rasters=[distance_path],
+        target_path=out_path,
+        target_dtype=numpy.uint8,
+        target_nodata=UINT8_NODATA)
 
 
 def raster_average(raster_path, radius, kernel_path, out_path):
@@ -1223,19 +1132,12 @@ def raster_average(raster_path, radius, kernel_path, out_path):
     Returns:
         None
     """
-    search_kernel = make_search_kernel(raster_path, radius)
-
-    srs = osr.SpatialReference()
-    srs.ImportFromEPSG(3857)
-    projection_wkt = srs.ExportToWkt()
-    pygeoprocessing.numpy_array_to_raster(
-        # float32 here to avoid pygeoprocessing bug issue #180
-        search_kernel.astype(numpy.float32),
-        FLOAT_NODATA,
-        (20, -20),
-        (0, 0),
-        projection_wkt,
-        kernel_path)
+    pixel_radius = radius / abs(pygeoprocessing.get_raster_info(
+        raster_path)['pixel_size'][0])
+    pygeoprocessing.kernels.dichotomous_kernel(
+        target_kernel_path=kernel_path,
+        max_distance=pixel_radius,
+        normalize=False)
 
     # convolve the signal (input raster) with the kernel and normalize
     # this is equivalent to taking an average of each pixel's neighborhood
diff --git a/src/natcap/invest/urban_cooling_model.py b/src/natcap/invest/urban_cooling_model.py
index c40a0593b0..1f1f86b535 100644
--- a/src/natcap/invest/urban_cooling_model.py
+++ b/src/natcap/invest/urban_cooling_model.py
@@ -9,6 +9,7 @@
 
 import numpy
 import pygeoprocessing
+import pygeoprocessing.kernels
 import rtree
 import shapely.prepared
 import shapely.wkb
@@ -530,8 +531,11 @@ def execute(args):
     area_kernel_path = os.path.join(
         intermediate_dir, f'area_kernel{file_suffix}.tif')
     area_kernel_task = task_graph.add_task(
-        func=flat_disk_kernel,
-        args=(green_area_decay_kernel_distance, area_kernel_path),
+        func=pygeoprocessing.kernels.dichotomous_kernel,
+        kwargs=dict(
+            target_kernel_path=area_kernel_path,
+            max_distance=green_area_decay_kernel_distance,
+            normalize=False),
         target_path_list=[area_kernel_path],
         task_name='area kernel')
 
@@ -604,10 +608,11 @@ def execute(args):
         LOGGER.info('Calculating Cooling Coefficient using '
                     'building intensity')
         cc_task = task_graph.add_task(
-            func=pygeoprocessing.raster_calculator,
-            args=([(task_path_prop_map['building_intensity'][1], 1)],
-                  calc_cc_op_intensity, cc_raster_path,
-                  gdal.GDT_Float32, TARGET_NODATA),
+            func=pygeoprocessing.raster_map,
+            kwargs=dict(
+                op=calc_cc_op_intensity,
+                rasters=[task_path_prop_map['building_intensity'][1]],
+                target_path=cc_raster_path),
             target_path_list=[cc_raster_path],
             dependent_task_list=[
                 task_path_prop_map['building_intensity'][0]],
@@ -1184,7 +1189,7 @@ def calc_t_air_nomix_op(t_ref_val, hm_array, uhi_max):
     result = numpy.empty(hm_array.shape, dtype=numpy.float32)
     result[:] = TARGET_NODATA
     # TARGET_NODATA should never be None
-    valid_mask = ~utils.array_equals_nodata(hm_array, TARGET_NODATA)
+    valid_mask = ~pygeoprocessing.array_equals_nodata(hm_array, TARGET_NODATA)
     result[valid_mask] = t_ref_val + (1-hm_array[valid_mask]) * uhi_max
     return result
 
@@ -1212,9 +1217,9 @@ def calc_cc_op_factors(
     result = numpy.empty(shade_array.shape, dtype=numpy.float32)
     result[:] = TARGET_NODATA
     valid_mask = ~(
-        utils.array_equals_nodata(shade_array, TARGET_NODATA) |
-        utils.array_equals_nodata(albedo_array, TARGET_NODATA) |
-        utils.array_equals_nodata(eti_array, TARGET_NODATA))
+        pygeoprocessing.array_equals_nodata(shade_array, TARGET_NODATA) |
+        pygeoprocessing.array_equals_nodata(albedo_array, TARGET_NODATA) |
+        pygeoprocessing.array_equals_nodata(eti_array, TARGET_NODATA))
     result[valid_mask] = (
         cc_weight_shade*shade_array[valid_mask] +
         cc_weight_albedo*albedo_array[valid_mask] +
@@ -1232,11 +1237,7 @@ def calc_cc_op_intensity(intensity_array):
         A numpy array of ``1 - intensity_array``.
 
     """
-    result = numpy.empty(intensity_array.shape, dtype=numpy.float32)
-    result[:] = TARGET_NODATA
-    valid_mask = ~utils.array_equals_nodata(intensity_array, TARGET_NODATA)
-    result[valid_mask] = 1 - intensity_array[valid_mask]
-    return result
+    return 1 - intensity_array
 
 
 def calc_eti_op(
@@ -1245,9 +1246,9 @@ def calc_eti_op(
     result = numpy.empty(kc_array.shape, dtype=numpy.float32)
     result[:] = target_nodata
     # kc intermediate output should always have a nodata value defined
-    valid_mask = ~utils.array_equals_nodata(kc_array, kc_nodata)
+    valid_mask = ~pygeoprocessing.array_equals_nodata(kc_array, kc_nodata)
     if et0_nodata is not None:
-        valid_mask &= ~utils.array_equals_nodata(et0_array, et0_nodata)
+        valid_mask &= ~pygeoprocessing.array_equals_nodata(et0_array, et0_nodata)
     result[valid_mask] = (
         kc_array[valid_mask] * et0_array[valid_mask] / et_max)
     return result
@@ -1271,106 +1272,12 @@ def calculate_wbgt(
 
     """
     LOGGER.info('Calculating WBGT')
-    t_air_nodata = pygeoprocessing.get_raster_info(
-        t_air_raster_path)['nodata'][0]
-
-    def wbgt_op(avg_rel_humidity, t_air_array):
-        wbgt = numpy.empty(t_air_array.shape, dtype=numpy.float32)
-
-        valid_mask = slice(None)
-        if t_air_nodata is not None:
-            valid_mask = ~utils.array_equals_nodata(t_air_array, t_air_nodata)
-        wbgt[:] = TARGET_NODATA
-        t_air_valid = t_air_array[valid_mask]
-        e_i = (
+    pygeoprocessing.raster_map(
+        op=lambda t_air: 0.567 * t_air + 0.393 * (
             (avg_rel_humidity / 100) * 6.105 * numpy.exp(
-                17.27 * (t_air_valid / (237.7 + t_air_valid))))
-        wbgt[valid_mask] = 0.567 * t_air_valid + 0.393 * e_i + 3.94
-        return wbgt
-
-    pygeoprocessing.raster_calculator(
-        [(avg_rel_humidity, 'raw'), (t_air_raster_path, 1)],
-        wbgt_op, target_vapor_pressure_path, gdal.GDT_Float32,
-        TARGET_NODATA)
-
-
-def flat_disk_kernel(max_distance, kernel_filepath):
-    """Create a flat disk  kernel.
-
-    The raster created will be a tiled GeoTiff, with 256x256 memory blocks.
-
-    Args:
-        max_distance (int): The distance (in pixels) of the
-            kernel's radius.
-        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
-
-    """
-    LOGGER.info(f'Creating a disk kernel of distance {max_distance} at '
-                f'{kernel_filepath}')
-    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_Byte, options=[
-            'BIGTIFF=IF_SAFER', 'TILED=YES', 'BLOCKXSIZE=256',
-            'BLOCKYSIZE=256'])
-
-    # Make some kind of geotransform and SRS. It doesn't matter what, but
-    # having one 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(255)
-
-    cols_per_block, rows_per_block = kernel_band.GetBlockSize()
-
-    n_cols = kernel_dataset.RasterXSize
-    n_rows = kernel_dataset.RasterYSize
-
-    n_col_blocks = int(math.ceil(n_cols / float(cols_per_block)))
-    n_row_blocks = int(math.ceil(n_rows / float(rows_per_block)))
-
-    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
-
-        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
-
-            # 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)
-
-            kernel_index_distances = numpy.hypot(
-                row_indices, col_indices)
-            kernel = kernel_index_distances < max_distance
-
-            kernel_band.WriteArray(kernel, xoff=col_offset,
-                                   yoff=row_offset)
-
-    # Need to flush the kernel's cache to disk before opening up a new Dataset
-    # object in interblocks()
-    kernel_dataset.FlushCache()
+                17.27 * (t_air / (237.7 + t_air)))) + 3.94,
+        rasters=[t_air_raster_path],
+        target_path=target_vapor_pressure_path)
 
 
 def hm_op(cc_array, green_area_sum, cc_park_array, green_area_threshold):
@@ -1392,8 +1299,8 @@ def hm_op(cc_array, green_area_sum, cc_park_array, green_area_threshold):
     """
     result = numpy.empty(cc_array.shape, dtype=numpy.float32)
     result[:] = TARGET_NODATA
-    valid_mask = ~(utils.array_equals_nodata(cc_array, TARGET_NODATA) &
-                   utils.array_equals_nodata(cc_park_array, TARGET_NODATA))
+    valid_mask = ~(pygeoprocessing.array_equals_nodata(cc_array, TARGET_NODATA) &
+                   pygeoprocessing.array_equals_nodata(cc_park_array, TARGET_NODATA))
     cc_mask = ((cc_array >= cc_park_array) |
                (green_area_sum < green_area_threshold))
     result[cc_mask & valid_mask] = cc_array[cc_mask & valid_mask]
@@ -1420,33 +1327,24 @@ def map_work_loss(
     """
     LOGGER.info(
         f'Calculating work loss using thresholds: {work_temp_threshold_array}')
-    byte_target_nodata = 255
 
     def classify_to_percent_op(temperature_array):
         result = numpy.empty(temperature_array.shape)
-        result[:] = byte_target_nodata
-        valid_mask = ~utils.array_equals_nodata(
-            temperature_array, TARGET_NODATA)
-        result[
-            valid_mask &
-            (temperature_array < work_temp_threshold_array[0])] = 0
+        result[temperature_array < work_temp_threshold_array[0]] = 0
         result[
-            valid_mask &
             (temperature_array >= work_temp_threshold_array[0]) &
             (temperature_array < work_temp_threshold_array[1])] = 25
         result[
-            valid_mask &
             (temperature_array >= work_temp_threshold_array[1]) &
             (temperature_array < work_temp_threshold_array[2])] = 50
-        result[
-            valid_mask &
-            (temperature_array >= work_temp_threshold_array[2])] = 75
+        result[temperature_array >= work_temp_threshold_array[2]] = 75
         return result
 
-    pygeoprocessing.raster_calculator(
-        [(temperature_raster_path, 1)], classify_to_percent_op,
-        work_loss_raster_path, gdal.GDT_Byte,
-        nodata_target=byte_target_nodata)
+    pygeoprocessing.raster_map(
+        op=classify_to_percent_op,
+        rasters=[temperature_raster_path],
+        target_path=work_loss_raster_path,
+        target_dtype=numpy.uint8)
 
 
 def _invoke_timed_callback(
@@ -1496,8 +1394,10 @@ def convolve_2d_by_exponential(
         dir=os.path.dirname(target_convolve_raster_path))
     exponential_kernel_path = os.path.join(
         temporary_working_dir, 'exponential_decay_kernel.tif')
-    utils.exponential_decay_kernel_raster(
-        decay_kernel_distance, exponential_kernel_path)
+    pygeoprocessing.kernels.exponential_decay_kernel(
+        target_kernel_path=exponential_kernel_path,
+        max_distance=decay_kernel_distance * 5,
+        expected_distance=decay_kernel_distance)
     pygeoprocessing.convolve_2d(
         (signal_raster_path, 1), (exponential_kernel_path, 1),
         target_convolve_raster_path, working_dir=temporary_working_dir,
diff --git a/src/natcap/invest/urban_flood_risk_mitigation.py b/src/natcap/invest/urban_flood_risk_mitigation.py
index f543b374d6..7703dfcd1e 100644
--- a/src/natcap/invest/urban_flood_risk_mitigation.py
+++ b/src/natcap/invest/urban_flood_risk_mitigation.py
@@ -722,7 +722,7 @@ def _flood_vol_op(
     """
     result = numpy.empty(q_pi_array.shape, dtype=numpy.float32)
     result[:] = target_nodata
-    valid_mask = ~utils.array_equals_nodata(q_pi_array, q_pi_nodata)
+    valid_mask = ~pygeoprocessing.array_equals_nodata(q_pi_array, q_pi_nodata)
     # 0.001 converts mm (quickflow) to m (pixel area units)
     result[valid_mask] = (
         q_pi_array[valid_mask] * pixel_area * 0.001)
@@ -747,7 +747,7 @@ def _runoff_retention_vol_op(
     """
     result = numpy.empty(runoff_retention_array.shape, dtype=numpy.float32)
     result[:] = target_nodata
-    valid_mask = ~utils.array_equals_nodata(
+    valid_mask = ~pygeoprocessing.array_equals_nodata(
         runoff_retention_array, runoff_retention_nodata)
     # the 1e-3 converts the mm of p_value to meters.
     result[valid_mask] = (
@@ -773,7 +773,7 @@ def _runoff_retention_op(q_pi_array, p_value, q_pi_nodata, result_nodata):
     result[:] = result_nodata
     valid_mask = numpy.ones(q_pi_array.shape, dtype=bool)
     if q_pi_nodata is not None:
-        valid_mask[:] = ~utils.array_equals_nodata(q_pi_array, q_pi_nodata)
+        valid_mask[:] = ~pygeoprocessing.array_equals_nodata(q_pi_array, q_pi_nodata)
     result[valid_mask] = 1 - (q_pi_array[valid_mask] / p_value)
     return result
 
@@ -798,7 +798,7 @@ def _q_pi_op(p_value, s_max_array, s_max_nodata, result_nodata):
     zero_mask = (p_value <= lam * s_max_array)
     non_nodata_mask = numpy.ones(s_max_array.shape, dtype=bool)
     if s_max_nodata is not None:
-        non_nodata_mask[:] = ~utils.array_equals_nodata(
+        non_nodata_mask[:] = ~pygeoprocessing.array_equals_nodata(
             s_max_array, s_max_nodata)
 
     # valid if not nodata and not going to be set to 0.
@@ -828,7 +828,7 @@ def _s_max_op(cn_array, cn_nodata, result_nodata):
     zero_mask = cn_array == 0
     valid_mask = ~zero_mask
     if cn_nodata is not None:
-        valid_mask[:] &= ~utils.array_equals_nodata(cn_array, cn_nodata)
+        valid_mask[:] &= ~pygeoprocessing.array_equals_nodata(cn_array, cn_nodata)
     result[valid_mask] = 25400 / cn_array[valid_mask] - 254
     result[zero_mask] = 0
     return result
@@ -856,10 +856,10 @@ def _lu_to_cn_op(
     result[:] = cn_nodata
     valid_mask = numpy.ones(lucode_array.shape, dtype=bool)
     if lucode_nodata is not None:
-        valid_mask[:] &= ~utils.array_equals_nodata(
+        valid_mask[:] &= ~pygeoprocessing.array_equals_nodata(
             lucode_array, lucode_nodata)
     if soil_type_nodata is not None:
-        valid_mask[:] &= ~utils.array_equals_nodata(
+        valid_mask[:] &= ~pygeoprocessing.array_equals_nodata(
             soil_type_array, soil_type_nodata)
 
     # this is an array where each column represents a valid landcover
diff --git a/src/natcap/invest/urban_nature_access.py b/src/natcap/invest/urban_nature_access.py
index 99a9d26cf5..f02a443618 100644
--- a/src/natcap/invest/urban_nature_access.py
+++ b/src/natcap/invest/urban_nature_access.py
@@ -9,6 +9,7 @@
 import numpy
 import numpy.testing
 import pygeoprocessing
+import pygeoprocessing.kernels
 import pygeoprocessing.symbolic
 import shapely.ops
 import shapely.wkb
@@ -22,7 +23,6 @@
 from . import utils
 from . import validation
 from .model_metadata import MODEL_METADATA
-from .ndr import ndr
 from .spec_utils import u
 
 LOGGER = logging.getLogger(__name__)
@@ -717,8 +717,6 @@ def execute(args):
         os.path.join(args['workspace_dir'], 'taskgraph_cache'), n_workers)
 
     kernel_creation_functions = {
-        KERNEL_LABEL_DICHOTOMY: _kernel_dichotomy,
-        KERNEL_LABEL_EXPONENTIAL: _kernel_exponential,
         KERNEL_LABEL_GAUSSIAN: _kernel_gaussian,
         KERNEL_LABEL_DENSITY: _kernel_density,
         # Use the user-provided beta args parameter if the user has provided
@@ -730,11 +728,6 @@ def execute(args):
     # kernel_creation_functions[KERNEL_LABEL_POWER].__name__ = (
     #     'functools_partial_decay_power')
 
-    # Since we have these keys defined in two places, I want to be super sure
-    # that the labels match.
-    assert sorted(kernel_creation_functions.keys()) == (
-        sorted(MODEL_SPEC['args']['decay_function']['options']))
-
     decay_function = args['decay_function']
     LOGGER.info(f'Using decay function {decay_function}')
 
@@ -970,17 +963,41 @@ def execute(args):
         kernel_path = os.path.join(
             intermediate_dir, f'kernel_{search_radius_m}{suffix}.tif')
         kernel_paths[search_radius_m] = kernel_path
+
+        if decay_function == KERNEL_LABEL_DICHOTOMY:
+            kernel_func = pygeoprocessing.kernels.dichotomous_kernel
+            kernel_kwargs = dict(
+                target_kernel_path=kernel_path,
+                max_distance=search_radius_in_pixels,
+                normalize=False)
+        elif decay_function == KERNEL_LABEL_EXPONENTIAL:
+            kernel_func = pygeoprocessing.kernels.exponential_decay_kernel
+            kernel_kwargs = dict(
+                target_kernel_path=kernel_path,
+                max_distance=math.ceil(expected_distance) * 2 + 1,
+                expected_distance=search_radius_in_pixels,
+                normalize=False)
+        elif decay_function in [KERNEL_LABEL_GAUSSIAN, KERNEL_LABEL_DENSITY]:
+            kernel_func = pygeoprocessing.kernels.create_distance_decay_kernel
+
+            def decay_func(dist_array):
+                return kernel_creation_functions[decay_function](
+                    dist_array, max_distance=search_radius_in_pixels)
+
+            kernel_kwargs = dict(
+                target_kernel_path=kernel_path,
+                distance_decay_function=decay_func,
+                max_distance=search_radius_in_pixels,
+                normalize=False)
+        else:
+            raise ValueError('Invalid kernel creation option selected')
+
         kernel_tasks[search_radius_m] = graph.add_task(
-            _create_kernel_raster,
-            kwargs={
-                'kernel_function': kernel_creation_functions[decay_function],
-                'expected_distance': search_radius_in_pixels,
-                'kernel_filepath': kernel_path,
-                'normalize': False},  # Model math calls for un-normalized
+            kernel_func,
+            kwargs=kernel_kwargs,
             task_name=(
                 f'Create {decay_function} kernel - {search_radius_m}m'),
-            target_path_list=[kernel_path]
-        )
+            target_path_list=[kernel_path])
 
     # Search radius mode 1: the same search radius applies to everything
     if args['search_radius_mode'] == RADIUS_OPT_UNIFORM:
@@ -1037,10 +1054,11 @@ def execute(args):
             intermediate_dir,
             f'urban_nature_population_ratio{suffix}.tif')
         urban_nature_population_ratio_task = graph.add_task(
-            _calculate_urban_nature_population_ratio,
-            args=(urban_nature_pixels_path,
-                  decayed_population_path,
-                  urban_nature_population_ratio_path),
+            func=pygeoprocessing.raster_map,
+            kwargs=dict(
+                op=_urban_nature_population_ratio,
+                rasters=[urban_nature_pixels_path, decayed_population_path],
+                target_path=urban_nature_population_ratio_path),
             task_name=(
                 '2SFCA: Calculate R_j urban nature/population ratio - '
                 f'{search_radius_m}'),
@@ -1127,10 +1145,11 @@ def execute(args):
                 intermediate_dir,
                 f'urban_nature_population_ratio_lucode_{lucode}{suffix}.tif')
             urban_nature_population_ratio_task = graph.add_task(
-                _calculate_urban_nature_population_ratio,
-                args=(urban_nature_pixels_path,
-                      decayed_population_paths[search_radius_m],
-                      urban_nature_population_ratio_path),
+                func=pygeoprocessing.raster_map,
+                kwargs=dict(
+                    op=_urban_nature_population_ratio,
+                    rasters=[urban_nature_pixels_path, decayed_population_paths[search_radius_m]],
+                    target_path=urban_nature_population_ratio_path),
                 task_name=(
                     '2SFCA: Calculate R_j urban nature/population ratio - '
                     f'{search_radius_m}'),
@@ -1161,13 +1180,11 @@ def execute(args):
                     urban_nature_population_ratio_task]))
 
         urban_nature_supply_percapita_task = graph.add_task(
-            ndr._sum_rasters,
-            kwargs={
-                'raster_path_list': partial_urban_nature_supply_percapita_paths,
-                'target_nodata': FLOAT32_NODATA,
-                'target_result_path':
-                    file_registry['urban_nature_supply_percapita'],
-            },
+            func=pygeoprocessing.raster_map,
+            kwargs=dict(
+                op=_sum_op,
+                rasters=partial_urban_nature_supply_percapita_paths,
+                target_path=file_registry['urban_nature_supply_percapita']),
             task_name='2SFCA - urban nature supply total',
             target_path_list=[file_registry['urban_nature_supply_percapita']],
             dependent_task_list=partial_urban_nature_supply_percapita_tasks
@@ -1238,22 +1255,22 @@ def execute(args):
             intermediate_dir,
             f'distance_weighted_population_all_groups{suffix}.tif')
         sum_of_decayed_population_task = graph.add_task(
-            ndr._sum_rasters,
-            kwargs={
-                'raster_path_list': decayed_population_in_group_paths,
-                'target_nodata': FLOAT32_NODATA,
-                'target_result_path': sum_of_decayed_population_path,
-            },
+            func=pygeoprocessing.raster_map,
+            kwargs=dict(
+                op=_sum_op,
+                rasters=decayed_population_in_group_paths,
+                target_path=sum_of_decayed_population_path),
             task_name='2SFCA - urban nature supply total',
             target_path_list=[sum_of_decayed_population_path],
             dependent_task_list=decayed_population_in_group_tasks
         )
 
         urban_nature_population_ratio_task = graph.add_task(
-            _calculate_urban_nature_population_ratio,
-            args=(urban_nature_pixels_path,
-                  sum_of_decayed_population_path,
-                  file_registry['urban_nature_population_ratio']),
+            func=pygeoprocessing.raster_map,
+            kwargs=dict(
+                op=_urban_nature_population_ratio,
+                rasters=[urban_nature_pixels_path, sum_of_decayed_population_path],
+                target_path=file_registry['urban_nature_population_ratio']),
             task_name=(
                 '2SFCA: Calculate R_j urban nature/population ratio - '
                 f'{search_radius_m}'),
@@ -1326,18 +1343,14 @@ def execute(args):
             urban_nature_balance_totalpop_by_group_paths[
                 pop_group] = urban_nature_balance_totalpop_by_group_path
             urban_nature_balance_totalpop_by_group_tasks.append(graph.add_task(
-                pygeoprocessing.raster_calculator,
-                kwargs={
-                    'base_raster_path_band_const_list': [
-                        (per_cap_urban_nature_balance_pop_group_path, 1),
-                        (proportional_pop_path, 1)
+                pygeoprocessing.raster_map,
+                kwargs=dict(
+                    op=_urban_nature_balance_totalpop_op,
+                    rasters=[
+                        per_cap_urban_nature_balance_pop_group_path,
+                        proportional_pop_path
                     ],
-                    'local_op': _urban_nature_balance_totalpop_op,
-                    'target_raster_path': (
-                        urban_nature_balance_totalpop_by_group_path),
-                    'datatype_target': gdal.GDT_Float32,
-                    'nodata_target': FLOAT32_NODATA
-                },
+                    target_path=urban_nature_balance_totalpop_by_group_path),
                 task_name='Calculate per-capita urban nature supply-demand',
                 target_path_list=[
                     urban_nature_balance_totalpop_by_group_path],
@@ -1414,14 +1427,11 @@ def execute(args):
             ])
 
         urban_nature_balance_totalpop_task = graph.add_task(
-            ndr._sum_rasters,
-            kwargs={
-                'raster_path_list':
-                    list(urban_nature_balance_totalpop_by_group_paths.values()),
-                'target_nodata': FLOAT32_NODATA,
-                'target_result_path':
-                    file_registry['urban_nature_balance_totalpop'],
-            },
+            func=pygeoprocessing.raster_map,
+            kwargs=dict(
+                op=_sum_op,
+                rasters=list(urban_nature_balance_totalpop_by_group_paths.values()),
+                target_path=file_registry['urban_nature_balance_totalpop']),
             task_name='2SFCA - urban nature - total population',
             target_path_list=[
                 file_registry['urban_nature_balance_totalpop']],
@@ -1478,18 +1488,15 @@ def execute(args):
 
         # This is "SUP_DEMi" from the user's guide
         urban_nature_balance_totalpop_task = graph.add_task(
-            pygeoprocessing.raster_calculator,
-            kwargs={
-                'base_raster_path_band_const_list': [
-                    (file_registry['urban_nature_balance_percapita'], 1),
-                    (file_registry['masked_population'], 1)
+            pygeoprocessing.raster_map,
+            kwargs=dict(
+                op=_urban_nature_balance_totalpop_op,
+                rasters=[
+                    file_registry['urban_nature_balance_percapita'],
+                    file_registry['masked_population']
                 ],
-                'local_op': _urban_nature_balance_totalpop_op,
-                'target_raster_path': (
-                    file_registry['urban_nature_balance_totalpop']),
-                'datatype_target': gdal.GDT_Float32,
-                'nodata_target': FLOAT32_NODATA
-            },
+                target_path=file_registry['urban_nature_balance_totalpop']
+            ),
             task_name='Calculate urban nature balance for the total population',
             target_path_list=[
                 file_registry['urban_nature_balance_totalpop']],
@@ -1567,6 +1574,10 @@ def execute(args):
     LOGGER.info('Finished Urban Nature Access Model')
 
 
+# Sum a list of arrays element-wise
+def _sum_op(*array_list): return numpy.sum(array_list, axis=0)
+
+
 def _geometries_overlap(vector_path):
     """Check if the geometries of the vector's first layer overlap.
 
@@ -1669,8 +1680,8 @@ def _weight_and_sum(*args):
         for source_array, weight_array, source_nodata, weight_nodata in zip(
                 pixel_arrays, weight_arrays, raster_nodata_list, weight_nodata_list):
             valid_pixels = (
-                ~utils.array_equals_nodata(source_array, source_nodata) &
-                ~utils.array_equals_nodata(weight_array, weight_nodata))
+                ~pygeoprocessing.array_equals_nodata(source_array, source_nodata) &
+                ~pygeoprocessing.array_equals_nodata(weight_array, weight_nodata))
             touched_pixels |= valid_pixels
             target_array[valid_pixels] += (
                 source_array[valid_pixels] * weight_array[valid_pixels])
@@ -1726,9 +1737,9 @@ def _reclassify_and_multiply(
         supply_block = supply_band.ReadAsArray(**block_info)
 
         valid_mask = (
-            ~utils.array_equals_nodata(
+            ~pygeoprocessing.array_equals_nodata(
                 pop_group_proportion_block, pop_group_nodata) &
-            ~utils.array_equals_nodata(supply_block, supply_nodata))
+            ~pygeoprocessing.array_equals_nodata(supply_block, supply_nodata))
         target_block = numpy.full(supply_block.shape, FLOAT32_NODATA,
                                   dtype=numpy.float32)
         target_block[valid_mask] = (
@@ -2145,112 +2156,69 @@ def _urban_nature_balance_totalpop_op(urban_nature_balance, population):
         A ``numpy.array`` of the area (in square meters) of urban nature
         supplied to each individual in each pixel.
     """
-    supply_demand = numpy.full(
-        urban_nature_balance.shape, FLOAT32_NODATA, dtype=numpy.float32)
-    valid_pixels = (
-        ~numpy.isclose(urban_nature_balance, FLOAT32_NODATA) &
-        ~numpy.isclose(population, FLOAT32_NODATA))
-    supply_demand[valid_pixels] = (
-        urban_nature_balance[valid_pixels] * population[valid_pixels])
-    return supply_demand
+    return urban_nature_balance * population
 
 
-def _calculate_urban_nature_population_ratio(
-        urban_nature_area_raster_path, convolved_population_raster_path,
-        target_ratio_raster_path):
+def _urban_nature_population_ratio(urban_nature_area, convolved_population):
     """Calculate the urban nature-population ratio R_j.
 
     Args:
-        urban_nature_area_raster_path (string): The path to a raster representing
-            the area of the pixel that represents urban nature.  Pixel values
-            will be ``0`` if there is no urban nature.
-        convolved_population_raster_path (string): The path to a raster
-            representing population counts that have been convolved over some
-            search kernel and perhaps weighted.
-        target_ratio_raster_path (string): The path to where the target
-            urban nature-population raster should be written.
+        urban_nature_area (numpy.array): A numpy array representing the area
+            of urban nature in the pixel.  Pixel values will be ``0`` if
+            there is no urban nature.  Pixel values may also match
+            ``urban_nature_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 weighted by a search kernel.
 
     Returns:
-        ``None``.
+        A numpy array with the ratio ``R_j`` representing the
+        urban nature-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-urban nature pixels will have a value of ``0`` in the
+              output matrix.
     """
-    urban_nature_nodata = pygeoprocessing.get_raster_info(
-        urban_nature_area_raster_path)['nodata'][0]
-    population_nodata = pygeoprocessing.get_raster_info(
-        convolved_population_raster_path)['nodata'][0]
-
-    def _urban_nature_population_ratio(urban_nature_area, convolved_population):
-        """Calculate the urban nature-population ratio R_j.
-
-        Args:
-            urban_nature_area (numpy.array): A numpy array representing the area
-                of urban nature in the pixel.  Pixel values will be ``0`` if
-                there is no urban nature.  Pixel values may also match
-                ``urban_nature_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 weighted by a search kernel.
-
-        Returns:
-            A numpy array with the ratio ``R_j`` representing the
-            urban nature-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-urban nature pixels will have a value of ``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(
-            urban_nature_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 urban nature pixels.
-        urban_nature_pixels = ~numpy.isclose(urban_nature_area, 0)
-        valid_pixels &= urban_nature_pixels
-        if population_nodata is not None:
-            valid_pixels &= ~utils.array_equals_nodata(
-                convolved_population, population_nodata)
-
-        if urban_nature_nodata is not None:
-            valid_pixels &= ~utils.array_equals_nodata(
-                urban_nature_area, urban_nature_nodata)
-
-        # The user's guide specifies that if the population in the search
-        # radius is numerically 0, the urban nature/population ratio should be
-        # set to the urban nature area.
-        # A consequence of this is that as the population approaches 0 from the
-        # positive side, the ratio will approach infinity.
-        # After checking with the science team, we decided that where the
-        # population is less than or equal to 1, the calculated
-        # urban nature/population ratio would be set to the available urban
-        # nature on that pixel.
-        population_close_to_zero = (convolved_population <= 1.0)
-        out_array[population_close_to_zero] = (
-            urban_nature_area[population_close_to_zero])
-        out_array[~urban_nature_pixels] = 0
-
-        valid_pixels_with_population = (
-            valid_pixels & (~population_close_to_zero))
-        out_array[valid_pixels_with_population] = (
-            urban_nature_area[valid_pixels_with_population] /
-            convolved_population[valid_pixels_with_population])
-
-        # eliminate pixel values < 0
-        out_array[valid_pixels & (out_array < 0)] = 0
-
-        return out_array
-
-    pygeoprocessing.raster_calculator(
-        [(urban_nature_area_raster_path, 1),
-         (convolved_population_raster_path, 1)],
-        _urban_nature_population_ratio, target_ratio_raster_path,
-        gdal.GDT_Float32, FLOAT32_NODATA)
+    # ASSUMPTION: population nodata value is not close to 0.
+    #  Shouldn't be if we're coming from convolution.
+    out_array = numpy.full(
+        urban_nature_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 urban nature pixels.
+    urban_nature_pixels = ~numpy.isclose(urban_nature_area, 0)
+    valid_pixels &= urban_nature_pixels
+
+    # The user's guide specifies that if the population in the search
+    # radius is numerically 0, the urban nature/population ratio should be
+    # set to the urban nature area.
+    # A consequence of this is that as the population approaches 0 from the
+    # positive side, the ratio will approach infinity.
+    # After checking with the science team, we decided that where the
+    # population is less than or equal to 1, the calculated
+    # urban nature/population ratio would be set to the available urban
+    # nature on that pixel.
+    population_close_to_zero = (convolved_population <= 1.0)
+    out_array[population_close_to_zero] = (
+        urban_nature_area[population_close_to_zero])
+    out_array[~urban_nature_pixels] = 0
+
+    valid_pixels_with_population = (
+        valid_pixels & (~population_close_to_zero))
+    out_array[valid_pixels_with_population] = (
+        urban_nature_area[valid_pixels_with_population] /
+        convolved_population[valid_pixels_with_population])
+
+    # eliminate pixel values < 0
+    out_array[valid_pixels & (out_array < 0)] = 0
+
+    return out_array
 
 
 def _convolve_and_set_lower_bound(
@@ -2284,7 +2252,7 @@ def _convolve_and_set_lower_bound(
     target_nodata = target_band.GetNoDataValue()
     for block_data, block in pygeoprocessing.iterblocks(
             (target_path, 1)):
-        valid_pixels = ~utils.array_equals_nodata(block, target_nodata)
+        valid_pixels = ~pygeoprocessing.array_equals_nodata(block, target_nodata)
         block[(block < 0) & valid_pixels] = 0
         target_band.WriteArray(
             block, xoff=block_data['xoff'], yoff=block_data['yoff'])
@@ -2368,7 +2336,6 @@ def _resample_population_raster(
     population_raster_info = pygeoprocessing.get_raster_info(
         source_population_raster_path)
     pixel_area = numpy.multiply(*population_raster_info['pixel_size'])
-    population_nodata = population_raster_info['nodata'][0]
 
     population_srs = osr.SpatialReference()
     population_srs.ImportFromWkt(population_raster_info['projection_wkt'])
@@ -2386,18 +2353,15 @@ def _convert_population_to_density(population):
 
         Returns:
             """
-        out_array = numpy.full(
-            population.shape, FLOAT32_NODATA, dtype=numpy.float32)
-        valid_mask = ~utils.array_equals_nodata(population, population_nodata)
-        out_array[valid_mask] = population[valid_mask] / population_pixel_area
-        return out_array
+        return population / population_pixel_area
 
     # Step 1: convert the population raster to population density per sq. km
     density_raster_path = os.path.join(tmp_working_dir, 'pop_density.tif')
-    pygeoprocessing.raster_calculator(
-        [(source_population_raster_path, 1)],
-        _convert_population_to_density,
-        density_raster_path, gdal.GDT_Float32, FLOAT32_NODATA)
+    pygeoprocessing.raster_map(
+        rasters=[source_population_raster_path],
+        op=_convert_population_to_density,
+        target_path=density_raster_path,
+        target_dtype=numpy.float32)
 
     # Step 2: align to the LULC
     warped_density_path = os.path.join(tmp_working_dir, 'warped_density.tif')
@@ -2432,60 +2396,16 @@ def _convert_density_to_population(density):
         # between multiple pixels.  So it's preserving an unrealistic degree of
         # precision, but that's probably OK because pixels are imprecise
         # measures anyways.
-        out_array = numpy.full(
-            density.shape, FLOAT32_NODATA, dtype=numpy.float32)
+        return density * target_pixel_area
 
-        # We already know that the nodata value is FLOAT32_NODATA
-        valid_mask = ~numpy.isclose(density, FLOAT32_NODATA)
-        out_array[valid_mask] = density[valid_mask] * target_pixel_area
-        return out_array
-
-    pygeoprocessing.raster_calculator(
-        [(warped_density_path, 1)],
-        _convert_density_to_population,
-        target_population_raster_path, gdal.GDT_Float32, FLOAT32_NODATA)
+    pygeoprocessing.raster_map(
+        op=_convert_density_to_population,
+        rasters=[warped_density_path],
+        target_path=target_population_raster_path)
 
     shutil.rmtree(tmp_working_dir, ignore_errors=True)
 
 
-def _kernel_dichotomy(distance, max_distance):
-    """Create a dichotomous kernel.
-
-    All pixels within ``max_distance`` have a value of 1.
-
-    Args:
-        distance (numpy.array): An array of euclidean distances (in pixels)
-            from the center of the kernel.
-        max_distance (float): The maximum distance of the kernel.  Pixels that
-            are more than this number of pixels will have a value of 0.
-
-    Returns:
-        ``numpy.array`` with dtype of numpy.float32 and same shape as
-        ``distance.
-    """
-    return (distance <= max_distance).astype(numpy.float32)
-
-
-def _kernel_exponential(distance, max_distance):
-    """Create an exponential-decay kernel.
-
-    Args:
-        distance (numpy.array): An array of euclidean distances (in pixels)
-            from the center of the kernel.
-        max_distance (float): The maximum distance of the kernel.  Pixels that
-            are more than this number of pixels will have a value of 0.
-
-    Returns:
-        ``numpy.array`` with dtype of numpy.float32 and same shape as
-        ``distance.
-    """
-    kernel = numpy.zeros(distance.shape, dtype=numpy.float32)
-    pixels_in_radius = (distance <= max_distance)
-    kernel[pixels_in_radius] = numpy.exp(
-        -distance[pixels_in_radius] / max_distance)
-    return kernel
-
-
 def _kernel_power(distance, max_distance, beta):
     """Create a power kernel with user-defined beta.
 
@@ -2551,95 +2471,6 @@ def _kernel_density(distance, max_distance):
     return kernel
 
 
-def _create_kernel_raster(
-        kernel_function, expected_distance, kernel_filepath, normalize=False):
-    """Create a raster distance-weighted decay kernel from a function.
-
-    Args:
-        kernel_function (callable): The kernel function to use.
-        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
-    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])
-
-        kernel = kernel_function(distance=pixel_dist_from_center,
-                                 max_distance=expected_distance)
-        if normalize:
-            running_sum += kernel.sum()
-
-        kernel_band.WriteArray(
-            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 _create_valid_pixels_nodata_mask(raster_list, target_mask_path):
     """Create a valid pixels mask across a stack of aligned rasters.
 
@@ -2661,7 +2492,7 @@ def _create_valid_pixels_nodata_mask(raster_list, target_mask_path):
     def _create_mask(*raster_arrays):
         valid_pixels_mask = numpy.ones(raster_arrays[0].shape, dtype=bool)
         for nodata, array in zip(nodatas, raster_arrays):
-            valid_pixels_mask &= ~utils.array_equals_nodata(array, nodata)
+            valid_pixels_mask &= ~pygeoprocessing.array_equals_nodata(array, nodata)
 
         return valid_pixels_mask
 
diff --git a/src/natcap/invest/utils.py b/src/natcap/invest/utils.py
index e70893395e..d845444e0a 100644
--- a/src/natcap/invest/utils.py
+++ b/src/natcap/invest/utils.py
@@ -326,193 +326,6 @@ def make_suffix_string(args, suffix_key):
     return file_suffix
 
 
-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.
-
-    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.
-        normalize=True (bool): Whether to divide the kernel values by the sum
-            of all values in the kernel.
-
-    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)
-
-    cols_per_block, rows_per_block = kernel_band.GetBlockSize()
-
-    n_cols = kernel_dataset.RasterXSize
-    n_rows = kernel_dataset.RasterYSize
-
-    n_col_blocks = int(math.ceil(n_cols / float(cols_per_block)))
-    n_row_blocks = int(math.ceil(n_rows / float(rows_per_block)))
-
-    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
-
-        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
-
-            # 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)
-
-            kernel_index_distances = numpy.hypot(
-                row_indices, col_indices)
-            kernel = numpy.where(
-                kernel_index_distances > max_distance, 0.0,
-                numpy.exp(-kernel_index_distances / expected_distance))
-            integration += numpy.sum(kernel)
-
-            kernel_band.WriteArray(kernel, xoff=col_offset,
-                                   yoff=row_offset)
-
-    # Need to flush the kernel's cache to disk before opening up a new Dataset
-    # object in interblocks()
-    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 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_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.
 
@@ -1003,30 +816,6 @@ def reclassify_raster(
         raise ValueError(error_message)
 
 
-def array_equals_nodata(array, nodata):
-    """Check for the presence of ``nodata`` values in ``array``.
-
-    The comparison supports ``numpy.nan`` and unset (``None``) nodata values.
-
-    Args:
-        array (numpy array): the array to mask for nodata values.
-        nodata (number): the nodata value to check for. Supports ``numpy.nan``.
-
-    Returns:
-        A boolean numpy array with values of 1 where ``array`` is equal to
-        ``nodata`` and 0 otherwise.
-    """
-    # If nodata is undefined, nothing matches nodata.
-    if nodata is None:
-        return numpy.zeros(array.shape, dtype=bool)
-
-    # comparing an integer array against numpy.nan works correctly and is
-    # faster than using numpy.isclose().
-    if numpy.issubdtype(array.dtype, numpy.integer):
-        return array == nodata
-    return numpy.isclose(array, nodata, equal_nan=True)
-
-
 def matches_format_string(test_string, format_string):
     """Assert that a given string matches a given format string.
 
diff --git a/src/natcap/invest/wave_energy.py b/src/natcap/invest/wave_energy.py
index 2cf4acc1a2..dc0a9a7ad6 100644
--- a/src/natcap/invest/wave_energy.py
+++ b/src/natcap/invest/wave_energy.py
@@ -1719,7 +1719,7 @@ def _create_percentile_rasters(base_raster_path, target_raster_path,
 
     def _mask_below_start_value(array):
         valid_mask = (
-            ~utils.array_equals_nodata(array, base_nodata) &
+            ~pygeoprocessing.array_equals_nodata(array, base_nodata) &
             (array >= float(start_value)))
         result = numpy.empty_like(array)
         result[:] = base_nodata
@@ -1760,19 +1760,13 @@ def _mask_below_start_value(array):
     value_ranges.append('Greater than %s' % rounded_percentiles[-1])
     LOGGER.debug('Range_values : %s', value_ranges)
 
-    def raster_percentile(band):
-        """Group the band pixels together based on _PERCENTILES, starting from 1.
-        """
-        valid_data_mask = ~utils.array_equals_nodata(band, base_nodata)
-        band[valid_data_mask] = numpy.searchsorted(
-            percentile_values, band[valid_data_mask]) + 1
-        band[~valid_data_mask] = target_nodata
-        return band
-
     # Classify the pixels of raster_dataset into groups and write to output
-    pygeoprocessing.raster_calculator([(base_raster_path, 1)],
-                                      raster_percentile, target_raster_path,
-                                      gdal.GDT_Byte, target_nodata)
+    pygeoprocessing.raster_map(
+        op=lambda band: numpy.searchsorted(percentile_values, band) + 1,
+        rasters=[base_raster_path],
+        target_path=target_raster_path,
+        target_dtype=numpy.uint8,
+        target_nodata=target_nodata)
 
     # Create percentile groups of how percentile ranges are classified
     percentile_groups = numpy.arange(1, len(percentile_values) + 2)
diff --git a/src/natcap/invest/wind_energy.py b/src/natcap/invest/wind_energy.py
index cc5b7d5c22..bdb33ede60 100644
--- a/src/natcap/invest/wind_energy.py
+++ b/src/natcap/invest/wind_energy.py
@@ -1445,7 +1445,7 @@ def _calculate_npv_levelized_rasters(
             pygeoprocessing.iterblocks((base_dist_raster_path, 1))):
 
         target_arr_shape = harvest_block_data.shape
-        target_nodata_mask = utils.array_equals_nodata(
+        target_nodata_mask = pygeoprocessing.array_equals_nodata(
             harvest_block_data, _TARGET_NODATA)
 
         # Total cable distance converted to kilometers
@@ -1624,7 +1624,7 @@ def _depth_op(bath, min_depth, max_depth):
     out_array = numpy.full(
         bath.shape, _TARGET_NODATA, dtype=numpy.float32)
     valid_pixels_mask = ((bath >= max_depth) & (bath <= min_depth) &
-                         ~utils.array_equals_nodata(bath, _TARGET_NODATA))
+                         ~pygeoprocessing.array_equals_nodata(bath, _TARGET_NODATA))
     out_array[
         valid_pixels_mask] = bath[valid_pixels_mask]
     return out_array
@@ -1644,7 +1644,7 @@ def _add_avg_dist_op(tmp_dist, avg_grid_distance):
     """
     out_array = numpy.full(
         tmp_dist.shape, _TARGET_NODATA, dtype=numpy.float32)
-    valid_pixels_mask = ~utils.array_equals_nodata(tmp_dist, _TARGET_NODATA)
+    valid_pixels_mask = ~pygeoprocessing.array_equals_nodata(tmp_dist, _TARGET_NODATA)
     out_array[valid_pixels_mask] = (
         tmp_dist[valid_pixels_mask] + avg_grid_distance)
     return out_array
@@ -1711,7 +1711,7 @@ def _mask_out_depth_dist(*rasters):
     out_array = numpy.full(rasters[0].shape, _TARGET_NODATA, dtype=numpy.float32)
     nodata_mask = numpy.full(rasters[0].shape, False, dtype=bool)
     for array in rasters:
-        nodata_mask = nodata_mask | utils.array_equals_nodata(
+        nodata_mask = nodata_mask | pygeoprocessing.array_equals_nodata(
                 array, _TARGET_NODATA)
     out_array[~nodata_mask] = rasters[0][~nodata_mask]
     return out_array
@@ -1731,7 +1731,7 @@ def _calculate_carbon_op(harvested_arr, carbon_coef):
     """
     out_array = numpy.full(
         harvested_arr.shape, _TARGET_NODATA, dtype=numpy.float32)
-    valid_pixels_mask = ~utils.array_equals_nodata(
+    valid_pixels_mask = ~pygeoprocessing.array_equals_nodata(
         harvested_arr, _TARGET_NODATA)
 
     # The energy value converted from MWhr/yr (Mega Watt hours as output
@@ -1877,7 +1877,7 @@ def _dist_mask_op(dist_arr):
         """Mask distance values by min/max values."""
         out_array = numpy.full(dist_arr.shape, out_nodata, dtype=numpy.float32)
         valid_pixels_mask = (
-            ~utils.array_equals_nodata(dist_arr, raster_nodata) &
+            ~pygeoprocessing.array_equals_nodata(dist_arr, raster_nodata) &
             (dist_arr >= min_dist) & (dist_arr <= max_dist))
         out_array[
             valid_pixels_mask] = dist_arr[valid_pixels_mask]
diff --git a/tests/test_coastal_vulnerability.py b/tests/test_coastal_vulnerability.py
index ac92b240fc..9f3f7b51d2 100644
--- a/tests/test_coastal_vulnerability.py
+++ b/tests/test_coastal_vulnerability.py
@@ -747,7 +747,7 @@ def test_aggregate_raster_edge_cases(self):
             ogr_geom_type=ogr.wkbPoint)
 
         coastal_vulnerability._aggregate_raster_values_in_radius(
-            simple_points_path, raster_path, search_radius,
+            simple_points_path, raster_path, search_radius, self.workspace_dir,
             target_pickle_path, _mean_op)
 
         with open(target_pickle_path, 'rb') as file:
@@ -800,7 +800,7 @@ def test_nodata_raster_aggregation(self):
 
         coastal_vulnerability._aggregate_raster_values_in_radius(
             simple_points_path, raster_path, sample_distance,
-            target_pickle_path, _mean_op)
+            self.workspace_dir, target_pickle_path, _mean_op)
 
         with open(target_pickle_path, 'rb') as file:
             actual_values = pickle.load(file)
diff --git a/tests/test_habitat_quality.py b/tests/test_habitat_quality.py
index fc9ceab06d..560793d2c6 100644
--- a/tests/test_habitat_quality.py
+++ b/tests/test_habitat_quality.py
@@ -1117,7 +1117,7 @@ def test_habitat_quality_missing_lucodes_in_table(self):
 
         scenarios = ['_bas_', '_cur_', '_fut_']
         for lulc_val, scenario in enumerate(scenarios, start=1):
-            lulc_array = numpy.ones((100, 100), dtype=numpy.int8)
+            lulc_array = numpy.ones((100, 100), dtype=numpy.uint8)
             lulc_array[50:, :] = lulc_val
             path = os.path.join(
                 args['workspace_dir'], 'lc_samp' + scenario + 'b.tif')
diff --git a/tests/test_pollination.py b/tests/test_pollination.py
index fdb2a4edd8..1897e9bff4 100644
--- a/tests/test_pollination.py
+++ b/tests/test_pollination.py
@@ -198,7 +198,7 @@ def test_pollination_no_farm_regression(self):
             result_sum += numpy.sum(data_block)
         # the number below is just what the sum rounded to two decimal places
         # when I manually inspected a run that appeared to be correct.
-        self.assertAlmostEqual(result_sum, 58.669518, places=2)
+        self.assertAlmostEqual(result_sum, 58.407155, places=2)
 
     def test_pollination_constant_abundance(self):
         """Pollination: regression testing when abundance is all 1."""
@@ -223,7 +223,7 @@ def test_pollination_constant_abundance(self):
             result_sum += numpy.sum(data_block)
         # the number below is just what the sum rounded to two decimal places
         # when I manually inspected a run that appeared to be correct.
-        self.assertAlmostEqual(result_sum, 68.44777, places=2)
+        self.assertAlmostEqual(result_sum, 68.14167, places=2)
 
     def test_pollination_bad_guild_headers(self):
         """Pollination: testing that model detects bad guild headers."""
diff --git a/tests/test_sdr.py b/tests/test_sdr.py
index 436072fec1..3c324934b0 100644
--- a/tests/test_sdr.py
+++ b/tests/test_sdr.py
@@ -427,7 +427,8 @@ def test_ls_factor(self):
                                  target_ls_factor_path)
 
         ls = pygeoprocessing.raster_to_numpy_array(target_ls_factor_path)
+        nodata = float(numpy.finfo(numpy.float32).max)
         expected_ls = numpy.array(
-            [[0.253996, 0.657229, 1.345856, 1.776729, 49.802994, -1]],
+            [[0.253996, 0.657229, 1.345856, 1.776729, 49.802994, nodata]],
             dtype=numpy.float32)
         numpy.testing.assert_allclose(ls, expected_ls, rtol=1e-6)
diff --git a/tests/test_seasonal_water_yield_regression.py b/tests/test_seasonal_water_yield_regression.py
index e7cb6ae4cd..899fc1446b 100644
--- a/tests/test_seasonal_water_yield_regression.py
+++ b/tests/test_seasonal_water_yield_regression.py
@@ -1127,42 +1127,6 @@ def test_monthly_quickflow_negative_values_set_to_zero(self):
             pygeoprocessing.raster_to_numpy_array(output_path),
             expected_quickflow_array, atol=1e-5)
 
-    def test_calculate_annual_qfi_different_nodata_areas(self):
-        """Test with qf rasters with different areas of nodata."""
-        from natcap.invest.seasonal_water_yield import seasonal_water_yield
-        qf_array_1 = numpy.array([
-            [10, 10],
-            [10, -1]], dtype=numpy.float32)
-        qf_array_2 = numpy.array([
-            [10, 10],
-            [-1, 10]], dtype=numpy.float32)
-        qf_array_3 = numpy.array([
-            [-1, 10],
-            [10, 10]], dtype=numpy.float32)
-
-        qf_1_path = os.path.join(self.workspace_dir, 'qf_1.tif')
-        qf_2_path = os.path.join(self.workspace_dir, 'qf_2.tif')
-        qf_3_path = os.path.join(self.workspace_dir, 'qf_3.tif')
-
-        srs = osr.SpatialReference()
-        srs.ImportFromEPSG(26910)  # UTM Zone 10N
-        project_wkt = srs.ExportToWkt()
-        output_path = os.path.join(self.workspace_dir, 'qf_sum.tif')
-
-        # write all the test arrays to raster files
-        for array, path in [(qf_array_1, qf_1_path),
-                            (qf_array_2, qf_2_path),
-                            (qf_array_3, qf_3_path)]:
-            pygeoprocessing.numpy_array_to_raster(
-                array, -1, (1, -1), (1180000, 690000), project_wkt, path)
-        seasonal_water_yield._calculate_annual_qfi(
-            [qf_1_path, qf_2_path, qf_3_path], output_path)
-        numpy.testing.assert_allclose(
-            pygeoprocessing.raster_to_numpy_array(output_path),
-            numpy.array([
-                [-1, 30],
-                [-1, -1]], dtype=numpy.float32))
-
     def test_local_recharge_undefined_nodata(self):
         """Test `calculate_local_recharge` with undefined nodata values"""
         from natcap.invest.seasonal_water_yield import \
diff --git a/tests/test_stormwater.py b/tests/test_stormwater.py
index ba5b2b9737..dbf6f7e442 100644
--- a/tests/test_stormwater.py
+++ b/tests/test_stormwater.py
@@ -629,50 +629,6 @@ def test_retention_value_op(self):
             replacement_cost)
         numpy.testing.assert_allclose(actual, expected)
 
-    def test_adjust_op(self):
-        """Stormwater: test adjust_op function."""
-        from natcap.invest import stormwater
-
-        ratio_array = numpy.array([
-            [0,   0.0001, stormwater.FLOAT_NODATA],
-            [0.5, 0.9,    1]], dtype=numpy.float32)
-        # these are obviously not averages from the above array but
-        # it doesn't matter for this test
-        avg_ratio_array = numpy.array([
-            [0.5, 0.5, 0.5],
-            [0.5, stormwater.FLOAT_NODATA, 0.5]], dtype=numpy.float32)
-        near_connected_lulc_array = numpy.array([
-            [0, 0, 1],
-            [stormwater.UINT8_NODATA, 0, 1]], dtype=numpy.uint8)
-        near_road_centerline_array = numpy.array([
-            [1, 1, 1],
-            [0, 0, 0]], dtype=numpy.uint8)
-
-        out = stormwater.adjust_op(
-            ratio_array,
-            avg_ratio_array,
-            near_connected_lulc_array,
-            near_road_centerline_array)
-        for y in range(ratio_array.shape[0]):
-            for x in range(ratio_array.shape[1]):
-                if (ratio_array[y, x] == stormwater.FLOAT_NODATA or
-                    avg_ratio_array[y, x] == stormwater.FLOAT_NODATA or
-                    near_connected_lulc_array[y, x] == stormwater.UINT8_NODATA or
-                        near_road_centerline_array[y, x] == stormwater.UINT8_NODATA):
-                    numpy.testing.assert_allclose(
-                        out[y, x], stormwater.FLOAT_NODATA)
-                else:
-                    # equation 2-4: Radj_ij = R_ij + (1 - R_ij) * C_ij
-                    adjust_factor = (
-                        0 if (
-                            near_connected_lulc_array[y, x] or
-                            near_road_centerline_array[y, x]
-                        ) else avg_ratio_array[y, x])
-                    adjusted = (ratio_array[y, x] +
-                                (1 - ratio_array[y, x]) * adjust_factor)
-                    numpy.testing.assert_allclose(out[y, x], adjusted,
-                                                  rtol=1e-6)
-
     def test_is_near(self):
         """Stormwater: test is_near function."""
         from natcap.invest import stormwater
@@ -710,36 +666,6 @@ def test_is_near(self):
             actual = pygeoprocessing.raster_to_numpy_array(out_path)
             numpy.testing.assert_equal(expected, actual)
 
-    def test_make_search_kernel(self):
-        """Stormwater: test make_search_kernel function."""
-        from natcap.invest import stormwater
-
-        array = numpy.zeros((10, 10))
-        path = os.path.join(self.workspace_dir, 'make_search_kernel.tif')
-        to_raster(array, path, pixel_size=(10, -10))
-
-        expected_5 = numpy.array([[1]], dtype=numpy.uint8)
-        actual_5 = stormwater.make_search_kernel(path, 5)
-        numpy.testing.assert_equal(expected_5, actual_5)
-
-        expected_9 = numpy.array([[1]], dtype=numpy.uint8)
-        actual_9 = stormwater.make_search_kernel(path, 9)
-        numpy.testing.assert_equal(expected_9, actual_9)
-
-        expected_10 = numpy.array([
-            [0, 1, 0],
-            [1, 1, 1],
-            [0, 1, 0]], dtype=numpy.uint8)
-        actual_10 = stormwater.make_search_kernel(path, 10)
-        numpy.testing.assert_equal(expected_10, actual_10)
-
-        expected_15 = numpy.array([
-            [1, 1, 1],
-            [1, 1, 1],
-            [1, 1, 1]], dtype=numpy.uint8)
-        actual_15 = stormwater.make_search_kernel(path, 15)
-        numpy.testing.assert_equal(expected_15, actual_15)
-
     def test_raster_average(self):
         """Stormwater: test raster_average function."""
         from natcap.invest import stormwater
diff --git a/tests/test_ucm.py b/tests/test_ucm.py
index a900bbf312..1696cd4e58 100644
--- a/tests/test_ucm.py
+++ b/tests/test_ucm.py
@@ -69,7 +69,7 @@ def test_ucm_regression_factors(self):
             'avg_cc': 0.222150472947109,
             'avg_tmp_v': 37.325275675470998,
             'avg_tmp_an': 2.325275675470998,
-            'avd_eng_cn': 3520213.280928277,
+            'avd_eng_cn': 3520217.313878,
             'avg_wbgt_v': 32.60417266705069,
             'avg_ltls_v': 75.000000000000000,
             'avg_hvls_v': 75.000000000000000,
@@ -90,7 +90,7 @@ def test_ucm_regression_factors(self):
 
         # Assert that the decimal value of the energy savings value is what we
         # expect.
-        expected_energy_sav = 3564034.496484185
+        expected_energy_sav = 3564038.678764
 
         energy_sav = 0.0
         n_nonetype = 0
@@ -206,9 +206,9 @@ def test_ucm_regression_intensity(self):
             'avg_cc': 0.428302583240327,
             'avg_tmp_v': 36.60869797039769,
             'avg_tmp_an': 1.608697970397692,
-            'avd_eng_cn': 7240015.1958200345,
+            'avd_eng_cn': 7239992.744486,
             'avg_wbgt_v': 31.91108630952381,
-            'avg_ltls_v': 28.744239631336406,
+            'avg_ltls_v': 28.73463901689708,
             'avg_hvls_v': 75.000000000000000,
         }
         try:
@@ -412,21 +412,6 @@ def test_bad_args(self):
                 header='column', header_name='green_area'))]
         self.assertEqual(result, expected)
 
-    def test_flat_disk_kernel(self):
-        """UCM: test flat disk kernel."""
-        import natcap.invest.urban_cooling_model
-
-        kernel_filepath = os.path.join(self.workspace_dir, 'kernel.tif')
-        natcap.invest.urban_cooling_model.flat_disk_kernel(
-            1000, kernel_filepath)
-
-        kernel_raster = gdal.OpenEx(kernel_filepath, gdal.OF_RASTER)
-        kernel_band = kernel_raster.GetRasterBand(1)
-        self.assertAlmostEqual(
-            numpy.sum(kernel_band.ReadAsArray())/1000,
-            numpy.ceil(1000**2*numpy.pi/1000),
-            places=0)
-
     def test_do_energy_valuation_option(self):
         """UCM: test separate valuation options."""
         import natcap.invest.urban_cooling_model
diff --git a/tests/test_urban_nature_access.py b/tests/test_urban_nature_access.py
index 3c4c3f6e3c..2cd43c1b29 100644
--- a/tests/test_urban_nature_access.py
+++ b/tests/test_urban_nature_access.py
@@ -169,132 +169,16 @@ def test_resample_population_raster(self):
                     population_array.sum(), resampled_population_array.sum(),
                     rtol=1e-3)
 
-    def test_dichotomous_decay_simple(self):
-        """UNA: Test dichotomous decay kernel on a simple case."""
+    def test_density_kernel(self):
+        """UNA: Test the density kernel."""
         from natcap.invest import urban_nature_access
 
-        expected_distance = 5
-        kernel_filepath = os.path.join(self.workspace_dir, 'kernel.tif')
-
-        urban_nature_access._create_kernel_raster(
-            urban_nature_access._kernel_dichotomy, 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_normalized(self):
-        """UNA: Test normalized dichotomous kernel."""
-        from natcap.invest import urban_nature_access
-
-        expected_distance = 5
-        kernel_filepath = os.path.join(self.workspace_dir, 'kernel.tif')
-
-        urban_nature_access._create_kernel_raster(
-            urban_nature_access._kernel_dichotomy,
-            expected_distance, kernel_filepath, normalize=True)
-
-        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.float32)
-        expected_array /= numpy.sum(expected_array)
-
-        extracted_kernel_array = pygeoprocessing.raster_to_numpy_array(
-            kernel_filepath)
-        numpy.testing.assert_allclose(
-            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._create_kernel_raster(
-            urban_nature_access._kernel_dichotomy,
-            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_simple(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._create_kernel_raster(
-            urban_nature_access._kernel_density,
-            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_density_decay_normalized(self):
-        """UNA: Test normalized 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._create_kernel_raster(
-            urban_nature_access._kernel_density,
-            expected_distance, kernel_filepath, normalize=True)
-
-        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(1, kernel_array.sum())
-        self.assertAlmostEqual(1.5915502e-05, kernel_array.max())
-        self.assertEqual(0, kernel_array.min())
+        max_distance = 3
+        distance = numpy.array([0, 1, 2, 3, 4])
+        kernel = urban_nature_access._kernel_density(distance, max_distance)
+        # These regression values are calculated by hand
+        expected_array = numpy.array([.75, 2/3, 5/12, 0, 0])
+        numpy.testing.assert_allclose(expected_array, kernel)
 
     def test_power_kernel(self):
         """UNA: Test the power kernel."""
@@ -310,20 +194,6 @@ def test_power_kernel(self):
         numpy.testing.assert_allclose(
             expected_array, kernel)
 
-    def test_exponential_kernel(self):
-        """UNA: Test the exponential decay kernel."""
-        from natcap.invest import urban_nature_access
-
-        max_distance = 3
-        distance = numpy.array([0, 1, 2, 3, 4])
-        kernel = urban_nature_access._kernel_exponential(
-            distance, max_distance)
-        # Regression values are calculated by hand
-        expected_array = numpy.array(
-            [1, 0.71653134, 0.5134171, 0.36787945, 0])
-        numpy.testing.assert_allclose(
-            expected_array, kernel)
-
     def test_gaussian_kernel(self):
         """UNA: Test the gaussian decay kernel."""
         from natcap.invest import urban_nature_access
@@ -361,14 +231,6 @@ def test_urban_nature_balance(self):
         numpy.testing.assert_allclose(
             urban_nature_budget, expected_urban_nature_budget)
 
-        supply_demand = urban_nature_access._urban_nature_balance_totalpop_op(
-            urban_nature_budget, population)
-        expected_supply_demand = numpy.array([
-            [nodata, 100 * 50.5],
-            [25 * 40.75, nodata]], dtype=numpy.float32)
-        numpy.testing.assert_allclose(
-            supply_demand, expected_supply_demand)
-
     def test_reclassify_and_multpliy(self):
         """UNA: test reclassification/multiplication function."""
         from natcap.invest import urban_nature_access
@@ -484,7 +346,7 @@ def test_core_model(self):
         accessible_urban_nature_array = pygeoprocessing.raster_to_numpy_array(
             os.path.join(args['workspace_dir'], 'output',
                          'accessible_urban_nature_suffix.tif'))
-        valid_mask = ~utils.array_equals_nodata(
+        valid_mask = ~pygeoprocessing.array_equals_nodata(
             accessible_urban_nature_array, urban_nature_access.FLOAT32_NODATA)
         valid_pixels = accessible_urban_nature_array[valid_mask]
         self.assertAlmostEqual(numpy.sum(valid_pixels), 6221004.41259766)
@@ -603,7 +465,7 @@ def test_split_population(self):
         def _read_and_sum_raster(path):
             array = pygeoprocessing.raster_to_numpy_array(path)
             nodata = pygeoprocessing.get_raster_info(path)['nodata'][0]
-            return numpy.sum(array[~utils.array_equals_nodata(array, nodata)])
+            return numpy.sum(array[~pygeoprocessing.array_equals_nodata(array, nodata)])
 
         intermediate_dir = os.path.join(args['workspace_dir'], 'intermediate')
         for (supply_type, supply_field), fieldname in itertools.product(
@@ -652,7 +514,7 @@ def _assert_urban_nature(self, path, sum_value, min_value, max_value):
 
         accessible_urban_nature_array = (
             pygeoprocessing.raster_to_numpy_array(path))
-        valid_mask = ~utils.array_equals_nodata(
+        valid_mask = ~pygeoprocessing.array_equals_nodata(
             accessible_urban_nature_array,
             urban_nature_access.FLOAT32_NODATA)
         valid_pixels = accessible_urban_nature_array[valid_mask]
@@ -713,8 +575,8 @@ def test_radii_by_pop_group(self):
             'Povr_adm': 0,
             'Povr_adm_female': 607.13671875,
             'Povr_adm_male': 477.0360107421875,
-            'SUP_DEMadm_cap': -17.90779987933412,
-            'SUP_DEMadm_cap_female': -17.907799675104435,
+            'SUP_DEMadm_cap': -17.907799109781322,
+            'SUP_DEMadm_cap_female': -17.90779830090304,
             'SUP_DEMadm_cap_male': -17.907800139262825,
         }
         self.assertEqual(
@@ -818,7 +680,7 @@ def test_modes_same_radii_same_results(self):
                 os.path.join(args['workspace_dir'], 'output',
                              f'urban_nature_demand_{suffix}.tif'))
             nodata = urban_nature_access.FLOAT32_NODATA
-            valid_pixels = ~utils.array_equals_nodata(population, nodata)
+            valid_pixels = ~pygeoprocessing.array_equals_nodata(population, nodata)
             numpy.testing.assert_allclose(
                 (population[valid_pixels].sum() *
                     float(args['urban_nature_demand'])),
diff --git a/tests/test_utils.py b/tests/test_utils.py
index 95d7ae5822..74f69821b4 100644
--- a/tests/test_utils.py
+++ b/tests/test_utils.py
@@ -180,93 +180,6 @@ def _norm_dict(path_dict):
         return result_dict
 
 
-class ExponentialDecayUtilsTests(unittest.TestCase):
-    """Tests for natcap.invest.utils.exponential_decay_kernel_raster."""
-
-    _REGRESSION_PATH = os.path.join(
-        os.path.dirname(__file__), '..', 'data', 'invest-test-data',
-        'exp_decay_kernel')
-
-    def setUp(self):
-        """Setup workspace."""
-        self.workspace_dir = tempfile.mkdtemp()
-
-    def tearDown(self):
-        """Delete workspace."""
-        shutil.rmtree(self.workspace_dir)
-
-    def test_exp_decay_kernel_raster(self):
-        """Utils: test exponential_decay_kernel_raster."""
-        from natcap.invest import utils
-        expected_distance = 100  # 10 pixels
-        kernel_filepath = os.path.join(self.workspace_dir, 'kernel_100.tif')
-        utils.exponential_decay_kernel_raster(
-            expected_distance, kernel_filepath)
-
-        model_array = pygeoprocessing.raster_to_numpy_array(
-            kernel_filepath)
-        reg_array = pygeoprocessing.raster_to_numpy_array(
-            os.path.join(
-                ExponentialDecayUtilsTests._REGRESSION_PATH,
-                'kernel_100.tif'))
-        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
-            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."""
 
@@ -1754,44 +1667,3 @@ def test_exception_raised_with_details(self):
             " the LULC raster but not the table are: [3].")
         self.assertTrue(
             expected_message in str(context.exception), str(context.exception))
-
-
-class ArrayEqualsNodataTests(unittest.TestCase):
-    """Tests for natcap.invest.utils.array_equals_nodata."""
-
-    def test_integer_array(self):
-        """Utils: test integer array is returned as expected."""
-        from natcap.invest import utils
-
-        nodata_values = [9, 9.0]
-
-        int_array = numpy.array(
-            [[4, 2, 9], [1, 9, 3], [9, 6, 1]], dtype=numpy.int16)
-
-        expected_array = numpy.array([[0, 0, 1], [0, 1, 0], [1, 0, 0]])
-
-        for nodata in nodata_values:
-            result_array = utils.array_equals_nodata(int_array, nodata)
-            numpy.testing.assert_array_equal(result_array, expected_array)
-
-    def test_nan_nodata_array(self):
-        """Utils: test array with numpy.nan nodata values."""
-        from natcap.invest import utils
-
-        array = numpy.array(
-            [[4, 2, numpy.nan], [1, numpy.nan, 3], [numpy.nan, 6, 1]])
-
-        expected_array = numpy.array([[0, 0, 1], [0, 1, 0], [1, 0, 0]])
-
-        result_array = utils.array_equals_nodata(array, numpy.nan)
-        numpy.testing.assert_array_equal(result_array, expected_array)
-
-    def test_none_nodata(self):
-        """Utils: if nodata is undefined (None), everything is valid."""
-        from natcap.invest import utils
-
-        array = numpy.array(
-            [[4, 2, numpy.nan], [1, numpy.nan, 3], [numpy.nan, 6, 1]])
-        result_array = utils.array_equals_nodata(array, None)
-        numpy.testing.assert_array_equal(
-            result_array, numpy.zeros(array.shape, dtype=bool))