Interactive Session: Exploring Urban Greenhouse Gas Emissions

About the Data

The Greenhouse gas And Air Pollutants Emissions System (GRA²PES) dataset at the GHG Center is an aggregated, regridded, monthly high-resolution (0.036 x 0.036°) data product with emissions of both greenhouse gases and air pollutants developed in a consistent framework. The dataset contains emissions over the contiguous United States covering major anthropogenic sectors, including energy, industrial fuel combustion and processes, commercial and residential combustion, oil and gas production, on-road and off-road transportation, etc. Carbon dioxide (CO₂) emissions are developed along with those of air pollutants including CO, NOₓ, SO₂, and PM2.5 with consistency in spatial and temporal distributions. Emissions by sectors are reported as column totals in units of metric tons per km² per month. Spatial-temporal surrogates are developed to distribute CO₂ emissions to grid cells to keep consistency between greenhouse gases and air quality species. The current version of GRA²PES is for 2021. Long-term emissions and more greenhouse gas species (e.g., methane) are under development and will be added in the future.

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 Greenhouse gas And Air Pollutants Emissions System (GRA²PES) data
  • How to visualize datasets using folium and perform zonal statistics
  • How to plot time series for Greenhouse gas And Air Pollutants Emissions System (GRA²PES) and analyze the results for a region of interest

Approach

  1. Identify available dates and temporal frequency of observations for the given collection using the GHGC API /stacendpoint. The collection processed in this notebook is the Wetland Methane Emissions, LPJ-wsl Model data product
  2. Pass the STAC item into the raster API /stac/tilejson.json endpoint
  3. Define the spatial region of interest
  4. Using plugins from folium to visualize two tiles (side-by-side), allowing time point comparison
  5. After the visualization, perform zonal statistics for a given polygon
  6. Plot monthly time series for GRA²PES

Setup

Import the required Python libraries by running the next cell.

from pystac_client import Client
import matplotlib.pyplot as plt
import pandas as pd
import stackstac
import os
os.environ['USE_PYGEOS'] = '0'
import datetime
import geopandas as gpd
import cartopy.crs as ccrs
import cartopy.feature as cf
import surfaceAreaGrid
# For reading COG to xarray
import stackstac
# For reading from S3 bucket
import rioxarray
import hvplot.xarray
pd.options.plotting.backend = 'holoviews'
# Provide STAC and RASTER API endpoints
STAC_API_URL = "http://earth.gov/ghgcenter/api/stac"
RASTER_API_URL = "https://earth.gov/ghgcenter/api/raster"
# Fetching the collection from STAC collections using appropriate endpoint
# The 'requests' library allows a HTTP request possible
# Open catalog, print info and available collections
catalog = Client.open(STAC_API_URL)

# print catalog and collections titles
print('>> Title: ',catalog.title)
print('>> Collections:')
for collection in catalog.get_collections():
    print(collection.id)
>> Title:  US GHG Center STAC API
>> Collections:
ct-ch4-monthgrid-v2023
emit-ch4plume-v1
lpjwsl-wetlandch4-monthgrid-v1
vulcan-ffco2-yeargrid-v4
goes-ch4plume-v1
odiac-ffco2-monthgrid-v2023
tm54dvar-ch4flux-monthgrid-v1
gra2pes-ghg-monthgrid-v1
micasa-carbonflux-daygrid-v1
sedac-popdensity-yeargrid5yr-v4.11
vulcan-ffco2-elc-res-yeargrid-v4
gosat-based-ch4budget-yeargrid-v1
tm54dvar-ch4flux-mask-monthgrid-v1
casagfed-carbonflux-monthgrid-v3
micasa-carbonflux-monthgrid-v1
lpjeosim-wetlandch4-daygrid-v1
lpjeosim-wetlandch4-monthgrid-v1
oco2geos-co2-daygrid-v10r
epa-ch4emission-yeargrid-v2express
lpjwsl-wetlandch4-daygrid-v1
oco2-mip-meanco2budget-yeargrid-v1
oco2-mip-co2budget-yeargrid-v1
eccodarwin-co2flux-monthgrid-v5
odiac-ffco2-monthgrid-v2022
# Select GRA2PES collection: browse item_assets, spatial and temporal extent
collection = catalog.get_collection('gra2pes-ghg-monthgrid-v1')
collection
# Get items from this collection, examine temporal extent of each
items = collection.get_items()
for item in items:
    print(item)
<Item id=gra2pes-ghg-monthgrid-v1-202112>
<Item id=gra2pes-ghg-monthgrid-v1-202111>
<Item id=gra2pes-ghg-monthgrid-v1-202110>
<Item id=gra2pes-ghg-monthgrid-v1-202109>
<Item id=gra2pes-ghg-monthgrid-v1-202108>
<Item id=gra2pes-ghg-monthgrid-v1-202107>
<Item id=gra2pes-ghg-monthgrid-v1-202106>
<Item id=gra2pes-ghg-monthgrid-v1-202105>
<Item id=gra2pes-ghg-monthgrid-v1-202104>
<Item id=gra2pes-ghg-monthgrid-v1-202103>
<Item id=gra2pes-ghg-monthgrid-v1-202102>
<Item id=gra2pes-ghg-monthgrid-v1-202101>
# Explore one item (one month!): Check out properties, assets
item = collection.get_item('gra2pes-ghg-monthgrid-v1-202112')
item

Science topic: Compare CO2 season cycles of different cities

Let’s compare CO2 emissions over the course of a year from two cities in different regions - for example, Chicago and Atlanta.

Read data from STAC into xarray Dataset

search = catalog.search(
    collections='gra2pes-ghg-monthgrid-v1',
    datetime=['2021-01-01T00:00:00Z','2021-12-31T00:00:00Z']
)
# Take a look at the items we found
for item in search.item_collection():
    print(item)
<Item id=gra2pes-ghg-monthgrid-v1-202112>
<Item id=gra2pes-ghg-monthgrid-v1-202111>
<Item id=gra2pes-ghg-monthgrid-v1-202110>
<Item id=gra2pes-ghg-monthgrid-v1-202109>
<Item id=gra2pes-ghg-monthgrid-v1-202108>
<Item id=gra2pes-ghg-monthgrid-v1-202107>
<Item id=gra2pes-ghg-monthgrid-v1-202106>
<Item id=gra2pes-ghg-monthgrid-v1-202105>
<Item id=gra2pes-ghg-monthgrid-v1-202104>
<Item id=gra2pes-ghg-monthgrid-v1-202103>
<Item id=gra2pes-ghg-monthgrid-v1-202102>
<Item id=gra2pes-ghg-monthgrid-v1-202101>
# Read data into an xarray Dataset
ds = stackstac.stack(search.item_collection(),epsg=4326).squeeze()
# For now, focus on CO2
# Later, you can change <band> to whichever species you like: co2, co, nox, so2, pm2.5, and repeat the analysis for another gas!
# You will likely have to change the colorbar range to correspond to the cocnentrations of a different gas.
band='co2'
# Set min and max values for colorbar
vmin = 0
vmax = 1000
# Specify colormap
cmap = 'Spectral_r'
# Generate map
ds.sel(band=band)[9,:,:].hvplot(x='x',y='y',cmap=cmap,coastline=True,clim=(vmin,vmax))
/srv/conda/envs/notebook/lib/python3.12/site-packages/cartopy/io/__init__.py:241: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/110m_physical/ne_110m_coastline.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)
# Take a look at the dimensions of the data
print(ds.sel(band=band).shape)
print(ds.sel(band=band).dims)
(12, 947, 2188)
('time', 'y', 'x')

Select areas of interest using GEOJSON shapes and rioxarray

Read in AOIs

GEOJSONs generated via geojson.io. You can create your own and drag + drop them into this directory to look at any AOIs you’d like.

aoi1 = gpd.read_file(f"{os.getenv('HOME')}/shared/agu-2024/activity-2/chicagoland.geojson").to_crs(epsg=4326)
name1 = 'Chicago'
aoi2 = gpd.read_file(f"{os.getenv('HOME')}/shared/agu-2024/activity-2/atlanta.geojson").to_crs(epsg=4326)
name2 = 'Atlanta'
import cartopy.feature as cf
cmap='Spectral_r'    # colormap definition
fig1=plt.figure(figsize=(8,5))
ax1 = fig1.add_subplot(projection=ccrs.PlateCarree())
# Add our data layer
# We can reuse our vmin, vmax, cmap from the previous section
p1 = ds.sel(band=band)[9,:,:].plot(vmin=vmin,vmax=vmax,cmap=cmap,ax=ax1)
# Now let's draw the shapes of our AOIs
# Shape 1 is green
color1 = '#00A676'
# Shape 2 is pink
color2 = '#eb7bc0'
aoi1.plot(ax=ax1,edgecolor=color1,linewidth=1.5,facecolor='none')
aoi2.plot(ax=ax1,edgecolor=color2,linewidth=1.5,facecolor='none')
ax1.set_xlim(-95.4,-78.5)
ax1.set_ylim(32.6,43.04)
ax1.add_feature(cf.STATES,linewidth=0.5)
/srv/conda/envs/notebook/lib/python3.12/site-packages/cartopy/io/__init__.py:241: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/10m_cultural/ne_10m_admin_1_states_provinces_lakes.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)

Clip data field to AOIs using rioxarray

rio.clip has a few options for selection of grid cells relative to the specified geometry. Here we will use the default functionality. See documentation for more info.

clip1 = ds.rio.clip(geometries=aoi1.geometry.values)
clip2 = ds.rio.clip(geometries=aoi2.geometry.values)
# We also need the latitude values reversed to utilize a function correctly in the next step.
clip1 = clip1.reindex(y=clip1.y[::-1])
clip2 = clip2.reindex(y=clip2.y[::-1])

Here’s what our selections look like:

import cartopy.feature as cf
fig1=plt.figure(figsize=(8,5))
ax1 = fig1.add_subplot(projection=ccrs.PlateCarree())
p1 = clip1.sel(band=band)[9,:,:].plot(vmin=vmin,vmax=vmax,cmap=cmap,ax=ax1)
p2 = clip2.sel(band=band)[9,:,:].plot(vmin=vmin,vmax=vmax,cmap=cmap,ax=ax1)
aoi1.plot(ax=ax1,edgecolor=color1,linewidth=1,facecolor='none')
aoi2.plot(ax=ax1,edgecolor=color2, linewidth=1,facecolor='none')
ax1.set_xlim(-95.4,-78.5)
ax1.set_ylim(32.6,43.04)
ax1.add_feature(cf.STATES,linewidth=0.5)

Calculate and apply area weights for lat/lon grid cells

Emissions are in tons km-2 month-1, but our grid cells are in lat/lon space - calculate how many km are in each grid cell, so we can calculate tons month-1 for each AOI.

# This function returns m2, and we need km2. Note the /1000./1000.
weights1 = surfaceAreaGrid.surfaceAreaGridd(lat_centers=clip1.y.values, \
    lon_centers=clip1.x.values,ret_area=True)/1000./1000.
weights2 = surfaceAreaGrid.surfaceAreaGridd(lat_centers=clip2.y.values,\
    lon_centers=clip2.x.values,ret_area=True)/1000./1000.     
# Apply these weights - simple as a multiplication
clip1_weighted = clip1.sel(band=band)*weights1
clip2_weighted = clip2.sel(band=band)*weights2
# Take the mean of emissions over AOI grid cells
clip1_means = clip1_weighted.mean(dim=['x','y'])
clip2_means = clip2_weighted.mean(dim=['x','y'])
# Reformat some times for plotting purposes
times = [sd[0:10] for sd in clip1_means.start_datetime.values]
times = [datetime.datetime.strptime(sd,'%Y-%m-%d').strftime('%b') for sd in times]

Plot emissions time series for 2021

fig,ax1 = plt.subplots()
ax1.plot(times[::-1],clip1_means.values[::-1],marker='o',color=color1,label=name1,lw=2.5)
ax1.plot(times[::-1],clip2_means.values[::-1],marker='o',color=color2,label=name2,lw=2.5)
#ax2.set_ylabel('Urban',color=urban_color,fontweight='bold')
plt.legend()
plt.title(f'Mean urban {band} emissions from GRA2PES, 2021')
Text(0.5, 1.0, 'Mean urban co2 emissions from GRA2PES, 2021')

Apart from their magnitudes, what other differences do you notice between the CO2 seasonal cycles for these two cities? What might explain these differences? Do think these difference are consistent across other cities, or other GHGs?

Follow-up activity: Do this same comparison between an urban area and a rural area and examine those differences!