raster_calculator stats to include percentage of valid pixels

HISTORY.rst

@@ -1,6 +1,12 @@
 Release History
 ===============
 
+.. Unreleased Changes
+* ``raster_calculator`` ``calc_raster_stats=True`` will count the number
+  of non-nodata pixels in the target raster and write
+  'STATISTICS_VALID_PERCENT' metadata along with other stats.
+  https://github.com/natcap/pygeoprocessing/issues/431
+
 2.4.8 (2025-05-02)
 ------------------
 * ``zonal_statistics`` raises an error if the vector is not a Polygon

docs/environment-rtd.yml

@@ -8,7 +8,7 @@ name: env-readthedocs
 channels:
     - conda-forge
 dependencies:
-    - python=3.8
+    - python=3.12
     - gdal>=3
     - numpy
     - rtree

src/pygeoprocessing/geoprocessing.py

@@ -435,6 +435,7 @@ def raster_calculator(
 
             pixels_processed = 0
             n_pixels = n_cols * n_rows
+            valid_pixel_count = 0
 
             # iterate over each block and calculate local_op
             for block_offset in block_offset_list:
@@ -489,7 +490,9 @@ def raster_calculator(
                 if stats_worker_queue:
                     # guard against an undefined nodata target
                     if nodata_target is not None:
-                        target_block = target_block[target_block != nodata_target]
+                        valid_mask = ~array_equals_nodata(target_block, nodata_target)
+                        valid_pixel_count += valid_mask.sum()
+                        target_block = target_block[valid_mask]
                     target_block = target_block.astype(numpy.float64).flatten()
 
                     if sys.version_info >= (3, 8) and use_shared_memory:
@@ -529,6 +532,9 @@ def raster_calculator(
                     target_band.SetStatistics(
                         float(target_min), float(target_max), float(target_mean),
                         float(target_stddev))
+                    target_band.SetMetadataItem(
+                        'STATISTICS_VALID_PERCENT',
+                        f'{(valid_pixel_count / n_pixels * 100):.2f}')
                     target_band.FlushCache()
         except Exception:
             LOGGER.exception('exception encountered in raster_calculator')

tests/test_geoprocessing.py

@@ -2373,6 +2373,36 @@ def test_raster_calculator(self):
                 pygeoprocessing.raster_to_numpy_array(base_path),
                 pygeoprocessing.raster_to_numpy_array(target_path)).all())
 
+    def test_raster_calculator_stats(self):
+        """PGP.geoprocessing: raster_calculator test stats are saved."""
+        blocksize = 5
+        pixel_matrix = numpy.ones((blocksize*2, blocksize*2), numpy.int16)
+        target_nodata = -1
+        pixel_matrix[:1] = target_nodata
+        base_path = os.path.join(self.workspace_dir, 'base.tif')
+        _array_to_raster(pixel_matrix, target_nodata, base_path)
+
+        target_path = os.path.join(self.workspace_dir, 'subdir', 'target.tif')
+        # Test with multiple blocks by limiting largest_block
+        pygeoprocessing.raster_calculator(
+            [(base_path, 1)], passthrough, target_path,
+            gdal.GDT_Int32, target_nodata, calc_raster_stats=True,
+            use_shared_memory=True, largest_block=blocksize)
+
+        raster = gdal.OpenEx(target_path)
+        band = raster.GetRasterBand(1)
+        metadata = band.GetMetadata()
+        raster = band = None
+        expected_metadata = {  # gdal metadata are represented as strings
+            'STATISTICS_MINIMUM': '1',
+            'STATISTICS_MAXIMUM': '1',
+            'STATISTICS_MEAN': '1',
+            'STATISTICS_STDDEV': '0',
+            'STATISTICS_VALID_PERCENT': '90.00'}
+        for key, value in expected_metadata.items():
+            self.assertIn(key, metadata)
+            self.assertEqual(metadata[key], value)
+
     def test_raster_calculator_multiprocessing(self):
         """PGP.geoprocessing: raster_calculator identity test."""
         pixel_matrix = numpy.ones((1024, 1024), numpy.int16)