Handle rasters without a defined NoData value in align_and_resize_raster_stack

HISTORY.rst

@@ -26,6 +26,11 @@ Unreleased Changes
 * Fixed a bug where ``convolve_2d`` was not handling ``nan`` NoData
   values correctly.
   https://github.com/natcap/pygeoprocessing/issues/473
+* Fixed a bug where ``align_and_resize_raster_stack`` would pad rasters that
+  were smaller than the target extent and didn't have a defined NoData value
+  with 0s. The function will now automatically assign an appropriate NoData
+  value to any input raster without one defined.
+  https://github.com/natcap/pygeoprocessing/issues/476
 
 2.4.10 (2026-01-13)
 -------------------

src/pygeoprocessing/geoprocessing.py

@@ -1102,6 +1102,44 @@ def align_and_resize_raster_stack(
                 n_pixels * align_pixel_size[index] +
                 align_bounding_box[index])
 
+    # Set a nodata value if the raster does not have one defined. This
+    # ensures that a raster that is smaller than the target extent will
+    # not get filled with 0s that could be confused with real data.
+    # We do this by creating a temporary VRT with an explicit nodata value
+    # for rasters that are missing one.
+    fixed_base_raster_path_list = list(base_raster_path_list)
+    temp_working_dir_nodata = None
+    for index, (raster_info, base_path) in enumerate(
+                zip(raster_info_list, base_raster_path_list)):
+        nodata_list = raster_info['nodata']
+        if None in nodata_list:
+            if temp_working_dir_nodata is None:
+                temp_working_dir_nodata = tempfile.mkdtemp(dir=working_dir,
+                                               prefix='align-nodata-')
+            raster = gdal.OpenEx(base_path, gdal.OF_RASTER)
+            chosen_nodata = choose_nodata(raster_info['numpy_type'])
+
+            if set(nodata_list) != {None}:
+                LOGGER.warning(
+                    "Input raster %s has some bands with a defined nodata "
+                    "value and some bands without. Setting %s as the nodata "
+                    "value for all bands in a temporary VRT before aligning "
+                    "and resizing.", base_path, chosen_nodata)
+            else:
+                LOGGER.warning(
+                    "Input raster %s has no defined nodata value. Setting "
+                    "%s as nodata in a temporary VRT before aligning and "
+                    "resizing.", base_path, chosen_nodata)
+
+            vrt_path = os.path.join(temp_working_dir_nodata,
+                                    f'input_with_nodata_{index}.vrt')
+            vrt = gdal.Translate(vrt_path, raster, format='VRT',
+                                 noData=chosen_nodata)
+            vrt = None
+            raster = None
+
+            fixed_base_raster_path_list[index] = vrt_path
+
     if mask_options:
         # Create a warped VRT.
         # This allows us to cheaply figure out the dimensions, projection, etc.
@@ -1110,7 +1148,7 @@ def align_and_resize_raster_stack(
         temp_working_dir = tempfile.mkdtemp(dir=working_dir, prefix='mask-')
         warped_vrt = os.path.join(temp_working_dir, 'warped.vrt')
         warp_raster(
-            base_raster_path=base_raster_path_list[0],
+            base_raster_path=fixed_base_raster_path_list[0],
             target_pixel_size=target_pixel_size,
             target_raster_path=warped_vrt,
             resample_method='near',
@@ -1142,7 +1180,7 @@ def align_and_resize_raster_stack(
                       else None))
 
     for index, (base_path, target_path, resample_method) in enumerate(zip(
-            base_raster_path_list, target_raster_path_list,
+            fixed_base_raster_path_list, target_raster_path_list,
             resample_method_list)):
         warp_raster(
             base_path, target_pixel_size, target_path, resample_method,
@@ -1160,6 +1198,8 @@ def align_and_resize_raster_stack(
 
     LOGGER.info("aligned all %d rasters.", n_rasters)
 
+    if temp_working_dir_nodata:
+        shutil.rmtree(temp_working_dir_nodata, ignore_errors=True)
     if mask_options:
         shutil.rmtree(temp_working_dir, ignore_errors=True)
 

tests/test_geoprocessing.py

@@ -2585,6 +2585,78 @@ def test_align_and_resize_raster_stack_bb(self):
         numpy.testing.assert_array_equal(
             numpy.ones((4, 4), pixel_a_matrix.dtype), target_array)
 
+    def test_align_and_resize_raster_stack_no_set_nodata_multiband(self):
+        """Test align_and_resize sets nodata if not defined for input."""
+        arr = numpy.arange(16).reshape(4, 4)
+
+        srs = osr.SpatialReference()
+        srs.ImportFromEPSG(26910)
+        input_raster_path = os.path.join(self.workspace_dir, 'test.tif')
+        pygeoprocessing.numpy_array_to_raster(
+            arr.astype(numpy.int16), target_nodata=None,
+            pixel_size=(10, -10), origin=(0, 0),
+            projection_wkt=srs.ExportToWkt(),
+            target_path=input_raster_path)
+
+        output_raster_path = os.path.join(self.workspace_dir, 'test_align.tif')
+        pygeoprocessing.align_and_resize_raster_stack(
+            [input_raster_path],
+            [output_raster_path],
+            ['near'],
+            (10, -10),
+            [-10, -20, 20, 0])
+            # this target bounding box extends beyond extent of input raster
+            # (input raster's bbox is [0, -40, 40, 0])
+
+        output_array = pygeoprocessing.raster_to_numpy_array(
+            output_raster_path)
+        expected_output_array = numpy.array([[32767, 0, 1], [32767, 4, 5]])
+        numpy.testing.assert_allclose(output_array, expected_output_array)
+
+    def test_align_and_resize_raster_stack_mixed_nodata_warns(self):
+        """Test warning if 1 band in multiband raster has no defined nodata.
+
+        Test that a warning is raised if a multiband raster with one band
+        that has a defined nodata value and one band that does not have a
+        defined nodata value is passed to align_and_resize_raster_stack.
+        """
+
+        raster_path = os.path.join(self.workspace_dir, 'multi.tif')
+        driver = gdal.GetDriverByName('GTiff')
+        raster = driver.Create(raster_path, 4, 4, 2, gdal.GDT_Int16)
+        srs = osr.SpatialReference()
+        srs.ImportFromEPSG(26910)
+        raster.SetProjection(srs.ExportToWkt())
+        raster.SetGeoTransform((0, 10, 0, 0, 0, -10))
+        arr = numpy.arange(16).reshape(4, 4)
+        band1 = raster.GetRasterBand(1)
+        band2 = raster.GetRasterBand(2)
+        band1.WriteArray(arr)
+        band2.WriteArray(arr)
+        raster = None
+
+        vrt_path = os.path.join(self.workspace_dir, 'multi.vrt')
+        gdal.BuildVRT(vrt_path, [raster_path])
+        vrt = gdal.Open(vrt_path, gdal.GA_Update)
+        vrt.GetRasterBand(2).SetNoDataValue(-9999)
+        # leave band 1 without nodata
+        vrt = None
+
+        output_raster_path = os.path.join(self.workspace_dir, 'out.tif')
+        with self.assertLogs(level='WARNING') as cm:
+            pygeoprocessing.align_and_resize_raster_stack(
+                [vrt_path],
+                [output_raster_path],
+                ['near'],
+                (10, -10),
+                [-10, -20, 20, 0]
+            )
+
+        self.assertTrue(
+            any("some bands with a defined nodata value and some bands without"
+                in msg for msg in cm.output),
+            f"Expected warning not found in logs: {cm.output}")
+
     def test_raster_calculator(self):
         """PGP.geoprocessing: raster_calculator identity test."""
         pixel_matrix = numpy.ones((5, 5), numpy.int16)