Generic raster resampling (#3988)

- Adds resampling method to already existing raster extent matcher
- Closes SFMS: Write generic grid resampling functions #3665

api/app/auto_spatial_advisory/elevation.py

@@ -20,7 +20,7 @@
 from app.db.models.auto_spatial_advisory import AdvisoryElevationStats, AdvisoryTPIStats
 from app.auto_spatial_advisory.hfi_filepath import get_raster_filepath, get_raster_tif_filename
 from app.utils.s3 import get_client
-from app.utils.geospatial import raster_mul, warp_to_match_extent
+from app.utils.geospatial import raster_mul, warp_to_match_raster
 
 
 logger = logging.getLogger(__name__)
@@ -253,7 +253,7 @@ async def process_tpi_by_firezone(run_type: RunType, run_date: date, for_date: d
     hfi_source: gdal.Dataset = gdal.Open(hfi_key, gdal.GA_ReadOnly)
 
     warped_mem_path = f"/vsimem/warp_{hfi_raster_filename}"
-    resized_hfi_source: gdal.Dataset = warp_to_match_extent(hfi_source, tpi_source, warped_mem_path)
+    resized_hfi_source: gdal.Dataset = warp_to_match_raster(hfi_source, tpi_source, warped_mem_path)
     hfi_masked_tpi = raster_mul(tpi_source, resized_hfi_source)
     resized_hfi_source = None
     hfi_source = None

api/app/tests/utils/test_geospatial.py

@@ -3,7 +3,7 @@
 from osgeo import gdal
 import numpy as np
 
-from app.utils.geospatial import raster_mul, warp_to_match_extent
+from app.utils.geospatial import raster_mul, warp_to_match_raster
 
 fixture_path = os.path.join(os.path.dirname(__file__), "snow_masked_hfi20240810.tif")
 
@@ -96,7 +96,7 @@ def test_warp_to_match_dimension():
     driver = gdal.GetDriverByName("MEM")
     out_dataset: gdal.Dataset = driver.Create("memory", hfi_ds.RasterXSize, hfi_ds.RasterYSize, 1, gdal.GDT_Byte)
 
-    warp_to_match_extent(tpi_ds, hfi_ds, out_dataset)
+    warp_to_match_raster(tpi_ds, hfi_ds, out_dataset)
     output_data = out_dataset.GetRasterBand(1).ReadAsArray()
     hfi_data = hfi_ds.GetRasterBand(1).ReadAsArray()
 

api/app/utils/geospatial-interpolation.md

@@ -0,0 +1,49 @@
+# Raster Resampling Methods
+
+When resampling or reprojecting a raster dataset using GDAL, different interpolation methods can be applied based on the use case. The interpolation method determines how pixel values are calculated when transforming a raster.
+
+## Resampling Methods
+
+### 1. **NEAREST_NEIGHBOUR (`gdal.GRA_NearestNeighbour`)**
+
+#### Description
+
+- Nearest neighbour interpolation takes the value from the closest pixel to the new pixel location without any modification.
+
+#### Use Cases
+
+- **Discrete data**: Best for categorical or discrete datasets (e.g., land cover classification, thematic maps).
+- **Maintains original values**: Since it doesn't modify pixel values, it's ideal when you need to preserve the integrity of original values (e.g., for classes of a fuel grid raster).
+
+---
+
+### 2. **BILINEAR (`gdal.GRA_Bilinear`)**
+
+#### Description
+
+- Bilinear interpolation calculates the new pixel value by taking a weighted average of the four nearest neighboring pixels.
+
+#### Use Cases
+
+- **Continuous data**: Appropriate for resampling continuous variables such as elevation or weather data (temp, precip, rh, wind speed).
+- **Alters original values**: Alters the original values (which can be undesirable for discrete datasets).
+
+---
+
+### 3. **CUBIC (`gdal.GRA_Cubic`)**
+
+#### Description
+
+- Cubic interpolation uses 16 surrounding pixels to calculate a new pixel value using cubic convolution. This produces smoother results than bilinear interpolation.
+
+#### Use Cases
+
+- **Continuous data**: Ideal for datasets where smoothness and visual quality are important (e.g., satellite imagery, elevation data).
+- **Alters original values**: Alters the original values (which can be undesirable for discrete datasets).
+
+---
+
+## How to Choose the Right Resampling Method
+
+- **For categorical or discrete data** (e.g., land cover, fuel grid classification): Use **NEAREST_NEIGHBOUR** to preserve the integrity of the original pixel values.
+- **For continuous data** (e.g., elevation, weather data): We use **BILINEAR** to limit the smoothing and "spreading" of weather data from the original pixel. **CUBIC** is another option that could be explored.

api/app/utils/geospatial.py

@@ -1,3 +1,4 @@
+from enum import Enum
 import logging
 from typing import Tuple
 from osgeo import gdal, ogr, osr
@@ -6,13 +7,27 @@
 logger = logging.getLogger(__name__)
 
 
-def warp_to_match_extent(source_ds: gdal.Dataset, ds_to_match: gdal.Dataset, output_path: str) -> gdal.Dataset:
+class GDALResamplingMethod(Enum):
     """
-    Warp the source dataset to match the extent and projection of the other dataset.
+    See api/app/utils/geospatial-interpolation.md for information about which interpolation method to use for your use case
+
+    """
+
+    NEAREST_NEIGHBOUR = gdal.GRA_NearestNeighbour
+    BILINEAR = gdal.GRA_Bilinear
+    CUBIC = gdal.GRA_Cubic
+
+
+def warp_to_match_raster(
+    source_ds: gdal.Dataset, ds_to_match: gdal.Dataset, output_path: str, resample_method: GDALResamplingMethod = GDALResamplingMethod.NEAREST_NEIGHBOUR
+) -> gdal.Dataset:
+    """
+    Warp the source dataset to match the extent, pixel size, and projection of the other dataset.
 
     :param source_ds: the dataset raster to warp
     :param ds_to_match: the reference dataset raster to match the source against
     :param output_path: output path of the resulting raster
+    :param resample_method: gdal resampling algorithm
     :return: warped raster dataset
     """
     source_geotransform = ds_to_match.GetGeoTransform()
@@ -25,7 +40,7 @@ def warp_to_match_extent(source_ds: gdal.Dataset, ds_to_match: gdal.Dataset, out
     extent = [minx, miny, maxx, maxy]
 
     # Warp to match input option parameters
-    return gdal.Warp(output_path, source_ds, dstSRS=ds_to_match.GetProjection(), outputBounds=extent, xRes=x_res, yRes=y_res, resampleAlg=gdal.GRA_NearestNeighbour)
+    return gdal.Warp(output_path, source_ds, dstSRS=ds_to_match.GetProjection(), outputBounds=extent, xRes=x_res, yRes=y_res, resampleAlg=resample_method.value)
 
 
 def raster_mul(tpi_ds: gdal.Dataset, hfi_ds: gdal.Dataset, chunk_size=256) -> gdal.Dataset: