Cif 314 add ghg layer and metrics

city_metrix/layers/__init__.py

@@ -26,3 +26,4 @@
 from .glad_lulc import LandCoverHabitatGlad
 from .glad_lulc import LandCoverHabitatChangeGlad
 from .cams import Cams
+from .cams_ghg import CamsGhg

city_metrix/layers/cams_ghg.py

@@ -0,0 +1,80 @@
+import xarray as xr
+import ee
+
+from .layer import Layer
+
+class CamsGhg(Layer):
+    # Returned units: metric tonnes of GHG species or tonnes CO2e
+    GWP = {
+        'co2': 1,  # by definition
+        'ch4': 28,  # https://www.ipcc.ch/site/assets/uploads/2018/02/WG1AR5_Chapter08_FINAL.pdf
+        'chlorinated-hydrocarbons': 400,  # Table A2 https://onlinelibrary.wiley.com/doi/pdf/10.1002/0470865172.app2
+        'n2o': 265  # https://www.ipcc.ch/site/assets/uploads/2018/02/WG1AR5_Chapter08_FINAL.pdf
+    }
+    SUPPORTED_SPECIES = ['co2', 'ch4', 'n2o', 'chlorinated-hydrocarbons']
+    SUPPORTED_YEARS = [2010, 2015, 2020, 2023]
+
+    def __init__(self, species=None, sector='sum', co2e=True, year=2023, **kwargs):
+        super().__init__(**kwargs)
+        if species is not None and not species in self.SUPPORTED_SPECIES:
+            raise Exception(f'Unsupported species: {species}')
+        if not year in self.SUPPORTED_YEARS:
+            raise Exception(f'Unsupported year: {year}')
+        if species is None and co2e==False:
+            raise Exception('If sector is unspecified, all supported species will be summed and co2e must be True.')
+        if species is None and sector != 'sum':
+            raise Exception('If sector is unspecified, sector must be \"sum.\"')
+        if species is not None and sector != 'sum':
+            data_ic = ee.ImageCollection(f'projects/wri-datalab/cams-glob-ant/{species}').filter('year', year)
+            sectors = data_ic.aggregate_array('sector').getInfo()
+            if not sector in sectors:
+                raise Exception(f'Sector \"{sector}\" not available for {species} in {year}.')
+        self.species = species  # None means all, summed
+        self.sector = sector
+        self.co2e = co2e  # Want results in CO2e? If so, multiplies by 100-year GWP.
+        self.year = year  # Currently supported: 2010, 2015, 2020, 2023
+
+    def get_data(self, bbox):
+        bbox = list(bbox)
+
+        if self.species is not None:
+            data_ic = ee.ImageCollection(f'projects/wri-datalab/cams-glob-ant/{self.species}')
+            data_im = data_ic.filter(ee.Filter.eq('year', self.year)).filter(ee.Filter.eq('sector', self.sector)).first()
+            scale = data_im.projection().getInfo()['transform'][0]
+            data_im = data_im.multiply(self.GWP[self.species])
+            if self.co2e:
+                data_im = data_im.multiply(self.GWP[self.species])
+            data_im = data_im.multiply(1000000)  # Tg to tonne
+            ic = ee.ImageCollection(data_im)
+            scale = data_im.projection().getInfo()['transform'][0]
+            ds = xr.open_dataset(
+                ic,
+                engine = 'ee',
+                geometry = bbox,
+                crs = 'EPSG:4326',
+                scale = scale
+            )
+        else:  # Sum over all species
+            allrasters_list = []
+            for species in self.SUPPORTED_SPECIES:
+                data_ic = ee.ImageCollection(f'projects/wri-datalab/cams-glob-ant/{species}')
+                data_im = data_ic.filter(ee.Filter.eq('year', self.year)).filter(ee.Filter.eq('sector', 'sum')).first()
+                data_im = data_im.multiply(self.GWP[species])
+                data_im = data_im.multiply(1000000)  # Tg to tonne
+                allrasters_list.append(data_im)
+            allrasters_ic = ee.ImageCollection(allrasters_list)
+            sum_im = allrasters_ic.sum()
+            scale = data_im.projection().getInfo()['transform'][0]
+            ic = ee.ImageCollection(sum_im)
+            ds = xr.open_dataset(
+                ic,
+                engine = 'ee',
+                geometry = bbox,
+                crs = 'EPSG:4326',
+                scale = scale
+            )
+        ds = ds.transpose('time', 'lat', 'lon')
+        ds = ds.squeeze(['time'])
+        ds = ds.rio.set_spatial_dims('lon', 'lat')
+        return ds
+

city_metrix/metrics/ghg_emissions.py

@@ -0,0 +1,38 @@
+import ee
+import geopandas as gpd
+import numpy as np
+from city_metrix.layers import Layer, CamsGhg
+
+SUPPORTED_SPECIES = CamsGhg.SUPPORTED_SPECIES
+SUPPORTED_YEARS = CamsGhg.SUPPORTED_YEARS
+GWP = CamsGhg.GWP
+
+def ghg_emissions(zones, years=[2023]):
+    # supported years: 2010, 2015, 2020, 2023
+    if not isinstance(years, (list, set, tuple)):
+        raise Exception('Parameter years must be list, set, or tuple.')
+    for year in years:
+        if not year in SUPPORTED_YEARS:
+            raise Exception(f'Unsupported year {year}')
+
+    zones_copy = zones.copy(deep=True)
+
+    centroids = [zones.iloc[[i]].centroid for i in range(len(zones))]
+    geoms = [ee.Geometry.Point(c.x[i], c.y[i]) for i, c in enumerate(centroids)]
+    for species in SUPPORTED_SPECIES:
+        data_ic = ee.ImageCollection(f'projects/wri-datalab/cams-glob-ant/{species}')
+        for year in years:
+            year_results = data_ic.filter(ee.Filter.eq('year', year))
+            year_total = gpd.GeoDataFrame([0] * len(zones))
+            sectors = year_results.aggregate_array('sector').getInfo()
+            sectors.sort()
+            for sector in sectors:
+                results = year_results.filter(ee.Filter.eq('sector', sector)).first().reduceRegions(geoms, ee.Reducer.mean())
+                result_Tg = gpd.GeoDataFrame([i['properties']['mean'] for i in results.getInfo()['features']])
+                result_tonne = result_Tg * 1000000
+                zones_copy[f'{species}_{sector}_{year}_tonnes'] = result_tonne
+                zones_copy[f'{species}_{sector}_{year}_tonnes-CO2e'] = (result_tonne * GWP[species])
+                if sector == 'sum':
+                    year_total += (result_tonne * GWP[species])
+            zones_copy[f'all-species_{year}_tonnes-CO2e'] = year_total
+    return zones_copy

tests/test_layers.py

@@ -6,6 +6,7 @@
     AverageNetBuildingHeight,
     BuiltUpHeight,
     Cams,
+    CamsGhg,
     Era5HottestDay,
     EsaWorldCover,
     EsaWorldCoverClass,
@@ -58,6 +59,10 @@ def test_cams():
     data = Cams().get_data(BBOX)
     assert np.size(data) > 0
 
+def test_camsghg():
+    data = CamsGhg().get_data(BBOX)
+    assert np.size(data) > 0
+
 @pytest.mark.skip(reason="CDS API needs personal access token file to run")
 def test_era_5_hottest_day():
     data = Era5HottestDay().get_data(BBOX)