From 87a99d8fb95c1cec43bef5c5299558adf5438de2 Mon Sep 17 00:00:00 2001 From: Brett Edwards Date: Wed, 9 Oct 2024 11:52:56 -0700 Subject: [PATCH] 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 | 4 +- api/app/tests/utils/test_geospatial.py | 4 +- api/app/utils/geospatial-interpolation.md | 49 ++++++++++++++++++++++ api/app/utils/geospatial.py | 21 ++++++++-- 4 files changed, 71 insertions(+), 7 deletions(-) create mode 100644 api/app/utils/geospatial-interpolation.md diff --git a/api/app/auto_spatial_advisory/elevation.py b/api/app/auto_spatial_advisory/elevation.py index b75f5d344..74d7e2540 100644 --- a/api/app/auto_spatial_advisory/elevation.py +++ b/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 diff --git a/api/app/tests/utils/test_geospatial.py b/api/app/tests/utils/test_geospatial.py index dc82ec3a4..7a88b3f70 100644 --- a/api/app/tests/utils/test_geospatial.py +++ b/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() diff --git a/api/app/utils/geospatial-interpolation.md b/api/app/utils/geospatial-interpolation.md new file mode 100644 index 000000000..61d0b5189 --- /dev/null +++ b/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. diff --git a/api/app/utils/geospatial.py b/api/app/utils/geospatial.py index a80ebca5d..f28fe2662 100644 --- a/api/app/utils/geospatial.py +++ b/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: