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

…-LAYER-Global-30m-Height-Above-the-Nearest-Drainage

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

@@ -2,25 +2,29 @@
 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):
     """
     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
     """
 
-    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")
 
@@ -34,19 +38,33 @@ def get_data(self, bbox):
         def mask_and_count_clouds(s2wc, geom):
             s2wc = ee.Image(s2wc)
             geom = ee.Geometry(geom.geometry())
-            is_cloud = ee.Image(s2wc.get('cloud_mask')).gt(MAX_CLOUD_PROB).rename('is_cloud')
+            is_cloud = (ee.Image(s2wc.get('cloud_mask'))
+                        .gt(MAX_CLOUD_PROB)
+                        .rename('is_cloud')
+                        )
+
             nb_cloudy_pixels = is_cloud.reduceRegion(
                 reducer=ee.Reducer.sum().unweighted(),
                 geometry=geom,
                 scale=self.spatial_resolution,
                 maxPixels=1e9
             )
-            return s2wc.updateMask(is_cloud.eq(0)).set('nb_cloudy_pixels',
-                                                       nb_cloudy_pixels.getNumber('is_cloud')).divide(10000)
+            mask = (s2wc
+                    .updateMask(is_cloud.eq(0))
+                    .set('nb_cloudy_pixels',nb_cloudy_pixels.getNumber('is_cloud'))
+                    .divide(10000)
+                    )
+
+            return mask
 
         def mask_clouds_and_rescale(im):
-            clouds = ee.Image(im.get('cloud_mask')).select('probability')
-            return im.updateMask(clouds.lt(MAX_CLOUD_PROB)).divide(10000)
+            clouds = ee.Image(im.get('cloud_mask')
+                              ).select('probability')
+            mask = im.updateMask(clouds
+                                 .lt(MAX_CLOUD_PROB)
+                                 ).divide(10000)
+
+            return mask
 
         def get_masked_s2_collection(roi, start, end):
             criteria = (ee.Filter.And(
@@ -55,23 +73,29 @@ def get_masked_s2_collection(roi, start, end):
             ))
             s2 = S2.filter(criteria)  # .select('B2','B3','B4','B8','B11','B12')
             s2c = S2C.filter(criteria)
-            s2_with_clouds = (ee.Join.saveFirst('cloud_mask').apply(**{
-                'primary': ee.ImageCollection(s2),
-                'secondary': ee.ImageCollection(s2c),
-                'condition': ee.Filter.equals(**{'leftField': 'system:index', 'rightField': 'system:index'})
-            }))
+            s2_with_clouds = (
+                ee.Join.saveFirst('cloud_mask')
+                .apply(**{
+                    'primary': ee.ImageCollection(s2),
+                    'secondary': ee.ImageCollection(s2c),
+                    'condition': ee.Filter.equals(**{'leftField': 'system:index', 'rightField': 'system:index'})
+                })
+            )
 
             def _mcc(im):
                 return mask_and_count_clouds(im, roi)
                 # s2_with_clouds=ee.ImageCollection(s2_with_clouds).map(_mcc)
 
             # s2_with_clouds=s2_with_clouds.limit(image_limit,'nb_cloudy_pixels')
-            s2_with_clouds = ee.ImageCollection(s2_with_clouds).map(
-                mask_clouds_and_rescale)  # .limit(image_limit,'CLOUDY_PIXEL_PERCENTAGE')
-            return ee.ImageCollection(s2_with_clouds)
+            s2_with_clouds = (ee.ImageCollection(s2_with_clouds)
+                              .map(mask_clouds_and_rescale)
+                              )  # .limit(image_limit,'CLOUDY_PIXEL_PERCENTAGE')
 
-        # calculate albedo for images
+            s2_with_clouds_ic = ee.ImageCollection(s2_with_clouds)
 
+            return s2_with_clouds_ic
+
+        # calculate albedo for images
         # weights derived from
         # S. Bonafoni and A. Sekertekin, "Albedo Retrieval From Sentinel-2 by New Narrow-to-Broadband Conversion Coefficients," in IEEE Geoscience and Remote Sensing Letters, vol. 17, no. 9, pp. 1618-1622, Sept. 2020, doi: 10.1109/LGRS.2020.2967085.
         def calc_s2_albedo(image):
@@ -89,20 +113,35 @@ def calc_s2_albedo(image):
                 'SWIR1': image.select('B11'),
                 'SWIR2': image.select('B12')
             }
-            return image.expression(S2_ALBEDO_EQN, config).double().rename('albedo')
+
+            albedo = image.expression(S2_ALBEDO_EQN, config).double().rename('albedo')
+
+            return albedo
 
         ## 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())
 
-        data = (get_image_collection(
-            ee.ImageCollection(albedo_mean), bbox, self.spatial_resolution, "albedo")
-                .albedo_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(
+            albedo_mean_ic,
+            bbox,
+            self.spatial_resolution,
+            "albedo"
+        ).albedo_mean
 
         if self.threshold is not None:
             return data.where(data < self.threshold)
 
         return data
-
-

city_metrix/layers/alos_dsm.py

@@ -2,26 +2,45 @@
 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: 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')
+            .map(lambda x:
+                 set_resampling_for_continuous_raster(x,
+                                                      self.resampling_method,
+                                                      self.spatial_resolution,
+                                                      bbox
+                                                      )
+                 )
+            .mean()
+        )
+
 
-    def get_data(self, bbox):
-        dataset = ee.ImageCollection("JAXA/ALOS/AW3D30/V3_2")
-        alos_dsm = ee.ImageCollection(dataset
-                                      .filterBounds(ee.Geometry.BBox(*bbox))
-                                      .select('DSM')
-                                      .mean()
-                                      )
-        data = get_image_collection(alos_dsm, bbox, self.spatial_resolution, "ALOS DSM").DSM
+        data = get_image_collection(
+            alos_dsm_ic,
+            bbox,
+            self.spatial_resolution,
+            "ALOS DSM"
+        ).DSM
 
         return data

city_metrix/layers/average_net_building_height.py

@@ -5,6 +5,7 @@
 
 from .layer import Layer, get_utm_zone_epsg, get_image_collection
 
+
 class AverageNetBuildingHeight(Layer):
     """
     Attributes:
@@ -23,8 +24,13 @@ def get_data(self, bbox):
         # GLOBE - ee.Image("projects/wri-datalab/GHSL/GHS-BUILT-H-ANBH_GLOBE_R2023A")
 
         anbh = ee.Image("projects/wri-datalab/GHSL/GHS-BUILT-H-ANBH_GLOBE_R2023A")
-        data = (get_image_collection(
-            ee.ImageCollection(anbh), bbox, self.spatial_resolution, "average net building height")
-                .b1)
-
+
+        anbh_ic = ee.ImageCollection(anbh)
+        data = get_image_collection(
+            anbh_ic,
+            bbox,
+            self.spatial_resolution,
+            "average net building height"
+        ).b1
+
         return data

city_metrix/layers/built_up_height.py

@@ -22,8 +22,15 @@ def get_data(self, bbox):
         # ANBH is the average height of the built surfaces, USE THIS
         # AGBH is the amount of built cubic meters per surface unit in the cell
         # ee.ImageCollection("projects/wri-datalab/GHSL/GHS-BUILT-H-ANBH_R2023A")
-        
+
         built_height = ee.Image("JRC/GHSL/P2023A/GHS_BUILT_H/2018")
-        data = get_image_collection(ee.ImageCollection(built_height), bbox, self.spatial_resolution, "built up height")
-        return data.built_height
 
+        built_height_ic = ee.ImageCollection(built_height)
+        data = get_image_collection(
+            built_height_ic,
+            bbox,
+            self.spatial_resolution,
+            "built up height"
+        ).built_height
+
+        return data

city_metrix/layers/esa_world_cover.py

@@ -39,19 +39,18 @@ def __init__(self, land_cover_class=None, year=2020, spatial_resolution=10, **kw
 
     def get_data(self, bbox):
         if self.year == 2020:
-            data = get_image_collection(
-                ee.ImageCollection("ESA/WorldCover/v100"),
-                bbox,
-                self.spatial_resolution,
-                "ESA world cover"
-            ).Map
+            esa_data_ic = ee.ImageCollection("ESA/WorldCover/v100")
         elif self.year == 2021:
-            data = get_image_collection(
-                ee.ImageCollection("ESA/WorldCover/v200"),
-                bbox,
-                self.spatial_resolution,
-                "ESA world cover"
-            ).Map
+            esa_data_ic = ee.ImageCollection("ESA/WorldCover/v200")
+        else:
+            raise ValueError(f'Specified year ({self.year}) is not currently supported')
+
+        data = get_image_collection(
+            esa_data_ic,
+            bbox,
+            self.spatial_resolution,
+            "ESA world cover"
+        ).Map
 
         if self.land_cover_class:
             data = data.where(data == self.land_cover_class.value)

city_metrix/layers/glad_lulc.py

@@ -17,8 +17,9 @@ def __init__(self, year=2020, spatial_resolution=30, **kwargs):
         self.spatial_resolution = spatial_resolution
 
     def get_data(self, bbox):
+        lcluc_ic = ee.ImageCollection(ee.Image(f'projects/glad/GLCLU2020/LCLUC_{self.year}'))
         data = get_image_collection(
-            ee.ImageCollection(ee.Image(f'projects/glad/GLCLU2020/LCLUC_{self.year}')),
+            lcluc_ic,
             bbox,
             self.spatial_resolution,
             "GLAD Land Cover"
@@ -80,7 +81,8 @@ def __init__(self, year=2020, spatial_resolution=30, **kwargs):
         self.spatial_resolution = spatial_resolution
 
     def get_data(self, bbox):
-        simplified_glad = LandCoverSimplifiedGlad(year=self.year, spatial_resolution=self.spatial_resolution).get_data(bbox)
+        simplified_glad = (LandCoverSimplifiedGlad(year=self.year, spatial_resolution=self.spatial_resolution)
+                           .get_data(bbox))
         # Copy the original data
         data = simplified_glad.copy(deep=True)
 
@@ -108,8 +110,10 @@ def __init__(self, start_year=2000, end_year=2020, spatial_resolution=30, **kwar
         self.spatial_resolution = spatial_resolution
 
     def get_data(self, bbox):
-        habitat_glad_start = LandCoverHabitatGlad(year=self.start_year, spatial_resolution=self.spatial_resolution).get_data(bbox)
-        habitat_glad_end = LandCoverHabitatGlad(year=self.end_year, spatial_resolution=self.spatial_resolution).get_data(bbox)
+        habitat_glad_start = (LandCoverHabitatGlad(year=self.start_year, spatial_resolution=self.spatial_resolution)
+                              .get_data(bbox))
+        habitat_glad_end = (LandCoverHabitatGlad(year=self.end_year, spatial_resolution=self.spatial_resolution)
+                            .get_data(bbox))
 
         # Class 01: Became habitat between start year and end year
         class_01 = ((habitat_glad_start == 0) & (habitat_glad_end == 1))

city_metrix/layers/high_land_surface_temperature.py

@@ -1,12 +1,10 @@
 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):
     """
     Attributes:
@@ -28,6 +26,7 @@ def get_data(self, bbox):
         end_date = (hottest_date + datetime.timedelta(days=45)).strftime("%Y-%m-%d")
 
         lst = LandSurfaceTemperature(start_date, end_date, self.spatial_resolution).get_data(bbox)
+
         lst_mean = lst.mean(dim=['x', 'y'])
         high_lst = lst.where(lst >= (lst_mean + self.THRESHOLD_ADD))
         return high_lst
@@ -36,12 +35,17 @@ def get_hottest_date(self, bbox):
         centroid = box(*bbox).centroid
 
         dataset = ee.ImageCollection("ECMWF/ERA5/DAILY")
-        AirTemperature = (dataset
-                          .filter(ee.Filter.And(
-            ee.Filter.date(self.start_date, self.end_date),
-            ee.Filter.bounds(ee.Geometry.BBox(*bbox))))
-                          .select(['maximum_2m_air_temperature'], ['tasmax'])
-                          )
+
+        AirTemperature = (
+            dataset
+            .filter(
+                ee.Filter
+                .And(ee.Filter.date(self.start_date, self.end_date),
+                     ee.Filter.bounds(ee.Geometry.BBox(*bbox))
+                     )
+            )
+            .select(['maximum_2m_air_temperature'], ['tasmax'])
+        )
 
         # add date as a band to image collection
         def addDate(image):
@@ -56,8 +60,17 @@ def addDate(image):
 
         # reduce composite to get the hottest date for centroid of ROI
         resolution = dataset.first().projection().nominalScale()
-        hottest_date = str(
-            ee.Number(hottest.reduceRegion(ee.Reducer.firstNonNull(), ee.Geometry.Point([centroid.x, centroid.y]), resolution).get('date')).getInfo())
+        hottest_date = (
+            ee.Number(
+                hottest.reduceRegion(ee.Reducer.firstNonNull(),
+                                     ee.Geometry.Point([centroid.x, centroid.y]),
+                                     resolution
+                                     ).get('date')
+            )
+            .getInfo()
+        )
 
         # convert to date object
-        return datetime.datetime.strptime(hottest_date, "%Y%m%d").date()
+        formated_hottest_data = datetime.datetime.strptime(str(hottest_date), "%Y%m%d").date()
+
+        return formated_hottest_data