Interactive Session: Human Anthropogenic Emissions (ODIAC)

About the Data

The Open-Data Inventory for Anthropogenic Carbon dioxide (ODIAC) is a high-spatial resolution global emission data product of CO₂ emissions from fossil fuel combustion (Oda and Maksyutov, 2011). ODIAC pioneered the combined use of space-based nighttime light data and individual power plant emission/location profiles to estimate the global spatial extent of fossil fuel CO₂ emissions. With the innovative emission modeling approach, ODIAC achieved the fine picture of global fossil fuel CO₂ emissions at a 1x1km.

Requirements

Set up Python Environment - See setup_instructions.md in the /setup/ folder

Learning Objectives

How to use U.S. GHG Center STAC Catalogto access ‘ODIAC Fossil Fuel CO₂ Emissions’ and OCO-2 GEOS Column CO₂ Concentrations data.
How to use VEDA STAC Catalog to access NO₂ Emissions data.
How to visualize datasets using folium and perfom zonal statistics over South Africa.
How to plot time series plot for ODIAC and analyze the results.

Approach

Identify available dates and temporal frequency of observations for the given collection using the GHGC API /stac endpoint. Collection processed in this notebook is ODIAC CO₂ emissions version 2022.
Pass the STAC item into raster API /stac/tilejson.json endpoint
We’ll visualize two tiles (side-by-side) allowing for comparison of each of the time points using folium.plugins.DualMap
After the visualization, we’ll perform zonal statistics for a given polygon.

Setup

Import the required Python libraries.

# Import required libraries

import earthaccess
import warnings
import requests
import pandas as pd
import os
os.environ['USE_PYGEOS'] = '0'
import geopandas
import folium
import folium.plugins
import seaborn as sns
import glob
import numpy as np
import matplotlib.pyplot as plt
import branca.colormap as cm

from folium import Map, TileLayer 
from branca.element import Figure
from pystac_client import Client 
from pyproj import Geod
from shapely import wkt
from datetime import datetime
from folium.plugins import MousePosition

Querying the STAC API

Search for ODIAC Fossil Fuel Co2 Emissions, OCO-2 and NO2 Emissions dataset

# Provide STAC and RASTER API endpoints
STAC_API_URL = "http://ghg.center/api/stac"
RASTER_API_URL = "https://ghg.center/api/raster"

STAC_API_URL_veda = "https://staging-stac.delta-backend.com"
RASTER_API_URL_veda = "https://staging-raster.delta-backend.com"

#Please use the collection name similar to the one used in STAC collection.
# Name of the collection for ODIAC dataset. 
collection_name_odiac = "odiac-ffco2-monthgrid-v2022"
collection_name_oco = "oco2geos-co2-daygrid-v10r"
collection_name_no2 = "no2-monthly"

# Fetching the collection from STAC collections using appropriate endpoint.
collection_odiac = requests.get(f"{STAC_API_URL}/collections/{collection_name_odiac}").json()
collection_odiac

collection_oco2 = requests.get(f"{STAC_API_URL}/collections/{collection_name_oco}").json()
collection_oco2

collection_no2 = requests.get(f"{STAC_API_URL_veda}/collections/{collection_name_no2}").json()
collection_no2

Examining the contents of our collection under summaries we see that the data is available from January 2000 to December 2021. By looking at the dashboard:time density we observe that the periodic frequency of these observations is monthly.

# Create a function that would search for the number of items in above data collection in the STAC API

def get_item_count(STAC_API_URL, collection_id):
    count = 0
    items_url = f"{STAC_API_URL}/collections/{collection_id}/items"

    while True:
        response = requests.get(items_url)

        if not response.ok:
            print("error getting items")
            exit()

        stac = response.json()
        count += int(stac["context"].get("returned", 0))
        next = [link for link in stac["links"] if link["rel"] == "next"]

        if not next:
            break
        items_url = next[0]["href"]

    return count

# Check total number of items available
number_of_items_odiac = get_item_count(STAC_API_URL,collection_name_odiac)
items_odiac = requests.get(f"{STAC_API_URL}/collections/{collection_name_odiac}/items?limit={number_of_items_odiac}").json()["features"]
print(f"Found {len(items_odiac)} odiac items")

number_of_items_oco2 = get_item_count(STAC_API_URL,collection_name_oco)
items_oco2 = requests.get(f"{STAC_API_URL}/collections/{collection_name_oco}/items?limit={number_of_items_oco2}").json()["features"]
print(f"Found {len(items_oco2)} oco2 items")

number_of_items_no2 = get_item_count(STAC_API_URL_veda,collection_name_no2)
items_no2 = requests.get(f"{STAC_API_URL_veda}/collections/{collection_name_no2}/items?limit={number_of_items_no2}").json()["features"]
print(f"Found {len(items_no2)} no2 items")

This makes sense as there are 22 years between 2000 - 2021, with 12 months per year, meaning 264 records in total.

Below, we are entering the minimum and maximum values to provide our upper and lower bounds in rescale_values.

Exploring Changes in Carbon Dioxide (CO₂) levels using the Raster API

We will explore changes in fossil fuel emissions in urban egions. In this notebook, we’ll explore the impacts of these emissions and explore these changes over time. We’ll then visualize the outputs on a map using folium.

# to access the year value from each item more easily, this will let us query more explicity by year and month (e.g., 2020-02)
# Set the asset values to the appropriate parameter

items_odiac = {item["properties"]["start_datetime"][:7]: item for item in items_odiac} 
asset_name = "co2-emissions"

items_oco2 = {item["properties"]["datetime"][:10]: item for item in items_oco2} 
asset_name1 = "xco2"

items_no2 = {item["properties"]["start_datetime"][:10]: item for item in items_no2}

# Below, we are entering the minimum and maximum values to provide our upper and lower bounds 
rescale_values_odiac = {"max":items_odiac[list(items_odiac.keys())[0]]["assets"][asset_name]["raster:bands"][0]["histogram"]["max"], "min":items_odiac[list(items_odiac.keys())[0]]["assets"][asset_name]["raster:bands"][0]["histogram"]["min"]}
rescale_values_oco2 = {'max':415 , 'min': 412}
rescale_values_no2 = {'max': 9050373673124971, 'min': 0}

Now we will pass the item id, collection name, and rescaling_factor to the Raster API endpoint. We will do this twice, once for January 2020 and again for January 2000, so that we can visualize each event independently.

Opening and Exploring Data Using the Raster API

In this notebook, we will explore the temporal impacts of methane emissions. We will visualize the outputs on a map using ‘folium’.

# Set a colormap for the granule
# Please refer to matplotlib library if you'd prefer choosing a different color ramp (https://matplotlib.org/stable/users/explain/colors/colormaps.html)
color_map = "rainbow" # please select the color ramp from matplotlib library.

# Extract and display the January 2020 tile using the appropriate ID, colormap, rescale values, and datetime (items['2020-01']) 
january_2020_tile = requests.get(
    f"{RASTER_API_URL}/stac/tilejson.json?collection={items_odiac['2020-01']['collection']}&item={items_odiac['2020-01']['id']}"
    f"&assets={asset_name}"
    f"&color_formula=gamma+r+1.05&colormap_name={color_map}"
    f"&rescale={rescale_values_odiac['min']},{rescale_values_odiac['max']}", 
).json()
january_2020_tile

color_map1 = "magma"
# Extract and display the OCO-2 January 2020 tile using the appropriate ID, colormap, rescale values, and datetime (items['2020-01-20']) 

oco2_1 = requests.get(
    f"{RASTER_API_URL}/stac/tilejson.json?collection={items_oco2['2020-01-20']['collection']}&item={items_oco2['2020-01-20']['id']}"
    f"&assets={asset_name1}"
    f"&color_formula=gamma+r+1.05&colormap_name={color_map1}"
    f"&rescale={rescale_values_oco2['min']},{rescale_values_oco2['max']}", 
).json()
oco2_1

# Extract and display the NO-2 January 2020 tile using the appropriate ID, colormap, rescale values, and datetime (items['2020-01-01']) 

no2_1 = requests.get(
    f"{RASTER_API_URL_veda}/stac/tilejson.json?collection={items_no2['2020-01-01']['collection']}&item={items_no2['2020-01-01']['id']}"
    f"&assets=cog_default"
    f"&color_formula=gamma+r+1.05&colormap_name={color_map1}"
    f"&rescale={rescale_values_no2['min']},{rescale_values_no2['max']}", 
).json()
no2_1

Define Spatial Region of Interest

For this example, our spatial region of interest (ROI) will be the a region in South Africa.

sa_aoi = {
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "properties": {},
      "geometry": {
        "coordinates": [
          [
            [
              25.519052777398997,
              -24.8470086420499
            ],
            [
              25.519052777398997,
              -28.145634397543844
            ],
            [
              30.29637465013832,
              -28.145634397543844
            ],
            [
              30.29637465013832,
              -24.8470086420499
            ],
            [
              25.519052777398997,
              -24.8470086420499
            ]
          ]
        ],
        "type": "Polygon"
      }
    }
  ]
}

Visualizing CO₂ emissions



map_ = folium.Map(location=(-25.943840, 29.789560), zoom_start=7)

map_layer_2020 = TileLayer(
    tiles=january_2020_tile["tiles"][0],
    attr="GHG",
    name ="ODIAC",
    opacity=0.5,
)
map_layer_2020_2 = TileLayer(
    tiles=no2_1["tiles"][0],
    attr="GHG",
    name = "NO2",
    opacity=0.5,
)
sbs = folium.plugins.SideBySideLayers(layer_left=map_layer_2020, layer_right=map_layer_2020_2)
map_layer_2020.add_to(map_)
map_layer_2020_2.add_to(map_)

# Load the GeoJSON file for South Africa
folium.GeoJson(sa_aoi, name="South Africa").add_to(map_)
sbs.add_to(map_)
MousePosition().add_to(map_)
# visualising the map
map_


# Set initial zoom and center of map
map_ = folium.Map(location=(-25.943840, 29.789560), zoom_start=7)

# January 2020
map_layer_2020_odiac = TileLayer(
    tiles=january_2020_tile["tiles"][0],
    attr="GHG",
    name ="ODIAC",
    opacity=0.5,
)
map_layer_2020_odiac.add_to(map_)

# OCO-2 2020
map_layer_2020_oco2 = TileLayer(
    tiles=oco2_1["tiles"][0],
    attr="GHG",
    name = "OCO2",
    opacity=0.5,
)
map_layer_2020_oco2.add_to(map_)

# NO-2 2020
map_layer_2020_no2 = TileLayer(
    tiles=no2_1["tiles"][0],
    attr="GHG",
    name = "NO2",
    opacity=0.5,
)
map_layer_2020_no2.add_to(map_)

folium.GeoJson(sa_aoi, name="South Africa").add_to(map_)
folium.LayerControl(collapsed=False,position='bottomleft').add_to(map_)

# visualising the map
map_

Generate Statistics and Time Series Lineplots for CO2 Emission

# Check total number of items available
items = requests.get(
    f"{STAC_API_URL}/collections/{collection_name_odiac}/items?limit={number_of_items}"
).json()["features"]
print(f"Found {len(items)} items")

# Explore one item to see what it contains
items

# the bounding box should be passed to the geojson param as a geojson Feature or FeatureCollection
def generate_stats(item, geojson):
    result = requests.post(
        f"{RASTER_API_URL}/cog/statistics",
        params={"url": item["assets"][asset_name]["href"]},
        json=geojson,
    ).json()['features']
    return {
        **result[0]["properties"],
        "start_datetime": item["properties"]["start_datetime"][:7],
    }

With the function above we can generate the statistics for the AOI.

%%time
stats = [generate_stats(item,sa_aoi) for item in items]

stats[0]

# Manipulating and cleaning the stats output from the previous cell
def clean_stats(stats_json) -> pd.DataFrame:
    df = pd.json_normalize(stats_json)
    df.columns = [col.replace("statistics.b1.", "") for col in df.columns]
    df["date"] = pd.to_datetime(df["start_datetime"])
    return df

df = clean_stats(stats)
df.head(5)

Visualizing the Data as a Time Series

We can now explore the ODIAC fossil fuel emission time series available (January 2000 -December 2021) for the Texas, Dallas area of USA. We can plot the data set using the code below:

import matplotlib.pyplot as plt

fig = plt.figure(figsize=(20, 10))

# Set the plot elements 
plt.plot(
    df["date"],
    df["max"],
    color="red",
    linestyle="-",
    linewidth=0.5,
    label="Max monthly CO₂ emissions",
)

# Set the labels for the X and Y axis and assign a title for the plot 
plt.legend()
plt.xlabel("Years")
plt.ylabel("CO2 emissions gC/m2/d")
plt.title("CO2 emission Values for South Africa (2000-2021)")