GOSAT-based Top-down Total and Natural Methane Emissions

Total and natural methane emissions for 2019 summed to a 1° resolution grid. Methane values for both prior to and after inclusion of GOSAT data to the GEOS-Chem global chemistry transport model, version 1.0
Author

Siddharth Chaudhary, Vishal Gaur

Published

September 21, 2023

Approach

  1. Identify available dates and temporal frequency of observations for the given collection using the GHGC API /stac endpoint. The collection processed in this notebook is the gridded methane emissions data product.
  2. Pass the STAC item into the raster API /stac/tilejson.jsonendpoint.
  3. Using folium.plugins.DualMap, we will visualize two tiles (side-by-side), allowing us to compare time points.
  4. After the visualization, we will perform zonal statistics for a given polygon.

About the Data

The NASA Carbon Monitoring System Flux (CMS-Flux) team analyzed remote sensing observations from Japan’s Greenhouse gases Observing SATellite (GOSAT) to produce the global Committee on Earth Observation Satellites (CEOS) CH₄ Emissions data product. They used an analytic Bayesian inversion approach and the GEOS-Chem global chemistry transport model to quantify annual methane (CH₄) emissions and their uncertainties at a spatial resolution of 1° by 1° and then projected these to each country for 2019.

For more information regarding this dataset, please visit the GOSAT-based Top-down Total and Natural Methane Emissions data overview page.

Install the Required Libraries

Required libraries are pre-installed on the GHG Center Hub. If you need to run this notebook elsewhere, please install them with this line in a code cell:

%pip install requests folium rasterstats pystac_client pandas matplotlib –quiet

Querying the STAC API

First, we are going to import the required libraries. Once imported, they allow better executing a query in the GHG Center Spatio Temporal Asset Catalog (STAC) Application Programming Interface (API) where the granules for this collection are stored.

# Import the following libraries
import requests
import folium
import folium.plugins
from folium import Map, TileLayer
from pystac_client import Client
import branca
import pandas as pd
import matplotlib.pyplot as plt
# Provide STAC and RASTER API endpoints
STAC_API_URL = "http://ghg.center/api/stac"
RASTER_API_URL = "https://ghg.center/api/raster"

# Please use the collection name similar to the one used in STAC collection.

# Name of the collection for gosat budget methane. 
collection_name = "gosat-based-ch4budget-yeargrid-v1"
# Fetching the collection from STAC collections using appropriate endpoint.
collection = requests.get(f"{STAC_API_URL}/collections/{collection_name}").json()
collection

Examining the contents of our collection under the temporal variable, we see that the data is available from January 2012 to December 2018. By looking at the dashboard:time density, we observe that the data is available for only one year, i.e. 2019.

def get_item_count(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 = get_item_count(collection_name)
items = requests.get(f"{STAC_API_URL}/collections/{collection_name}/items?limit={number_of_items}").json()["features"]
print(f"Found {len(items)} items")
# Examining the first item in the collection
items[0]

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

Exploring Changes in GOSAT Methane budgets (CH4) Levels Using the Raster API

In this notebook, we will explore the impacts of methane emissions and by examining changes over time in urban regions. We will 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)
items = {item["properties"]["start_datetime"][:10]: item for item in items} 
asset_name = "prior-total"
# Fetching the min and max values for a specific item
rescale_values = {"max":items[list(items.keys())[0]]["assets"][asset_name]["raster:bands"][0]["histogram"]["max"], "min":items[list(items.keys())[0]]["assets"][asset_name]["raster:bands"][0]["histogram"]["min"]}
items.keys()

Now, we will pass the item id, collection name, and rescaling_factor to the Raster API endpoint. We will do this for first January 2019.

color_map = "rainbow" # please select the color ramp from matplotlib library.
january_2019_tile = requests.get(
    f"{RASTER_API_URL}/stac/tilejson.json?collection={items['2019-01-01']['collection']}&item={items['2019-01-01']['id']}"
    f"&assets={asset_name}"
    f"&color_formula=gamma+r+1.05&colormap_name={color_map}"
    f"&rescale={rescale_values['min']},{rescale_values['max']}", 
).json()
january_2019_tile

Visualizing CH₄ Emissions

# Set initial zoom and center of map for CH₄ Layer
# Centre of map [latitude,longitude]
map_ = folium.Map(location=(34, -118), zoom_start=6)

# January 2019
map_layer_2019 = TileLayer(
    tiles=january_2019_tile["tiles"][0],
    attr="GHG",
    opacity=0.7,
)
map_layer_2019.add_to(map_)

# # January 2012
# map_layer_2012 = TileLayer(
#     tiles=january_2012_tile["tiles"][0],
#     attr="GHG",
#     opacity=0.7,
# )
# map_layer_2012.add_to(map_.m2)

# visualising the map
map_

Calculating Zonal Statistics

To perform zonal statistics, first we need to create a polygon. In this use case we are creating a polygon in Texas (USA).

# Texas, USA
texas_aoi = {
    "type": "Feature",
    "properties": {},
    "geometry": {
        "coordinates": [
            [
                # [13.686159004559698, -21.700046934333145],
                # [13.686159004559698, -23.241974326585833],
                # [14.753560168039911, -23.241974326585833],
                # [14.753560168039911, -21.700046934333145],
                # [13.686159004559698, -21.700046934333145],
                [-95, 29],
                [-95, 33],
                [-104, 33],
                [-104,29],
                [-95, 29]
            ]
        ],
        "type": "Polygon",
    },
}
# We will plug in the coordinates for a location inside the the polygon and a zoom level

aoi_map = Map(
    tiles="OpenStreetMap",
    location=[
        30,-100
    ],
    zoom_start=6,
)

folium.GeoJson(texas_aoi, name="Texas, USA").add_to(aoi_map)
aoi_map
# Check total number of items available
items = requests.get(
    f"{STAC_API_URL}/collections/{collection_name}/items?limit=300"
).json()["features"]
print(f"Found {len(items)} items")
# Explore the first item
items[0]
# 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()
    return {
        **result["properties"],
        "datetime": item["properties"]["start_datetime"][:10],
    }

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

%%time
stats = [generate_stats(item, texas_aoi) for item in items]
stats[0]
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["datetime"])
    return df


df = clean_stats(stats)
df.head(5)
print(items[0]["properties"]["start_datetime"][:10])
tile_2016 = requests.get(
    f"{RASTER_API_URL}/stac/tilejson.json?collection={items[0]['collection']}&item={items[0]['id']}"
    f"&assets={asset_name}"
    f"&color_formula=gamma+r+1.05&colormap_name={color_map}"
    f"&rescale={rescale_values['min']},{rescale_values['max']}",
).json()
tile_2016
# Use bbox initial zoom and map
# Set up a map located w/in event bounds

aoi_map_bbox = Map(
    tiles="OpenStreetMap",
    location=[
        30,-100
    ],
    zoom_start=8,
)

map_layer = TileLayer(
    tiles=tile_2016["tiles"][0],
    attr="GHG", opacity = 0.5
)

map_layer.add_to(aoi_map_bbox)

aoi_map_bbox

Summary

In this notebook we have successfully completed the following steps for the STAC collection for the GOSAT-based Top-down Total and Natural Methane Emissions dataset.

  1. Install and import the necessary libraries
  2. Fetch the collection from STAC collections using the appropriate endpoints
  3. Count the number of existing granules within the collection
  4. Map the methane emission levels
  5. Generate zonal statistics for the area of interest (AOI)

If you have any questions regarding this user notebook, please contact us using the feedback form.

Back to top