Cif 316 add bilinear interpolation for dem dsm and albedo layers

README.md

@@ -43,7 +43,7 @@ To run the module,
 
   1. You need access to Google Earth Engine
   2. Install <https://cloud.google.com/sdk/docs/install>
-  3. If you want to use the ERA5 layer, you need to install the [Climate Data Store (CDS) Application Program Interface (API)](https://cds.climate.copernicus.eu/api-how-to)
+  3. If you want to use the ERA5 layer, you need to install the [Climate Data Store (CDS) Application Program Interface (API)](https://cds.climate.copernicus.eu/how-to-api)
 
 ### Interactive development
 

city_metrix/layers/albedo.py

@@ -1,21 +1,25 @@
 import ee
 
-from .layer import Layer, get_image_collection
+from .layer import Layer, get_image_collection, set_resampling_method
+
 
 class Albedo(Layer):
     """
     Attributes:
         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. Albedo 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=10,
+                 resampling_method='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):
@@ -115,7 +119,16 @@ 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())
+
+        if self.resampling_method is not None:
+            albedo_mean = (s2_albedo
+                           .map(lambda x: set_resampling_method(x, self.resampling_method))
+                           .reduce(ee.Reducer.mean())
+                           )
+        else:
+            albedo_mean = (s2_albedo
+                           .reduce(ee.Reducer.mean())
+                           )
 
         albedo_mean_ic = ee.ImageCollection(albedo_mean)
         data = get_image_collection(
@@ -129,3 +142,5 @@ def calc_s2_albedo(image):
             return data.where(data < self.threshold)
 
         return data
+
+

city_metrix/layers/alos_dsm.py

@@ -1,26 +1,38 @@
 import ee
 
-from .layer import Layer, get_image_collection
+from .layer import Layer, get_image_collection, set_resampling_method
 
 
 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. AlosDSM default is 'bilinear'. All options are: ('bilinear', 'bicubic', None).
     """
 
-    def __init__(self, spatial_resolution=30, **kwargs):
+    def __init__(self, spatial_resolution=30, resampling_method='bilinear', **kwargs):
         super().__init__(**kwargs)
         self.spatial_resolution = spatial_resolution
+        self.resampling_method = resampling_method
 
     def get_data(self, bbox):
         alos_dsm = ee.ImageCollection("JAXA/ALOS/AW3D30/V3_2")
 
-        alos_dsm_ic = ee.ImageCollection(alos_dsm
-                                         .filterBounds(ee.Geometry.BBox(*bbox))
-                                         .select('DSM')
-                                         .mean()
-                                         )
+        if self.resampling_method is not None:
+            alos_dsm_ic = ee.ImageCollection(
+                alos_dsm
+                .filterBounds(ee.Geometry.BBox(*bbox))
+                .map(lambda x: set_resampling_method(x, self.resampling_method), )
+                .select('DSM')
+                .mean()
+            )
+        else:
+            alos_dsm_ic = ee.ImageCollection(
+                alos_dsm
+                .filterBounds(ee.Geometry.BBox(*bbox))
+                .select('DSM')
+                .mean()
+            )
 
         data = get_image_collection(
             alos_dsm_ic,

city_metrix/layers/layer.py

@@ -323,6 +323,18 @@ def get_stats_funcs(stats_func):
         return [stats_func]
 
 
+def set_resampling_method(image: ee.Image, resampling_method: str):
+    valid_raster_resampling_methods = ['bilinear', 'bicubic']
+
+    if resampling_method not in valid_raster_resampling_methods:
+        raise ValueError(f"Invalid resampling method ('{resampling_method}'). "
+                         f"Valid methods: {valid_raster_resampling_methods}")
+
+    data = image.resample(resampling_method)
+
+    return data
+
+
 def get_image_collection(
         image_collection: ImageCollection,
         bbox: Tuple[float],

city_metrix/layers/nasa_dem.py

@@ -1,26 +1,36 @@
 import ee
 
-from .layer import Layer, get_image_collection
+from .layer import Layer, get_image_collection, set_resampling_method
 
 
 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. NasaDEM default is 'bilinear'. All options are: ('bilinear', 'bicubic', None).
     """
 
-    def __init__(self, spatial_resolution=30, **kwargs):
+    def __init__(self, spatial_resolution=30, resampling_method='bilinear', **kwargs):
         super().__init__(**kwargs)
         self.spatial_resolution = spatial_resolution
+        self.resampling_method = resampling_method
 
     def get_data(self, bbox):
         nasa_dem = ee.Image("NASA/NASADEM_HGT/001")
 
-        nasa_dem_elev = (ee.ImageCollection(nasa_dem)
-                         .filterBounds(ee.Geometry.BBox(*bbox))
-                         .select('elevation')
-                         .mean()
-                         )
+        if self.resampling_method is not None:
+            nasa_dem_elev = (ee.ImageCollection(nasa_dem)
+                             .filterBounds(ee.Geometry.BBox(*bbox))
+                             .map(lambda x: set_resampling_method(x, self.resampling_method), )
+                             .select('elevation')
+                             .mean()
+                             )
+        else:
+            nasa_dem_elev = (ee.ImageCollection(nasa_dem)
+                             .filterBounds(ee.Geometry.BBox(*bbox))
+                             .select('elevation')
+                             .mean()
+                             )
 
         nasa_dem_elev_ic = ee.ImageCollection(nasa_dem_elev)
         data = get_image_collection(

tests/resources/bbox_constants.py

@@ -24,3 +24,13 @@
     -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,7 +3,7 @@
 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
 

tests/test_layer_metrics.py

@@ -47,8 +47,24 @@ 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.32)
+    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 +78,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