Merge branch 'task/1379' of https://github.com/emlys/invest into task…

…/1379

.github/workflows/auto-pr-from-main-into-releases.yml

@@ -39,7 +39,7 @@ jobs:
           PR_NUM=$(git log -1 --pretty=%B main | head -n1 | egrep -o '#[1-9][0-9]*' | sed 's|#||g')
 
           # Get the username of the person who opened up the last PR into `main`.
-          PR_USERNAME=$(hub pr show -f "%au" $PR_NUM)
+          PR_USERNAME=$(gh pr view --json author $PR_NUM | jq -j '.author.login')
           echo "Latest PR on main ($PR_NUM) was authored by $PR_USERNAME."
 
           if [ "$PR_USERNAME" = "github-actions" ]
@@ -70,13 +70,13 @@ jobs:
               # This PR will be assigned to $GITHUB_ACTOR, which should be
               # the person who merged the PR that caused this commit to be
               # created.  Others could of course be assigned later.
-              hub pull-request \
+              gh pr create \
                   --head $GITHUB_REPOSITORY:$SOURCE_BRANCH \
                   --base $GITHUB_REPOSITORY:$RELEASE_BRANCH \
                   --reviewer "$PR_USERNAME" \
-                  --assign "$PR_USERNAME" \
-                  --labels "auto" \
-                  --file $PR_BODY_FILE || ERRORSPRESENT=$(($ERRORSPRESENT | $?))
+                  --assignee "$PR_USERNAME" \
+                  --label "auto" \
+                  --body-file $PR_BODY_FILE || ERRORSPRESENT=$(($ERRORSPRESENT | $?))
           done
 
           if [[ $ERRORSPRESENT -gt 0 ]]

HISTORY.rst

@@ -52,6 +52,19 @@ Unreleased Changes
     * Validation of tables has been improved and standardized, which should
       result in more readable validation errors.
       (`#1379 <https://github.com/natcap/invest/issues/1379>`_)
+    * 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.
+    * Version metadata at import time is now fetched with
+      ``importlib.metadata`` instead of ``pkg_resources``.
+      (`#1442 <https://github.com/natcap/invest/issues/1442>`_)
+* Coastal Vulnerability
+    * Fixed a bug where the model would crash when processing a float type
+      bathymetry raster with no nodata value.
+      https://github.com/natcap/invest/issues/992
 * 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.
@@ -62,6 +75,9 @@ Unreleased Changes
     * Fixed an issue in NDR's effective retention where, on rasters with more
       than 2^31 pixels, the model would crash with an error relating to a
       negative (overflowed) index. https://github.com/natcap/invest/issues/1431
+* Pollination
+    * Replaced custom kernel implementation with ``pygeoprocessing.kernels``.
+      Convolution results may be slightly different (more accurate).
 * SDR
     * RKLS, USLE, avoided erosion, and avoided export rasters will now have
       nodata in streams (`#1415 <https://github.com/natcap/invest/issues/1415>`_)

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/__init__.py

@@ -1,19 +1,19 @@
 """init module for natcap.invest."""
+import importlib.metadata
 import logging
 import os
 import sys
 from gettext import translation
 
 import babel
-import pkg_resources
 
 LOGGER = logging.getLogger('natcap.invest')
 LOGGER.addHandler(logging.NullHandler())
 __all__ = ['local_dir', ]
 
 try:
-    __version__ = pkg_resources.get_distribution(__name__).version
-except pkg_resources.DistributionNotFound:
+    __version__ = importlib.metadata.version('natcap.invest')
+except importlib.metadata.PackageNotFoundError:
     # package is not installed.  Log the exception for debugging.
     LOGGER.exception('Could not load natcap.invest version information')
 

src/natcap/invest/annual_water_yield.py

@@ -734,12 +734,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')
@@ -765,24 +765,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,
@@ -867,6 +868,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.
 
@@ -978,55 +983,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):
@@ -1063,19 +1019,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
@@ -1122,30 +1078,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.