Merge branch 'main' of https://github.com/wri/cities-cif into CIF-234…

…-veg-water-change

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,23 +1,29 @@
 import ee
+import xarray
+from dask.diagnostics import ProgressBar
+
+from .layer import Layer, get_utm_zone_epsg, get_image_collection, set_resampling_for_continuous_raster
 
-from .layer import Layer, get_image_collection
 
 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. Default is 'bilinear'. All options are: ('bilinear', 'bicubic', None).
         threshold: threshold value for filtering the retrieval
     """
     MAX_CLOUD_PROB = 30
     S2_ALBEDO_EQN = '((B*Bw)+(G*Gw)+(R*Rw)+(NIR*NIRw)+(SWIR1*SWIR1w)+(SWIR2*SWIR2w))'
 
-    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
 
     ## METHODS
@@ -109,7 +115,17 @@ def calc_s2_albedo(image):
         # S2 MOSAIC AND ALBEDO
         dataset = self.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

@@ -1,26 +1,40 @@
 import ee
+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 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/average_net_building_height.py

@@ -1,6 +1,9 @@
+from dask.diagnostics import ProgressBar
+import xarray as xr
+import xee
 import ee
 
-from .layer import Layer, get_image_collection
+from .layer import Layer, get_utm_zone_epsg, get_image_collection
 
 
 class AverageNetBuildingHeight(Layer):

city_metrix/layers/built_up_height.py

@@ -1,6 +1,9 @@
+from dask.diagnostics import ProgressBar
+import xarray as xr
+import xee
 import ee
 
-from .layer import Layer, get_image_collection
+from .layer import Layer, get_utm_zone_epsg, get_image_collection
 
 
 class BuiltUpHeight(Layer):

city_metrix/layers/esa_world_cover.py

@@ -1,7 +1,9 @@
+from dask.diagnostics import ProgressBar
+from enum import Enum
+import xarray as xr
 import ee
 
-from enum import Enum
-from .layer import Layer, get_image_collection
+from .layer import Layer, get_utm_zone_epsg, get_image_collection
 
 
 class EsaWorldCoverClass(Enum):

city_metrix/layers/glad_lulc.py

@@ -1,7 +1,7 @@
 import xarray as xr
 import ee
 
-from .layer import Layer, get_image_collection
+from .layer import Layer, get_utm_zone_epsg, get_image_collection
 
 
 class LandCoverGlad(Layer):

city_metrix/layers/high_land_surface_temperature.py

@@ -1,9 +1,9 @@
-import datetime
-import ee
-
-from shapely.geometry import box
+from .landsat_collection_2 import LandsatCollection2
 from .land_surface_temperature import LandSurfaceTemperature
 from .layer import Layer
+from shapely.geometry import box
+import datetime
+import ee
 
 class HighLandSurfaceTemperature(Layer):
     """

city_metrix/layers/impervious_surface.py

@@ -1,6 +1,9 @@
+from dask.diagnostics import ProgressBar
+import xarray as xr
+import xee
 import ee
 
-from .layer import Layer, get_image_collection
+from .layer import Layer, get_utm_zone_epsg, get_image_collection
 
 
 class ImperviousSurface(Layer):

city_metrix/layers/land_surface_temperature.py

@@ -1,6 +1,8 @@
+from .landsat_collection_2 import LandsatCollection2
+from .layer import Layer, get_utm_zone_epsg, get_image_collection
+from dask.diagnostics import ProgressBar
 import ee
-
-from .layer import Layer, get_image_collection
+import xarray
 
 class LandSurfaceTemperature(Layer):
     """

city_metrix/layers/layer.py

@@ -2,6 +2,8 @@
 from abc import abstractmethod
 from typing import Union, Tuple
 from uuid import uuid4
+# This osgeo import is essential for proper functioning. Do not remove.
+from osgeo import gdal
 
 import ee
 import boto3
@@ -323,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],
@@ -393,4 +418,4 @@ def offset_meters_to_geographic_degrees(decimal_latitude, length_m):
     lon_degree_offset = abs((length_m / (earth_radius_m * math.cos(math.pi*decimal_latitude/180))) * rad)
     lat_degree_offset = abs((length_m / earth_radius_m) * rad)
 
-    return lon_degree_offset, lat_degree_offset
+    return lon_degree_offset, lat_degree_offset

city_metrix/layers/nasa_dem.py

@@ -1,24 +1,35 @@
 import ee
+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()
                          )
 
@@ -29,5 +40,5 @@ def get_data(self, bbox):
             self.spatial_resolution,
             "NASA DEM"
         ).elevation
-        
+
         return data

city_metrix/layers/natural_areas.py

@@ -1,3 +1,4 @@
+import xarray as xr
 from xrspatial.classify import reclassify
 
 from .layer import Layer

city_metrix/layers/smart_surface_lulc.py

@@ -1,7 +1,9 @@
 import xarray as xr
 import numpy as np
 import pandas as pd
+import geopandas as gpd
 from shapely.geometry import CAP_STYLE, JOIN_STYLE
+from shapely.geometry import box
 from exactextract import exact_extract
 from geocube.api.core import make_geocube
 import warnings

city_metrix/layers/tree_canopy_height.py

@@ -1,6 +1,9 @@
-import ee
+from .layer import Layer, get_utm_zone_epsg, get_image_collection
 
-from .layer import Layer, get_image_collection
+from dask.diagnostics import ProgressBar
+import xarray as xr
+import xee
+import ee
 
 class TreeCanopyHeight(Layer):
     """

city_metrix/layers/tree_cover.py

@@ -1,6 +1,10 @@
+from .layer import Layer, get_utm_zone_epsg, get_image_collection
+
+from dask.diagnostics import ProgressBar
+import xarray as xr
+import xee
 import ee
 
-from .layer import Layer, get_image_collection
 
 class TreeCover(Layer):
     """

city_metrix/layers/urban_land_use.py

@@ -1,6 +1,9 @@
+from dask.diagnostics import ProgressBar
+import xarray as xr
+import xee
 import ee
 
-from .layer import Layer, get_image_collection
+from .layer import Layer, get_utm_zone_epsg, get_image_collection
 
 
 class UrbanLandUse(Layer):

city_metrix/layers/world_pop.py

@@ -1,6 +1,9 @@
+from dask.diagnostics import ProgressBar
+import xarray as xr
+import xee
 import ee
 
-from .layer import Layer, get_image_collection
+from .layer import Layer, get_utm_zone_epsg, get_image_collection
 
 
 class WorldPop(Layer):

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