Merge branch 'main' into feature/transit_stops_osm

.github/requirements.txt

@@ -6,7 +6,7 @@ odc-stac==0.3.8
 pystac-client==0.7.5
 pytest==7.4.3
 xarray-spatial==0.3.7
-xee==0.0.3
+xee==0.0.15
 utm==0.7.0
 osmnx==1.9.3
 dask[complete]==2023.11.0
@@ -16,5 +16,6 @@ geemap==0.32.0
 pip==23.3.1
 boto3==1.34.124
 scikit-learn==1.5.1
+scikit-image==0.24.0
 overturemaps==0.6.0
 exactextract==0.2.0.dev252

city_metrix/layers/__init__.py

@@ -19,3 +19,4 @@
 from .alos_dsm import AlosDSM
 from .overture_buildings import OvertureBuildings
 from .nasa_dem import NasaDEM
+from .impervious_surface import ImperviousSurface

city_metrix/layers/albedo.py

@@ -4,12 +4,20 @@
 
 from .layer import Layer, get_utm_zone_epsg, get_image_collection
 
-
 class Albedo(Layer):
-    def __init__(self, start_date="2021-01-01", end_date="2022-01-01", threshold=None, **kwargs):
+    """
+    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)
+        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):
         super().__init__(**kwargs)
         self.start_date = start_date
         self.end_date = end_date
+        self.spatial_resolution = spatial_resolution
         self.threshold = threshold
 
     def get_data(self, bbox):
@@ -30,7 +38,7 @@ def mask_and_count_clouds(s2wc, geom):
             nb_cloudy_pixels = is_cloud.reduceRegion(
                 reducer=ee.Reducer.sum().unweighted(),
                 geometry=geom,
-                scale=10,
+                scale=self.spatial_resolution,
                 maxPixels=1e9
             )
             return s2wc.updateMask(is_cloud.eq(0)).set('nb_cloudy_pixels',
@@ -88,7 +96,9 @@ def calc_s2_albedo(image):
         s2_albedo = dataset.map(calc_s2_albedo)
         albedo_mean = s2_albedo.reduce(ee.Reducer.mean())
 
-        data = get_image_collection(ee.ImageCollection(albedo_mean), bbox, 10, "albedo").albedo_mean
+        data = (get_image_collection(
+            ee.ImageCollection(albedo_mean), bbox, self.spatial_resolution, "albedo")
+                .albedo_mean)
 
         if self.threshold is not None:
             return data.where(data < self.threshold)

city_metrix/layers/alos_dsm.py

@@ -6,8 +6,14 @@
 
 
 class AlosDSM(Layer):
-    def __init__(self, **kwargs):
+    """
+    Attributes:
+        spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster)
+    """
+
+    def __init__(self, spatial_resolution=30, **kwargs):
         super().__init__(**kwargs)
+        self.spatial_resolution = spatial_resolution
 
     def get_data(self, bbox):
         dataset = ee.ImageCollection("JAXA/ALOS/AW3D30/V3_2")
@@ -16,6 +22,6 @@ def get_data(self, bbox):
                                       .select('DSM')
                                       .mean()
                                       )
-        data = get_image_collection(alos_dsm, bbox, 30, "ALOS DSM").DSM
+        data = get_image_collection(alos_dsm, bbox, self.spatial_resolution, "ALOS DSM").DSM
 
         return data

city_metrix/layers/average_net_building_height.py

@@ -5,10 +5,15 @@
 
 from .layer import Layer, get_utm_zone_epsg, get_image_collection
 
-
 class AverageNetBuildingHeight(Layer):
-    def __init__(self, **kwargs):
+    """
+    Attributes:
+        spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster)
+    """
+
+    def __init__(self, spatial_resolution=100, **kwargs):
         super().__init__(**kwargs)
+        self.spatial_resolution = spatial_resolution
 
     def get_data(self, bbox):
         # https://ghsl.jrc.ec.europa.eu/ghs_buH2023.php
@@ -18,6 +23,8 @@ 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, 100, "average net building height").b1
+        data = (get_image_collection(
+            ee.ImageCollection(anbh), bbox, self.spatial_resolution, "average net building height")
+                .b1)
 
         return data

city_metrix/layers/built_up_height.py

@@ -7,8 +7,14 @@
 
 
 class BuiltUpHeight(Layer):
-    def __init__(self, **kwargs):
+    """
+    Attributes:
+        spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster)
+    """
+
+    def __init__(self, spatial_resolution=100, **kwargs):
         super().__init__(**kwargs)
+        self.spatial_resolution = spatial_resolution
 
     def get_data(self, bbox):
         # Notes for Heat project:
@@ -18,6 +24,6 @@ def get_data(self, bbox):
         # 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, 100, "built up height")
+        data = get_image_collection(ee.ImageCollection(built_height), bbox, self.spatial_resolution, "built up height")
         return data.built_height
 

city_metrix/layers/esa_world_cover.py

@@ -19,32 +19,39 @@ class EsaWorldCoverClass(Enum):
     MANGROVES = 95
     MOSS_AND_LICHEN = 100
 
-
 class EsaWorldCover(Layer):
+    """
+    Attributes:
+        land_cover_class: Enum value from EsaWorldCoverClass
+        year: year used for data retrieval
+        spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster)
+    """
+
     STAC_CATALOG_URI = "https://services.terrascope.be/stac/"
     STAC_COLLECTION_ID = "urn:eop:VITO:ESA_WorldCover_10m_2020_AWS_V1"
     STAC_ASSET_ID = "ESA_WORLDCOVER_10M_MAP"
 
-    def __init__(self, land_cover_class=None, year=2020, **kwargs):
+    def __init__(self, land_cover_class=None, year=2020, spatial_resolution=10, **kwargs):
         super().__init__(**kwargs)
         self.land_cover_class = land_cover_class
         self.year = year
+        self.spatial_resolution = spatial_resolution
 
     def get_data(self, bbox):
         if self.year == 2020:
             data = get_image_collection(
-                    ee.ImageCollection("ESA/WorldCover/v100"),
-                    bbox,
-                    10,
-                    "ESA world cover"
-                ).Map
+                ee.ImageCollection("ESA/WorldCover/v100"),
+                bbox,
+                self.spatial_resolution,
+                "ESA world cover"
+            ).Map
         elif self.year == 2021:
             data = get_image_collection(
-                    ee.ImageCollection("ESA/WorldCover/v200"),
-                    bbox,
-                    10,
-                    "ESA world cover"
-                ).Map
+                ee.ImageCollection("ESA/WorldCover/v200"),
+                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/high_land_surface_temperature.py

@@ -8,19 +8,26 @@
 
 
 class HighLandSurfaceTemperature(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)
+    """
     THRESHOLD_ADD = 3
 
-    def __init__(self, start_date="2013-01-01", end_date="2023-01-01", **kwargs):
+    def __init__(self, start_date="2013-01-01", end_date="2023-01-01", spatial_resolution=30, **kwargs):
         super().__init__(**kwargs)
         self.start_date = start_date
         self.end_date = end_date
+        self.spatial_resolution = spatial_resolution
 
     def get_data(self, bbox):
         hottest_date = self.get_hottest_date(bbox)
         start_date = (hottest_date - datetime.timedelta(days=45)).strftime("%Y-%m-%d")
         end_date = (hottest_date + datetime.timedelta(days=45)).strftime("%Y-%m-%d")
 
-        lst = LandSurfaceTemperature(start_date, end_date).get_data(bbox)
+        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

city_metrix/layers/impervious_surface.py

@@ -0,0 +1,23 @@
+from dask.diagnostics import ProgressBar
+import xarray as xr
+import xee
+import ee
+
+from .layer import Layer, get_utm_zone_epsg, get_image_collection
+
+
+class ImperviousSurface(Layer):
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+
+    def get_data(self, bbox):
+        # load impervious_surface
+        dataset = ee.ImageCollection(ee.Image("Tsinghua/FROM-GLC/GAIA/v10").gt(0))  # change_year_index is zero if permeable as of 2018
+        imperv_surf = ee.ImageCollection(dataset
+                                         .filterBounds(ee.Geometry.BBox(*bbox))
+                                         .select('change_year_index')
+                                         .sum()
+                                         )
+
+        data = get_image_collection(imperv_surf, bbox, 100, "imperv surf")
+        return data.change_year_index

city_metrix/layers/land_surface_temperature.py

@@ -5,12 +5,19 @@
 import ee
 import xarray
 
-
 class LandSurfaceTemperature(Layer):
-    def __init__(self, start_date="2013-01-01", end_date="2023-01-01", **kwargs):
+    """
+    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)
+    """
+
+    def __init__(self, start_date="2013-01-01", end_date="2023-01-01", spatial_resolution=30, **kwargs):
         super().__init__(**kwargs)
         self.start_date = start_date
         self.end_date = end_date
+        self.spatial_resolution = spatial_resolution
 
     def get_data(self, bbox):
         def cloud_mask(image):
@@ -31,5 +38,5 @@ def apply_scale_factors(image):
             .map(apply_scale_factors) \
             .reduce(ee.Reducer.mean())
 
-        data = get_image_collection(ee.ImageCollection(l8_st), bbox, 30, "LST").ST_B10_mean
+        data = get_image_collection(ee.ImageCollection(l8_st), bbox, self.spatial_resolution, "LST").ST_B10_mean
         return data

city_metrix/layers/landsat_collection_2.py

@@ -29,7 +29,8 @@ def get_data(self, bbox):
             fail_on_error=False,
         )
 
-        # TODO: Determine how to output xarray
-
         qa_lst = lc2.where((lc2.qa_pixel & 24) == 0)
         return qa_lst.drop_vars("qa_pixel")
+
+
+

city_metrix/layers/layer.py

@@ -274,9 +274,12 @@ def get_image_collection(
     :param name: optional name to print while reporting progress
     :return:
     """
+    if scale is None:
+        raise Exception("Spatial_resolution cannot be None.")
 
     crs = get_utm_zone_epsg(bbox)
 
+    # See link regarding bug in crs specification https://github.com/google/Xee/issues/118
     ds = xr.open_dataset(
         image_collection,
         engine='ee',

city_metrix/layers/nasa_dem.py

@@ -6,8 +6,14 @@
 
 
 class NasaDEM(Layer):
-    def __init__(self, **kwargs):
+    """
+    Attributes:
+        spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster)
+    """
+
+    def __init__(self, spatial_resolution=30, **kwargs):
         super().__init__(**kwargs)
+        self.spatial_resolution = spatial_resolution
 
     def get_data(self, bbox):
         dataset = ee.Image("NASA/NASADEM_HGT/001")
@@ -16,6 +22,6 @@ def get_data(self, bbox):
                                       .select('elevation')
                                       .mean()
                                       )
-        data = get_image_collection(nasa_dem, bbox, 30, "NASA DEM").elevation
+        data = get_image_collection(nasa_dem, bbox, self.spatial_resolution, "NASA DEM").elevation
 
         return data

city_metrix/layers/natural_areas.py

@@ -4,13 +4,18 @@
 from .layer import Layer
 from .esa_world_cover import EsaWorldCover, EsaWorldCoverClass
 
-
 class NaturalAreas(Layer):
-    def __init__(self, **kwargs):
+    """
+    Attributes:
+        spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster)
+    """
+
+    def __init__(self, spatial_resolution=10, **kwargs):
         super().__init__(**kwargs)
+        self.spatial_resolution = spatial_resolution
 
     def get_data(self, bbox):
-        esa_world_cover = EsaWorldCover().get_data(bbox)
+        esa_world_cover = EsaWorldCover(spatial_resolution=self.spatial_resolution).get_data(bbox)
         reclass_map = {
             EsaWorldCoverClass.TREE_COVER.value: 1,
             EsaWorldCoverClass.SHRUBLAND.value: 1,

city_metrix/layers/ndvi_sentinel2_gee.py

@@ -4,15 +4,18 @@
 class NdviSentinel2(Layer):
     """"
     NDVI = Sentinel-2 Normalized Difference Vegetation Index
-    param: year: The satellite imaging year.
+    Attributes:
+        year: The satellite imaging year.
+        spatial_resolution: raster resolution in meters (see https://github.com/stac-extensions/raster)
     return: a rioxarray-format DataArray
     Author of associated Jupyter notebook: Ted.Wong@wri.org
     Notebook: https://github.com/wri/cities-cities4forests-indicators/blob/dev-eric/scripts/extract-VegetationCover.ipynb
     Reference: https://en.wikipedia.org/wiki/Normalized_difference_vegetation_index
     """
-    def __init__(self, year=None, **kwargs):
+    def __init__(self, year=None, spatial_resolution=10, **kwargs):
         super().__init__(**kwargs)
         self.year = year
+        self.spatial_resolution = spatial_resolution
 
     def get_data(self, bbox):
         if self.year is None:
@@ -39,8 +42,7 @@ def calculate_ndvi(image):
         ndvi_mosaic = ndvi.qualityMosaic('NDVI')
 
         ic = ee.ImageCollection(ndvi_mosaic)
-        ndvi_data = get_image_collection(ic, bbox, 10, "NDVI")
+        ndvi_data = (get_image_collection(ic, bbox, self.spatial_resolution, "NDVI")
+                     .NDVI)
 
-        xdata = ndvi_data.to_dataarray()
-
-        return xdata
+        return ndvi_data

city_metrix/layers/sentinel_2_level_2.py

@@ -50,6 +50,4 @@ def get_data(self, bbox):
         cloud_masked = s2.where(s2 != 0).where(s2.scl != 3).where(s2.scl != 8).where(s2.scl != 9).where(
             s2.scl != 10)
 
-        # TODO: Determine how to output as an xarray
-
         return cloud_masked.drop_vars("scl")