Warn about geographic CRS in raster_band_percentile

HISTORY.rst

@@ -6,6 +6,9 @@ Release History
   to the target if doing so would raise a RuntimeError,
   such as when a string value cannot be represented by UTF-8.
   https://github.com/natcap/pygeoprocessing/issues/418
+* ``raster_band_percentile`` can now optionally log a warning if the raster
+  has a geographic CRS.
+  https://github.com/natcap/pygeoprocessing/issues/299
 
 2.4.7 (2025-01-23)
 ------------------

src/pygeoprocessing/geoprocessing_core.pyx

@@ -681,8 +681,8 @@ ctypedef FastFileIterator[double]* FastFileIteratorDoublePtr
 
 def raster_band_percentile(
         base_raster_path_band, working_sort_directory, percentile_list,
-        heap_buffer_size=2**28, ffi_buffer_size=2**10):
-    """Calculate percentiles of a raster band.
+        heap_buffer_size=2**28, ffi_buffer_size=2**10, geographic_crs_warn=False):
+    """Calculate percentiles of a raster band based on pixel values.
 
     Parameters:
         base_raster_path_band (tuple): raster path band tuple to a raster
@@ -700,6 +700,8 @@ def raster_band_percentile(
             disk.
         ffi_buffer_size (int): defines how many elements will be stored per
             heap file buffer for iteration.
+        geographic_crs_warn (boolean): defaults to False. If True, a warning will
+            be issued if the base raster has a geographic CRS.
 
     Returns:
         A list of len(percentile_list) elements long containing the
@@ -708,6 +710,18 @@ def raster_band_percentile(
         will select the next element higher than the percentile cutoff).
 
     """
+    if geographic_crs_warn:
+        try:
+            base_raster = gdal.OpenEx(base_raster_path_band[0], gdal.OF_RASTER)
+            srs = base_raster.GetSpatialRef()
+            if srs.IsGeographic():
+                LOGGER.warning(
+                    f'Raster {base_raster_path_band[0]} has a geographic CRS. '
+                    'Because `raster_band_percentile` calculates percentiles '
+                    'of pixel values, percentile results may be skewed.')
+        finally:
+            base_raster = None
+
     numpy_type = pygeoprocessing.get_raster_info(
         base_raster_path_band[0])['numpy_type']
     if numpy.issubdtype(numpy_type, numpy.integer):

tests/test_geoprocessing.py

@@ -4444,6 +4444,42 @@ def test_assert_is_valid_pixel_size(self):
         actual_message = str(cm.exception)
         self.assertIn(expected_message, actual_message)
 
+    def test_raster_band_percentile_warning(self):
+        """PGP: test raster_band_percentile geographic CRS warning."""
+        geo_raster_path = os.path.join(self.workspace_dir, 'geo_raster.tif')
+        _array_to_raster(
+            numpy.ones((10,10)), 0, geo_raster_path,
+            projection_epsg=4326) # Geographic CRS
+        proj_raster_path = os.path.join(self.workspace_dir, 'proj_raster.tif')
+        _array_to_raster(
+            numpy.ones((10,10)), 0, proj_raster_path,
+            projection_epsg=3116) # Projected CRS
+
+        percentile_cutoffs = [0.0]*5
+        working_dir = os.path.join(
+            self.workspace_dir, 'percentile_working_dir')
+
+        with capture_logging(
+                logging.getLogger('pygeoprocessing')) as log_messages:
+            pygeoprocessing.raster_band_percentile(
+                (geo_raster_path, 1), working_dir, percentile_cutoffs,
+                heap_buffer_size=8, ffi_buffer_size=4,
+                geographic_crs_warn=True)
+        self.assertEqual(len(log_messages), 1)
+        self.assertEqual(log_messages[0].levelno, logging.WARNING)
+        self.assertIn('geographic CRS', log_messages[0].msg)
+
+        with capture_logging(
+                logging.getLogger('pygeoprocessing')) as log_messages:
+            pygeoprocessing.raster_band_percentile(
+                (proj_raster_path, 1), working_dir, percentile_cutoffs,
+                heap_buffer_size=8, ffi_buffer_size=4,
+                geographic_crs_warn=True)
+        self.assertEqual(len(log_messages), 0)
+
+        self.assertTrue(
+            not os.path.exists(working_dir), 'working dir was not deleted')
+
     def test_percentile_int_type(self):
         """PGP: test percentile with int type."""
         srs = osr.SpatialReference()