Interactive Session: Identifying and quantifying emissions from large point source methane emission events leveraging aircraft and satellite data

About the Data

The National Aeronautics and Space Administration (NASA) has a long history of developing advanced imaging spectrometers for new science and applications, and the Earth Surface Mineral Dust Source Investigation (EMIT) builds on this foundation. EMIT launched to the International Space Station (ISS) on July 14, 2022. The data shows high-confidence research grade methane plumes from point source emitters - updated as they are identified - in keeping with Jet Propulsion Laboratory (JPL) Open Science and Open Data policy.

Methane is a strong greenhouse gas that is invisible to the human eye. Large methane emissions, typically referred to as point source emissions, represent a significant proportion of total methane emissions from the production, transport, and processing of oil and natural gas, landfills, and other sources. By measuring the spectral fingerprint of methane, EMIT can map areas of high methane concentration over background levels in the atmosphere, identifying plume complexes, and estimating the methane enhancements.

Requirements

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

Learning Objectives

How to use U.S. GHG Center STAC Catalog to access EMIT Methane Point Source Plume Complexes data.
How to visualize ISS Spatial coverage and detect plumes along the path of the ISS.
How to visualize datasets using folium.

Approach

Visualize the spatial coverage of the ISS and the observed plumes along its path
Identify available dates and temporal frequency of observations for a given collection using the GHGC API /stac endpoint. The collection processed in this notebook is the Earth Surface Mineral Dust Source Investigation (EMIT) methane emission plumes data product
Pass the STAC item into the raster API /stac/tilejson.json endpoint
Display the plumes using folium.Map

Setup

Import the required Python libraries by running the next cell.

import os

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

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

NASA Earth Data Login Credentials

To download or stream NASA data you will need an Earthdata account, you can create one here https://urs.earthdata.nasa.gov/home. We will use the login function from the earthaccess library for authentication before downloading at the end of the notebook. This function can also be used to create a local .netrc file if it doesn’t exist or add your login info to an existing .netrc file. If no Earthdata Login credentials are found in the .netrc you’ll be prompted for them. This step is not necessary to conduct searches but is needed to download or stream data.

Define the Spatial Region of Interest

For this example, our spatial region of interest (ROI) will be the Permian Basin, USA. The Permian Basin is a regional depression expanding from the eastern side of New Mexico to West Texas, covering roughly 75 thousand square miles of land. This region is known to be “the largest hydrocarbon-producing basins in the United States and the world” U.S. Energy Information Administration. Consequently, the Permian Basin is prone to excessive air pollution and a suitable site for monitoring methane variability and concentrations using the EMIT data.

# Loading Permian Basin shape file as a region of interest
# User can pass any json or shape file here
# The "geopandas" library makes working with geospatial data easier by facilitating spatial operations on geometric types
file_path = f"{os.getenv('HOME')}/shared/data/use-case-3"   
permian_basin = geopandas.read_file(f"{file_path}/Permian.zip")
permian_basin = permian_basin.set_crs(4326,allow_override=True)

# The ISS spatial coverage data is called on here
coverage = geopandas.read_file(f'{file_path}/coverage.json')
coverage = coverage.set_crs(4326,allow_override=True)

# The ISS Date Time is converted to strings
coverage['new_start_time'] = coverage.start_time.astype('str')

# locations of where plumes are detected for the August 1, 2022 - present time period
metadata_json = geopandas.read_file(f'{file_path}/methane_metadata.json')
metadata_json = metadata_json.set_crs(4326,allow_override=True)
metadata_json['new_start_time'] = metadata_json['UTC Time Observed'].astype('str')
metadata_json['Scene FIDs'] = list(map(lambda x : x[0] ,metadata_json['Scene FIDs']))
target_mask = geopandas.read_file(f'{file_path}/target_mask.json')
target_mask = target_mask.set_crs(4326,allow_override=True)


# The ISS spatial coverage data is clipped to the ROI (Permian Basin)
result_coverage = geopandas.clip(coverage, permian_basin)
# result_metadata = geopandas.clip(metadata_json, permian_basin)

new_coverage = result_coverage.reset_index(drop=True)
# result_coverage['date'] = result_coverage[(result_coverage['new_start_time'][:10]<='2023-12-15') & (result_coverage['new_start_time'][:10]>='2023-07-01')]
new_coverage['date'] = new_coverage.start_time.dt.strftime('%Y-%m-%d')
new_coverage = new_coverage[(new_coverage['date']<='2023-12-15') & (new_coverage['date']>='2023-07-01')]
new_coverage.reset_index(drop=True, inplace=True)
new_coverage

Opening and Exploring the Spatial Coverage of the ISS and where plumes were detected

In this cell, we will visualize the EMIT-derived methane plumes over the Permian Basin using the folium library.

# Setting a colormap using the “branca.colormap” library 
colormap = cm.linear.YlGn_09.scale(
    coverage.start_time.min(), coverage.start_time.max()
)


# Initializing a map with a center at [43, -100] and a zoom level of 4. The use of tiles = None has been employed to remove the default basemap.
m_ = folium.Map(location=[43, -100], zoom_start=4, tiles=None)


# The TileLayer library helps in manipulating and displaying layers on a map
# For more information on other available tile layers, please visit https://leaflet-extras.github.io/leaflet-providers/preview/
folium.TileLayer(tiles='https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}.png', name='ESRI World Imagery', attr='Tiles &copy; Esri &mdash; Source: Esri, i-cubed, USDA, USGS, AEX, GeoEye, Getmapping, Aerogrid, IGN, IGP, UPR-EGP, and the GIS User Community',overlay='True').add_to(m_)


# From covergae.json fid, geometry has been selected
# map_layer_coverage_target = folium.GeoJson(target_mask, name = ' Target Mask', marker=folium.Marker(icon=folium.Icon(color='blue', icon='globe')))
map_layer_coverage_target = folium.GeoJson(target_mask, name = ' Target Mask', marker=folium.Circle(radius=3, fill_color='blue', color='blue'))
map_layer_coverage_metadata = folium.GeoJson(metadata_json[['Scene FIDs', 'geometry', 'new_start_time']], name = ' Metadata', marker=folium.Circle(radius=3, fill_color='red', color='red'))


map_layer_coverage_target.add_to(m_)
# map_layer_coverage_target.marker(icon=folium.Icon(color='blue', icon='star'))
map_layer_coverage_metadata.add_to(m_)


# We select “geometry” from the “fid” field in “coverage.json”
coverage[['fid', 'geometry', 'new_start_time']].explore(
    "fid",
    categorical=False,
    tooltip=[
        "fid",
        "new_start_time",
    ],
    popup=True,
    style_kwds=dict(fillOpacity=0.1, width=2),
    name="Coverage",
    m=m_,
    legend=False,
    show=False
)


# The “LayerControl” creates the toggle button for all the layers
folium.LayerControl(collapsed=False).add_to(m_)
# Display the map
m_

Opening and Exploring the EMIT Methane Point Source Plume Complexes Data

Search for the EMIT Methane Point Source Plume Complexes Data using the Raster API and its STAC collection name

# 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 methane emission plumes 
collection_name = "emit-ch4plume-v1"

# Fetching the collection from STAC collections using appropriate endpoint
# the 'requests' library allows a HTTP request possible
collection = requests.get(f"{STAC_API_URL}/collections/{collection_name}").json()
collection = {key: value for key, value in collection.items() if key != 'summaries'}
collection

Here we are examining the contents of our collection under the temporal variable. Please keep in mind that data is available from August 2022 to May 2023. By looking at the dashboard: time density, we can see that observations are collected on a daily-basis with plume emissions emerging from multiple places on the same date.

# Create a function that would search for the number of items in above data collection in the STAC API
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

# Apply the above function and check the total number of items available within the collection
number_of_items = get_item_count(collection_name)
plume_complexes = requests.get(f"{STAC_API_URL}/collections/{collection_name}/items?limit={number_of_items}").json()["features"]
parse_plume_complexes = plume_complexes
print(f"Found {len(plume_complexes)} items")

# Exploring the items in the collection
plume_complex_ids = list(map(lambda d: d.get('id', f"id not found in dictionary"), plume_complexes))
plume_complex_ids[:10]

# To access the year value from each item more easily, this will let us query more explicitly by year and month (e.g., 2020-02)
plume_complexes = {plume_complex["id"]: plume_complex for plume_complex in plume_complexes} 


# Set the asset value to the appropriate parameter 
asset_name = "ch4-plume-emissions"

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

# Fetching the min and max values for a specific item
rescale_values = {"max":plume_complexes[list(plume_complexes.keys())[0]]["assets"][asset_name]["raster:bands"][0]["histogram"]["max"], "min":plume_complexes[list(plume_complexes.keys())[0]]["assets"][asset_name]["raster:bands"][0]["histogram"]["min"]}

Now we will pass the item id, collection name, and rescaling_factor to the Raster API endpoint. We will do this for only one item so that we can visualize the event.

# 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 = "magma" 


# Select the item ID which you want to visualize. Item ID is in the format yyyymmdd followed by the timestamp. This ID can be extracted from the Cloud Optimized GeoTIFF (OG) name as well.
plume_complex_id_1 = "EMIT_L2B_CH4PLM_001_20230418T200118_000829"


# Extract the tile using the appropriate ID and colormap 
methane_plume_tile_1 = requests.get(
    f"{RASTER_API_URL}/stac/tilejson.json?collection={plume_complexes[plume_complex_id_1]['collection']}&item={plume_complexes[plume_complex_id_1]['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()
methane_plume_tile_1

# 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 = "magma" 


# Select the item ID which you want to visualize. Item ID is in the format yyyymmdd followed by the timestamp. This ID can be extracted from the COG name as well.
plume_complex_id_2 = "EMIT_L2B_CH4PLM_001_20230614T193706_000038"


# Extract the tile using the appropriate ID and colormap 
methane_plume_tile_2 = requests.get(
    f"{RASTER_API_URL}/stac/tilejson.json?collection={plume_complexes[plume_complex_id_2]['collection']}&item={plume_complexes[plume_complex_id_2]['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()
methane_plume_tile_2

Visualizing CH₄ Emission Plume

# 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)
colormap = "magma" 


#Defining the breaks in the colormap 
color_map = cm.LinearColormap(colors = ['#310597', '#4C02A1', '#6600A7', '#7E03A8', '#9511A1', '#AA2395', '#BC3587', '#CC4778', '#DA5A6A', '#E66C5C', '#F0804E', '#F89540','#FDAC33', '#FDC527', '#F8DF25'], vmin = 0, vmax = 1500 )


# Add an appropriate caption, in this case it would be Parts per million meter
color_map.caption = 'ppm-m'


# Select the item ID which you want to visualize. Item ID is in the format yyyymmdd followed by the timestamp. This ID can be extracted from the COG name as well.
# Set initial zoom and center of map for plume Layer
map_ = folium.Map(location=(methane_plume_tile_1["center"][1], methane_plume_tile_1["center"][0]), zoom_start=13, tiles=None, tooltip = 'test tool tip')
folium.TileLayer(tiles='https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}.png', name='ESRI World Imagery', attr='Tiles &copy; Esri &mdash; Source: Esri, i-cubed, USDA, USGS, AEX, GeoEye, Getmapping, Aerogrid, IGN, IGP, UPR-EGP, and the GIS User Community',overlay='True').add_to(map_)


# Use the 'TileLayer' library to display the raster layer, add an appropriate caption, and adjust the transparency of the layer on the map
map_layer = TileLayer(
    tiles=methane_plume_tile_1["tiles"][0],
    name='Plume Complex Landfill',
    overlay='True',
    attr="GHG",
    opacity=1,
)
map_layer.add_to(map_)


# Adjust map elements 
folium.LayerControl(collapsed=False, position='bottomleft').add_to(map_)
map_.add_child(color_map)
svg_style = '<style>svg#legend {font-size: 14px; background-color: white;}</style>'
map_.get_root().header.add_child(folium.Element(svg_style))


# Visualizing the map
map_

colormap = "magma" # please refer to matplotlib library if you'd prefer choosing a different color ramp.
# For more information on Colormaps in Matplotlib, please visit https://matplotlib.org/stable/users/explain/colors/colormaps.html

#Defining the breaks in the colormap 
color_map = cm.LinearColormap(colors = ['#310597', '#4C02A1', '#6600A7', '#7E03A8', '#9511A1', '#AA2395', '#BC3587', '#CC4778', '#DA5A6A', '#E66C5C', '#F0804E', '#F89540','#FDAC33', '#FDC527', '#F8DF25'], vmin = 0, vmax = 1500 )
color_map.caption = 'ppm-m'

# Set initial zoom and center of map for plume Layer
map_ = folium.Map(location=(methane_plume_tile_2["center"][1], methane_plume_tile_2["center"][0]), zoom_start=13, tiles=None, tooltip = 'test tool tip')
folium.TileLayer(tiles='https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}.png', name='ESRI World Imagery', attr='Tiles &copy; Esri &mdash; Source: Esri, i-cubed, USDA, USGS, AEX, GeoEye, Getmapping, Aerogrid, IGN, IGP, UPR-EGP, and the GIS User Community',overlay='True').add_to(map_)


map_layer = TileLayer(
    tiles=methane_plume_tile_2["tiles"][0],
    name='Plume Complex Agriculture',
    overlay='True',
    attr="GHG",
    opacity=1,
)
map_layer.add_to(map_)


# Adjust map elements 
folium.LayerControl(collapsed=False, position='bottomleft').add_to(map_)
map_.add_child(color_map)
svg_style = '<style>svg#legend {font-size: 14px; background-color: white;}</style>'
map_.get_root().header.add_child(folium.Element(svg_style))


# Visualizing the map
map_

# Creating a dataframe for each plume in the STAC and its corresponding area 
plume_df = pd.DataFrame()
for plume in parse_plume_complexes:
    temp_polygon = geopandas.GeoDataFrame.from_features([plume])['geometry'].values[0]
    geod = Geod(ellps='WGS84')
    ply = wkt.loads(str(temp_polygon))

    temp_dict = {'id':plume['id'], 'geometry':temp_polygon, 'area(km2)':(abs(geod.geometry_area_perimeter(ply)[0])/1e+6), 'date':plume['properties']['datetime'][:10]}
    plume_df = plume_df._append(temp_dict, ignore_index = True)

# Print the generated dataframe
plume_df
# plume['properties']['datetime'][:10]

Visualizing CH4 Emission Plumes in the Permian Basin

# Converting plume_df into a geo dataframe for further geospatial analysis
geo_temp_df = geopandas.GeoDataFrame(plume_df, geometry=plume_df.geometry, crs = 'EPSG:4326')

# Selecting all the targeted plume areas that are within the Permian Basin region
geo_temp_df = geopandas.clip(geo_temp_df, permian_basin)
# filtered_geo_temp_df = geo_temp_df[(geo_temp_df['date']<='2023-12-15') & (geo_temp_df['date']>='2023-07-01') ]
# filtered_geo_temp_df


colormap = "magma" # please refer to matplotlib library if you'd prefer choosing a different color ramp.
# For more information on Colormaps in Matplotlib, please visit https://matplotlib.org/stable/users/explain/colors/colormaps.html


# Defining the breaks in the colormap
color_map = cm.LinearColormap(colors = ['#310597', '#4C02A1', '#6600A7', '#7E03A8', '#9511A1', '#AA2395', '#BC3587', '#CC4778', '#DA5A6A', '#E66C5C', '#F0804E', '#F89540','#FDAC33', '#FDC527', '#F8DF25'], vmin = 0, vmax = 1500 )
color_map.caption = 'ppm-m'


# Set the initial zoom value and adjust the plume layer to be displayed at the 31.8 degree North and -101.7 degree West
map_ = folium.Map(location=[31.8,-101.7], zoom_start=11, tiles=None)
folium.TileLayer(tiles='https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}.png', name='ESRI World Imagery', attr='Tiles &copy; Esri &mdash; Source: Esri, i-cubed, USDA, USGS, AEX, GeoEye, Getmapping, Aerogrid, IGN, IGP, UPR-EGP, and the GIS User Community',overlay='True').add_to(map_)
new_coverage[['fid', 'geometry', 'new_start_time']].explore(
    "fid",
    categorical=True,
    tooltip=[
        "fid",
        "new_start_time",
    ],
    popup=True,
    style_kwds=dict(fillOpacity=0.1, width=2),
    name="Coverage",
    m=map_,
    legend=False,
    show=True
)

map_layer_permian = folium.GeoJson(permian_basin, name= 'Permian Shape')
map_layer_permian.add_to(map_)


# Extract and display the tile using the appropriate ID, colormap, asset name (parameter), and rescale values 
for tile_id in geo_temp_df['id']: 
    methane_plume_tile = requests.get(
        f"{RASTER_API_URL}/stac/tilejson.json?collection={plume_complexes[tile_id]['collection']}&item={plume_complexes[tile_id]['id']}"
        f"&assets={asset_name}"
        f"&color_formula=gamma+r+1.05&colormap_name={colormap}"
        f"&rescale={rescale_values['min']},{rescale_values['max']}", 
    ).json()
    methane_plume_tile

    map_layer = TileLayer(
        tiles=methane_plume_tile["tiles"][0],
        name=tile_id,
        show=True,
        overlay='True',
        attr="GHG",
        opacity=1,
    )
    map_layer.add_to(map_)


# Adjust map elements
folium.LayerControl(collapsed=True, position='bottomleft').add_to(map_)
map_.add_child(color_map)
svg_style = '<style>svg#legend {font-size: 14px; background-color: white;}</style>'
map_.get_root().header.add_child(folium.Element(svg_style))


# Visualizing the map
map_