[CIF-165] Add option to save layers when computing metrics

city_metrix/dashboard.py

@@ -0,0 +1,3 @@
+
+def get_data_lake_path(dataset_name, region_name, iso, admin_level):
+    return f"s3://cities-indicators/{dataset_name}/{iso}/{region_name}/{admin_level}/{iso}-{region_name}-{admin_level}-{dataset_name}.tif"

city_metrix/layers/esa_world_cover.py

@@ -25,6 +25,8 @@ class EsaWorldCover(Layer):
     STAC_COLLECTION_ID = "urn:eop:VITO:ESA_WorldCover_10m_2020_AWS_V1"
     STAC_ASSET_ID = "ESA_WORLDCOVER_10M_MAP"
 
+    name = "esa_world_cover_2020"
+
     def __init__(self, land_cover_class=None, **kwargs):
         super().__init__(**kwargs)
         self.land_cover_class = land_cover_class

city_metrix/layers/layer.py

@@ -13,13 +13,17 @@
 import utm
 import shapely.geometry as geometry
 import pandas as pd
+from uuid import uuid4
+import os
+import boto3
 
 
 MAX_TILE_SIZE = 0.5
 
 
 class Layer:
-    def __init__(self, aggregate=None, masks=[]):
+    def __init__(self, path=None, aggregate=None, masks=[]):
+        self.path = path
         self.aggregate = aggregate
         if aggregate is None:
             self.aggregate = self
@@ -41,7 +45,7 @@ def mask(self, *layers):
         :param layers: lis
         :return:
         """
-        return Layer(aggregate=self, masks=self.masks + list(layers))
+        return Layer(path=self.path, aggregate=self, masks=self.masks + list(layers))
 
     def groupby(self, zones, layer=None):
         """
@@ -102,8 +106,17 @@ def _zonal_stats_fishnet(self, stats_func):
 
     def _zonal_stats_tile(self, tile_gdf, stats_func):
         bbox = tile_gdf.total_bounds
+
         aggregate_data = self.aggregate.get_data(bbox)
-        mask_datum = [mask.get_data(bbox) for mask in self.masks]
+        if self.aggregate.path is not None:
+            write_layer(self.aggregate.path, aggregate_data)
+
+        mask_datum = []
+        for mask in self.masks:
+            mask_data = mask.get_data(bbox)
+
+            if mask.path is not None:
+                write_layer(self.mask.path, mask_data)
 
         # align to highest resolution raster, which should be the largest raster
         # since all are clipped to the extent
@@ -240,3 +253,21 @@ def get_image_collection(
 
     return data
 
+
+def write_layer(path, data):
+    if isinstance(data, xr.DataArray):
+        # for rasters, need to write to locally first then copy to cloud storage
+        if path.startswith("s3://"):
+            tmp_path = f"{uuid4()}.tif"
+            data.rio.to_raster(raster_path=tmp_path, driver="COG")
+
+            s3 = boto3.client('s3')
+            s3.upload_file(tmp_path, path.split('/')[2], '/'.join(path.split('/')[3:]))
+
+            os.remove(tmp_path)
+        else:
+            data.rio.to_raster(raster_path=path, driver="COG")
+    elif isinstance(data, gpd.GeoDataFrame):
+        data.to_file(path, driver="GeoJSON")
+    else:
+        raise NotImplementedError("Can only write DataArray or GeoDataFrame")

city_metrix/layers/tree_cover.py

@@ -13,6 +13,8 @@ class TreeCover(Layer):
 
     NO_DATA_VALUE = 255
 
+    name = "tropical_tree_cover"
+
     def __init__(self, min_tree_cover=None, **kwargs):
         super().__init__(**kwargs)
         self.min_tree_cover = min_tree_cover

city_metrix/metrics/built_land_without_tree_cover.py

@@ -3,15 +3,20 @@
 from city_metrix.layers import EsaWorldCoverClass, TreeCover, EsaWorldCover
 
 
-def built_land_without_tree_cover(zones: GeoDataFrame) -> GeoSeries:
+def built_land_without_tree_cover(
+        zones: GeoDataFrame,
+        tree_cover_path: str=None,
+        esa_path: str=None
+) -> GeoSeries:
     """
     Get percentage of built up land (using ESA world cover)
     with no tree cover (>0 WRI tropical tree cover).
     :param zones: GeoDataFrame with geometries to collect zonal stats on
     :return: Pandas Series of percentages
     """
-    built_up_land = EsaWorldCover(land_cover_class=EsaWorldCoverClass.BUILT_UP)
-    tree_cover = TreeCover(min_tree_cover=1)
+
+    built_up_land = EsaWorldCover(land_cover_class=EsaWorldCoverClass.BUILT_UP, path=esa_path)
+    tree_cover = TreeCover(min_tree_cover=1, path=tree_cover_path)
 
     built_land = built_up_land.groupby(zones).count()
     tree_cover_in_built_land = tree_cover.mask(built_up_land).groupby(zones).count()

city_metrix/metrics/mean_tree_cover.py

@@ -3,10 +3,10 @@
 from city_metrix.layers import TreeCover, EsaWorldCoverClass, EsaWorldCover
 
 
-def mean_tree_cover(zones: GeoDataFrame) -> GeoSeries:
+def mean_tree_cover(zones: GeoDataFrame, tree_cover_path: str=None) -> GeoSeries:
     """
     Get mean tree cover (WRI tropical tree cover).
     :param zones: GeoDataFrame with geometries to collect zonal stats on
     :return: Pandas Series of percentages
     """
-    return TreeCover().groupby(zones).mean().divide(100)
+    return TreeCover(path=tree_cover_path).groupby(zones).mean().divide(100)