Merge branch 'bugfix/254-watershed-delineation-requires-srs' of githu…

…b.com:phargogh/pygeoprocessing into bugfix/254-watershed-delineation-requires-srs

.github/workflows/pythonpackage.yml

@@ -62,19 +62,21 @@ jobs:
           fetch-depth: 0
 
       - name: setup-micromamba
-        uses: mamba-org/provision-with-micromamba@main
+        uses: mamba-org/setup-micromamba@v1
         with:
           # Grab requirements from pip-compatible requirements.txt
           environment-file: requirements.txt
-          extra-specs: |
+          condarc: |
+              channels:
+                - conda-forge
+          create-args: >-
               python=${{ matrix.python-version }}
               gdal=${{ matrix.gdal }}
               setuptools
               python-build
               flake8
               pytest
           environment-name: pyenv
-          channels: conda-forge
 
       - name: Lint with flake8
         run: |

HISTORY.rst

@@ -23,6 +23,20 @@ Unreleased Changes
 * ``pygeoprocessing.routing.delineate_watersheds_d8`` now handles the case
   where the input flow direction raster does not have a defined spatial
   reference. https://github.com/natcap/pygeoprocessing/issues/254
+* Updating internal documentation describing TauDEM flow directions, and adding
+  for how to convert from a flow direction raster from what TauDEM expects to
+  what pygeoprocessing expects.
+  https://github.com/natcap/pygeoprocessing/issues/255
+* Users may now specify the overview level to use when calling ``warp_raster``.
+  By default, ``pygeoprocessing`` will use the base layer.
+  https://github.com/natcap/pygeoprocessing/issues/326
+* Fixed a bug across ``pygeoprocessing`` where some valid resampling methods
+  would throw an exception because they were not recognized.  This was only
+  happening when ``pygeoprocessing`` was installed alongside GDAL < 3.4.
+* Fixing an issue with ``pygeoprocessing.multiprocessing.raster_calculator``
+  where the function would raise an Exception when the target raster path was
+  provided as a filename only, not within a directory, even though the parent
+  directory could be inferred. https://github.com/natcap/pygeoprocessing/issues/313
 
 2.4.0 (2023-03-03)
 ------------------

src/pygeoprocessing/geoprocessing.py

@@ -7,6 +7,7 @@
 import os
 import pprint
 import queue
+import re
 import shutil
 import sys
 import tempfile
@@ -30,10 +31,10 @@
 from osgeo import osr
 
 from . import geoprocessing_core
-from .geoprocessing_core import DEFAULT_GTIFF_CREATION_TUPLE_OPTIONS
 from .geoprocessing_core import DEFAULT_CREATION_OPTIONS
-from .geoprocessing_core import INT8_CREATION_OPTIONS
+from .geoprocessing_core import DEFAULT_GTIFF_CREATION_TUPLE_OPTIONS
 from .geoprocessing_core import DEFAULT_OSR_AXIS_MAPPING_STRATEGY
+from .geoprocessing_core import INT8_CREATION_OPTIONS
 
 # This is used to efficiently pass data to the raster stats worker if available
 if sys.version_info >= (3, 8):
@@ -97,19 +98,49 @@ def __init__(self, missing_values, raster_path, value_map):
     gdal.GDT_CFloat64: numpy.complex64,
 }
 
-# GDAL's python API recognizes certain strings but the only way to retrieve
-# those strings is to do this conversion of gdalconst.GRA_* types to the
-# human-readable labels via gdal.WarpOptions like so.
 _GDAL_WARP_ALGORITHMS = []
 for _warp_algo in (_attrname for _attrname in dir(gdalconst)
                    if _attrname.startswith('GRA_')):
-    # Appends ['-r', 'near'] to _GDAL_WARP_ALGORITHMS if _warp_algo is
-    # gdalconst.GRA_NearestNeighbor.  See gdal.WarpOptions for the dict
-    # defining this mapping.
-    gdal.WarpOptions(options=_GDAL_WARP_ALGORITHMS,
-                     resampleAlg=getattr(gdalconst, _warp_algo))
+    # In GDAL < 3.4, the options used were actually gdal.GRIORA_*
+    # instead of gdal.GRA_*, and gdal.WarpOptions would:
+    #   * return an integer for any gdal.GRIORA options it didn't recognize
+    #   * Return an abbreviated string (e.g. -rb instead of 'bilinear') for
+    #     several warp algorithms
+    # This behavior was changed in GDAL 3.4, where the correct name is
+    # returned in all cases.
+    #
+    # In GDAL < 3.4, an error would be logged if GDAL couldn't recognize the
+    # option.  Pushing a null error handler avoids this and cleans up the
+    # logging.
+    gdal.PushErrorHandler(lambda *args: None)
+    _options = []
+    # GDAL's python bindings only offer this one way of translating gdal.GRA_*
+    # options to their string names (but compatibility is limited in different
+    # GDAL versions, see notes above)
+    warp_opts = gdal.WarpOptions(
+        options=_options,  # SIDE EFFECT: Adds flags to this list.
+        resampleAlg=getattr(gdalconst, _warp_algo),  # use GRA_* name.
+        overviewLevel=None)  # Don't add overview parameters.
+    gdal.PopErrorHandler()
+
+    # _options is populated during the WarpOptions call above.
+    for _item in _options:
+        is_int = False
+        try:
+            int(_item)
+            is_int = True
+        except ValueError:
+            pass
+
+        if _item == '-r':
+            continue
+        elif (_item.startswith('-r') and len(_item) > 2) or is_int:
+            # (GDAL < 3.4) Translate shorthand params to name.
+            # (GDAL < 3.4) Translate unrecognized (int) codes to name.
+            _item = re.sub('^gra_', '', _warp_algo.lower())
+        _GDAL_WARP_ALGORITHMS.append(_item)
+
 _GDAL_WARP_ALGORITHMS = set(_GDAL_WARP_ALGORITHMS)
-_GDAL_WARP_ALGORITHMS.discard('-r')
 _GDAL_WARP_ALGOS_FOR_HUMAN_EYES = "|".join(_GDAL_WARP_ALGORITHMS)
 LOGGER.debug(
     f'Detected warp algorithms: {", ".join(_GDAL_WARP_ALGOS_FOR_HUMAN_EYES)}')
@@ -2356,7 +2387,7 @@ def warp_raster(
         base_raster_path, target_pixel_size, target_raster_path,
         resample_method, target_bb=None, base_projection_wkt=None,
         target_projection_wkt=None, n_threads=None, vector_mask_options=None,
-        gdal_warp_options=None, working_dir=None,
+        gdal_warp_options=None, working_dir=None, use_overview_level=-1,
         raster_driver_creation_tuple=DEFAULT_GTIFF_CREATION_TUPLE_OPTIONS,
         osr_axis_mapping_strategy=DEFAULT_OSR_AXIS_MAPPING_STRATEGY):
     f"""Resize/resample raster to desired pixel size, bbox and projection.
@@ -2402,6 +2433,15 @@ def warp_raster(
         working_dir (string): if defined uses this directory to make
             temporary working files for calculation. Otherwise uses system's
             temp directory.
+        use_overview_level=-1 (int/str): The overview level to use for warping.
+            A value of ``-1`` (the default) indicates that the base raster
+            should be used for the source pixels. A value of ``'AUTO'``
+            will make GDAL select the overview with the resolution that is
+            closest to the target pixel size and warp using that overview's
+            pixel values.  Any other integer indicates that that overview index
+            should be used.  For example, suppose the raster has overviews at
+            levels 2, 4 and 8.  To use level 2, set ``use_overview_level=0``.
+            To use level 8, set ``use_overview_level=2``.
         raster_driver_creation_tuple (tuple): a tuple containing a GDAL driver
             name string as the first element and a GDAL creation options
             tuple/list as the second. Defaults to a GTiff driver tuple
@@ -2540,6 +2580,7 @@ def warp_raster(
         dstSRS=target_projection_wkt,
         multithread=True if warp_options else False,
         warpOptions=warp_options,
+        overviewLevel=use_overview_level,
         creationOptions=raster_creation_options,
         callback=reproject_callback,
         callback_data=[target_raster_path])

src/pygeoprocessing/multiprocessing/raster_calculator.py

@@ -9,6 +9,9 @@
 import sys
 import time
 
+import numpy
+from osgeo import gdal
+
 from ..geoprocessing import _is_raster_path_band_formatted
 from ..geoprocessing import _LARGEST_ITERBLOCK
 from ..geoprocessing import _LOGGING_PERIOD
@@ -19,8 +22,6 @@
 from ..geoprocessing import iterblocks
 from ..geoprocessing import LOGGER
 from ..geoprocessing_core import DEFAULT_GTIFF_CREATION_TUPLE_OPTIONS
-from osgeo import gdal
-import numpy
 
 if sys.version_info >= (3, 8):
     import multiprocessing.shared_memory
@@ -451,6 +452,7 @@ def raster_calculator(
         for path_band in base_raster_path_band_const_list
         if _is_raster_path_band_formatted(path_band)]
 
+    target_raster_path = os.path.abspath(target_raster_path)
     _validate_raster_input(
         base_raster_path_band_const_list, raster_info_list, target_raster_path)
 

src/pygeoprocessing/routing/routing.pyx

@@ -10,12 +10,24 @@ incoming DEM type is an int64 type and values in that dem exceed 2^52 but GDAL
 does not support int64 rasters so no precision loss is possible with a
 float64.
 
-D8 float direction conventions follow TauDEM where each flow direction
-is encoded as::
+D8 flow direction conventions encode the flow direction as::
 
      3 2 1
      4 x 0
      5 6 7
+
+This is slightly different from how TauDEM encodes flow direction, which is as::
+    4 3 2
+    5 x 1
+    6 7 8
+
+To convert a TauDEM flow direction raster to a pygeoprocessing-compatible flow
+direction raster, the following ``raster_map`` call may be used::
+
+    taudem_flow_dir_path = 'taudem_d8_flow_dir.tif'
+    pygeoprocessing.raster_map(
+        lambda d: d+1, [taudem_flow_dir_path],
+        'pygeoprocessing_d8_flow_dir.tif')
 """
 import collections
 import logging

tests/test_geoprocessing.py

@@ -29,10 +29,11 @@
 from osgeo import gdal
 from osgeo import ogr
 from osgeo import osr
+from pygeoprocessing.geoprocessing_core import DEFAULT_CREATION_OPTIONS
+from pygeoprocessing.geoprocessing_core import \
+    DEFAULT_GTIFF_CREATION_TUPLE_OPTIONS
+from pygeoprocessing.geoprocessing_core import INT8_CREATION_OPTIONS
 from pygeoprocessing.geoprocessing_core import \
-    DEFAULT_CREATION_OPTIONS, \
-    INT8_CREATION_OPTIONS, \
-    DEFAULT_GTIFF_CREATION_TUPLE_OPTIONS, \
     INT8_GTIFF_CREATION_TUPLE_OPTIONS
 
 _DEFAULT_ORIGIN = (444720, 3751320)
@@ -1693,6 +1694,54 @@ def test_warp_raster(self):
                 pygeoprocessing.raster_to_numpy_array(
                     target_raster_path)).all())
 
+    def test_warp_raster_overview_level(self):
+        """PGP.geoprocessing: warp raster overview test."""
+        # This array is big enough that build_overviews will render several
+        # overview levels.
+        pixel_a_matrix = numpy.array([
+            [1, 2, 3, 4],
+            [1, 2, 3, 4],
+            [1, 2, 3, 4],
+            [1, 2, 3, 4],
+        ], dtype=numpy.float32)
+
+        # Using overview level 0 (the base raster), we should have the same
+        # output when we warp the array that has and does not have overviews
+        # present.
+        target_nodata = -1
+        base_a_path = os.path.join(self.workspace_dir, 'base_a.tif')
+        _array_to_raster(
+            pixel_a_matrix, target_nodata, base_a_path)
+        base_a_raster_info = pygeoprocessing.get_raster_info(base_a_path)
+        warped_a_path = os.path.join(self.workspace_dir, 'warped_a.tif')
+        pygeoprocessing.warp_raster(
+            base_a_path, base_a_raster_info['pixel_size'], warped_a_path,
+            'bilinear', use_overview_level=-1)
+
+        base_b_path = os.path.join(self.workspace_dir, 'base_b.tif')
+        _array_to_raster(
+            pixel_a_matrix, target_nodata, base_b_path)
+        warped_b_path = os.path.join(self.workspace_dir, 'warped_b.tif')
+        pygeoprocessing.build_overviews(
+            base_b_path, levels=[2, 4], resample_method='bilinear')
+        pygeoprocessing.warp_raster(
+            base_b_path, base_a_raster_info['pixel_size'], warped_b_path,
+            'bilinear', use_overview_level=-1)
+
+        warped_a_array = pygeoprocessing.raster_to_numpy_array(warped_a_path)
+        warped_b_array = pygeoprocessing.raster_to_numpy_array(warped_b_path)
+        numpy.testing.assert_allclose(warped_a_array, warped_b_array)
+
+        # Force warping using a higher overview level.
+        # Overview level 2 really means the 2nd overview in the stack, which is
+        # where 4 pixels are aggregated into 1 value.
+        target_raster_path = os.path.join(self.workspace_dir, 'target_c.tif')
+        pygeoprocessing.warp_raster(
+            base_b_path, base_a_raster_info['pixel_size'], target_raster_path,
+            'bilinear', n_threads=1, use_overview_level=1)
+        array = pygeoprocessing.raster_to_numpy_array(target_raster_path)
+        numpy.testing.assert_allclose(array, 2.5)
+
     def test_warp_raster_invalid_resample_alg(self):
         """PGP.geoprocessing: error on invalid resample algorithm."""
         pixel_a_matrix = numpy.ones((5, 5), numpy.int16)
@@ -2234,6 +2283,29 @@ def test_raster_calculator_mutiprocessing(self):
                 arithmetic_wrangle(pixel_matrix),
                 pygeoprocessing.raster_to_numpy_array(target_path)).all())
 
+    def test_raster_calculator_mutiprocessing_cwd(self):
+        """PGP.geoprocessing: raster_calculator identity test in cwd."""
+        pixel_matrix = numpy.ones((1024, 1024), numpy.int16)
+        target_nodata = -1
+        try:
+            cwd = os.getcwd()
+            os.chdir(self.workspace_dir)
+            base_path = 'base.tif'
+            _array_to_raster(pixel_matrix, target_nodata, base_path)
+
+            target_path = 'target.tif'
+            pygeoprocessing.multiprocessing.raster_calculator(
+                [(base_path, 1)], arithmetic_wrangle, target_path,
+                gdal.GDT_Int32, target_nodata, calc_raster_stats=True,
+                use_shared_memory=True)
+
+            self.assertTrue(
+                numpy.isclose(
+                    arithmetic_wrangle(pixel_matrix),
+                    pygeoprocessing.raster_to_numpy_array(target_path)).all())
+        finally:
+            os.chdir(cwd)
+
     def test_raster_calculator_bad_target_type(self):
         """PGP.geoprocessing: raster_calculator bad target type value."""
         pixel_matrix = numpy.ones((5, 5), numpy.int16)