Merge pull request natcap#1437 from emlys/task/update-pygeoprocessing…

…-2.4.2

Update to pygeoprocessing 2.4.2

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.

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

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,24 +764,25 @@ 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')
     dependent_tasks_for_watersheds_list.append(calculate_wyield_task)
 
     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.
 

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):