add resampling method

bcgov · Oct 7, 2024 · b8b88b5 · b8b88b5
1 parent 8aefeda
commit b8b88b5
Show file tree

Hide file tree

Showing 3 changed files with 12 additions and 8 deletions.
diff --git 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
@@ -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.py b/api/app/utils/geospatial.py
@@ -1,18 +1,22 @@
 import logging
-from typing import Tuple
+from typing import Tuple, Literal
 from osgeo import gdal, ogr, osr
 
 
 logger = logging.getLogger(__name__)
 
 
-def warp_to_match_extent(source_ds: gdal.Dataset, ds_to_match: gdal.Dataset, output_path: str) -> gdal.Dataset:
+GDALResamplingMethod = Literal[gdal.GRA_NearestNeighbour, gdal.GRA_Bilinear, gdal.GRA_Cubic]
+
+
+def warp_to_match_raster(source_ds: gdal.Dataset, ds_to_match: gdal.Dataset, output_path: str, resample_method: GDALResamplingMethod = gdal.GRA_NearestNeighbour) -> gdal.Dataset:  # type: ignore
     """
-    Warp the source dataset to match the extent and projection of the other 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 +29,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)
 
 
 def raster_mul(tpi_ds: gdal.Dataset, hfi_ds: gdal.Dataset, chunk_size=256) -> gdal.Dataset: