Merge pull request #101 from wri/CIF-316-Add-bilinear-interpolation-f…

…or-DEM-DSM-and-albedo-layers

Added resampling method for 3 layers and updated tests

city_metrix/layers/albedo.py

@@ -2,7 +2,7 @@
 import xarray
 from dask.diagnostics import ProgressBar
 
-from .layer import Layer, get_utm_zone_epsg, get_image_collection
+from .layer import Layer, get_utm_zone_epsg, get_image_collection, set_resampling_for_continuous_raster
 
 
 class Albedo(Layer):
@@ -11,17 +11,20 @@ class Albedo(Layer):
         start_date: starting date for data retrieval
         end_date: ending date for data retrieval
         spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster)
+        resampling_method: interpolation method used by Google Earth Engine. Default is 'bilinear'. All options are: ('bilinear', 'bicubic', None).
         threshold: threshold value for filtering the retrieval
     """
 
-    def __init__(self, start_date="2021-01-01", end_date="2022-01-01", spatial_resolution=10, threshold=None, **kwargs):
+    def __init__(self, start_date="2021-01-01", end_date="2022-01-01", spatial_resolution:int=10,
+                 resampling_method:str='bilinear', threshold=None, **kwargs):
         super().__init__(**kwargs)
         self.start_date = start_date
         self.end_date = end_date
         self.spatial_resolution = spatial_resolution
+        self.resampling_method = resampling_method
         self.threshold = threshold
 
-    def get_data(self, bbox):
+    def get_data(self, bbox: tuple[float, float, float, float]):
         S2 = ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED")
         S2C = ee.ImageCollection("COPERNICUS/S2_CLOUD_PROBABILITY")
 
@@ -118,7 +121,17 @@ def calc_s2_albedo(image):
         ## S2 MOSAIC AND ALBEDO
         dataset = get_masked_s2_collection(ee.Geometry.BBox(*bbox), self.start_date, self.end_date)
         s2_albedo = dataset.map(calc_s2_albedo)
-        albedo_mean = s2_albedo.reduce(ee.Reducer.mean())
+
+        albedo_mean = (s2_albedo
+                       .map(lambda x:
+                                    set_resampling_for_continuous_raster(x,
+                                                                         self.resampling_method,
+                                                                         self.spatial_resolution,
+                                                                         bbox
+                                                                         )
+                            )
+                       .reduce(ee.Reducer.mean())
+                       )
 
         albedo_mean_ic = ee.ImageCollection(albedo_mean)
         data = get_image_collection(

city_metrix/layers/alos_dsm.py

@@ -3,27 +3,38 @@
 import xarray as xr
 
 
-from .layer import Layer, get_image_collection
+from .layer import Layer, get_image_collection, set_resampling_for_continuous_raster
 
 
 class AlosDSM(Layer):
     """
     Attributes:
         spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster)
+        resampling_method: interpolation method used by Google Earth Engine. Default is 'bilinear'. All options are: ('bilinear', 'bicubic', None).
     """
 
-    def __init__(self, spatial_resolution=30, **kwargs):
+    def __init__(self, spatial_resolution:int=30, resampling_method:str='bilinear', **kwargs):
         super().__init__(**kwargs)
         self.spatial_resolution = spatial_resolution
+        self.resampling_method = resampling_method
 
-    def get_data(self, bbox):
+    def get_data(self, bbox: tuple[float, float, float, float]):
         alos_dsm = ee.ImageCollection("JAXA/ALOS/AW3D30/V3_2")
 
-        alos_dsm_ic = ee.ImageCollection(alos_dsm
-                                         .filterBounds(ee.Geometry.BBox(*bbox))
-                                         .select('DSM')
-                                         .mean()
-                                         )
+        alos_dsm_ic = ee.ImageCollection(
+            alos_dsm
+            .filterBounds(ee.Geometry.BBox(*bbox))
+            .select('DSM')
+            .map(lambda x:
+                 set_resampling_for_continuous_raster(x,
+                                                      self.resampling_method,
+                                                      self.spatial_resolution,
+                                                      bbox
+                                                      )
+                 )
+            .mean()
+        )
+
 
         data = get_image_collection(
             alos_dsm_ic,

city_metrix/layers/layer.py

@@ -325,6 +325,29 @@ def get_stats_funcs(stats_func):
         return [stats_func]
 
 
+def set_resampling_for_continuous_raster(image: ee.Image, resampling_method: str, resolution: int,
+                                         bbox: tuple[float, float, float, float]):
+    """
+    Function sets the resampling method on the GEE query dictionary for use on continuous raster layers.
+    GEE only supports bilinear and bicubic interpolation methods.
+    """
+    valid_raster_resampling_methods = ['bilinear', 'bicubic', None]
+
+    if resampling_method not in valid_raster_resampling_methods:
+        raise ValueError(f"Invalid resampling method ('{resampling_method}'). "
+                         f"Valid methods: {valid_raster_resampling_methods}")
+
+    if resampling_method is None:
+        data = image
+    else:
+        crs = get_utm_zone_epsg(bbox)
+        data = (image
+                .resample(resampling_method)
+                .reproject(crs=crs, scale=resolution))
+
+    return data
+
+
 def get_image_collection(
         image_collection: ImageCollection,
         bbox: Tuple[float],

city_metrix/layers/nasa_dem.py

@@ -2,25 +2,34 @@
 import xee
 import xarray as xr
 
-from .layer import Layer, get_image_collection
+from .layer import Layer, get_image_collection, set_resampling_for_continuous_raster
 
 
 class NasaDEM(Layer):
     """
     Attributes:
         spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster)
+        resampling_method: interpolation method used by Google Earth Engine. Default is 'bilinear'. All options are: ('bilinear', 'bicubic', None).
     """
 
-    def __init__(self, spatial_resolution=30, **kwargs):
+    def __init__(self, spatial_resolution:int=30, resampling_method:str='bilinear', **kwargs):
         super().__init__(**kwargs)
         self.spatial_resolution = spatial_resolution
+        self.resampling_method = resampling_method
 
-    def get_data(self, bbox):
+    def get_data(self, bbox: tuple[float, float, float, float]):
         nasa_dem = ee.Image("NASA/NASADEM_HGT/001")
 
         nasa_dem_elev = (ee.ImageCollection(nasa_dem)
                          .filterBounds(ee.Geometry.BBox(*bbox))
                          .select('elevation')
+                         .map(lambda x:
+                              set_resampling_for_continuous_raster(x,
+                                                                   self.resampling_method,
+                                                                   self.spatial_resolution,
+                                                                   bbox
+                                                                   )
+                              )
                          .mean()
                          )
 

tests/resources/bbox_constants.py

@@ -24,3 +24,12 @@
     -38.39993,-12.93239
 )
 
+BBOX_NLD_AMSTERDAM_TEST = (
+    4.9012,52.3720,
+    4.9083,52.3752
+)
+
+BBOX_NLD_AMSTERDAM_LARGE_TEST = (
+    4.884629880473071,52.34146514406914,
+    4.914180290924863,52.359560786247165
+)

tests/resources/layer_dumps_for_br_lauro_de_freitas/conftest.py → tests/resources/layer_dumps/conftest.py

@@ -4,7 +4,8 @@
 import shutil
 from collections import namedtuple
 
-from tests.resources.bbox_constants import BBOX_BRA_LAURO_DE_FREITAS_1
+from tests.resources.bbox_constants import BBOX_BRA_LAURO_DE_FREITAS_1, BBOX_NLD_AMSTERDAM_TEST, \
+    BBOX_NLD_AMSTERDAM_LARGE_TEST
 from tests.tools.general_tools import create_target_folder, is_valid_path
 
 # RUN_DUMPS is the master control for whether the writes and tests are executed
@@ -19,6 +20,8 @@
 # Both the tests and QGIS file are implemented for the same bounding box in Brazil.
 COUNTRY_CODE_FOR_BBOX = 'BRA'
 BBOX = BBOX_BRA_LAURO_DE_FREITAS_1
+# BBOX = BBOX_NLD_AMSTERDAM_TEST
+# BBOX = BBOX_NLD_AMSTERDAM_LARGE_TEST
 
 # Specify None to write to a temporary default folder otherwise specify a valid custom target path.
 CUSTOM_DUMP_DIRECTORY = None

tests/resources/layer_dumps_for_br_lauro_de_freitas/test_write_layers_using_fixed_resolution.py → tests/resources/layer_dumps/test_write_layers_using_fixed_resolution.py

@@ -3,20 +3,24 @@
 import pytest
 
 from city_metrix.layers import *
-from .conftest import RUN_DUMPS, prep_output_path, verify_file_is_populated, get_file_count_in_folder
+from .conftest import RUN_DUMPS, prep_output_path, verify_file_is_populated
 
 TARGET_RESOLUTION = 5
+TARGET_RESAMPLING_METHOD = 'bilinear'
+# TARGET_RESAMPLING_METHOD = None
 
 @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False')
 def test_write_albedo_fixed_res(target_folder, bbox_info, target_spatial_resolution_multiplier):
     file_path = prep_output_path(target_folder, f'albedo_{TARGET_RESOLUTION}m.tif')
-    Albedo(spatial_resolution=TARGET_RESOLUTION).write(bbox_info.bounds, file_path, tile_degrees=None)
+    (Albedo(spatial_resolution=TARGET_RESOLUTION, resampling_method=TARGET_RESAMPLING_METHOD)
+     .write(bbox_info.bounds, file_path, tile_degrees=None))
     assert verify_file_is_populated(file_path)
 
 @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False')
 def test_write_alos_dsm_fixed_res(target_folder, bbox_info, target_spatial_resolution_multiplier):
     file_path = prep_output_path(target_folder, f'alos_dsm_{TARGET_RESOLUTION}m.tif')
-    AlosDSM(spatial_resolution=TARGET_RESOLUTION).write(bbox_info.bounds, file_path, tile_degrees=None)
+    (AlosDSM(spatial_resolution=TARGET_RESOLUTION, resampling_method=TARGET_RESAMPLING_METHOD)
+     .write(bbox_info.bounds, file_path, tile_degrees=None))
     assert verify_file_is_populated(file_path)
 
 @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False')
@@ -64,7 +68,8 @@ def test_write_land_surface_temperature_fixed_res(target_folder, bbox_info, targ
 @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False')
 def test_write_nasa_dem_fixed_res(target_folder, bbox_info, target_spatial_resolution_multiplier):
     file_path = prep_output_path(target_folder, f'nasa_dem_{TARGET_RESOLUTION}m.tif')
-    NasaDEM(spatial_resolution=TARGET_RESOLUTION).write(bbox_info.bounds, file_path, tile_degrees=None)
+    (NasaDEM(spatial_resolution=TARGET_RESOLUTION, resampling_method=TARGET_RESAMPLING_METHOD)
+     .write(bbox_info.bounds, file_path, tile_degrees=None))
     assert verify_file_is_populated(file_path)
 
 @pytest.mark.skipif(RUN_DUMPS == False, reason='Skipping since RUN_DUMPS set to False')

tests/test_layer_metrics.py

@@ -47,8 +47,23 @@ def test_read_image_collection_scale():
         pytest.approx(expected_y_size, rel=EE_IMAGE_DIMENSION_TOLERANCE) == actual_y_size
     )
 
-def test_albedo_metrics():
-    data = Albedo().get_data(BBOX)
+def test_albedo_metrics_default_resampling():
+    # Default resampling_method is bilinear
+    data = Albedo(spatial_resolution=10).get_data(BBOX)
+
+    # Bounding values
+    expected_min_value = _convert_fraction_to_rounded_percent(0.03)
+    expected_max_value = _convert_fraction_to_rounded_percent(0.34)
+    actual_min_value = _convert_fraction_to_rounded_percent(data.values.min())
+    actual_max_value = _convert_fraction_to_rounded_percent(data.values.max())
+
+    # Value range
+    assert expected_min_value == actual_min_value
+    assert expected_max_value == actual_max_value
+
+
+def test_albedo_metrics_no_resampling():
+    data = Albedo(spatial_resolution=10, resampling_method= None).get_data(BBOX)
 
     # Bounding values
     expected_min_value = _convert_fraction_to_rounded_percent(0.03)
@@ -62,7 +77,7 @@ def test_albedo_metrics():
 
 
 def test_alos_dsm_metrics():
-    data = AlosDSM().get_data(BBOX)
+    data = AlosDSM(resampling_method=None).get_data(BBOX)
 
     # Bounding values
     expected_min_value = 16