Run matched filters in EMIT product¶
Author: Gonzalo Mateo-GarcΓa
Example from work: > V. RΕ―ΕΎiΔka, G. Mateo-Garcia, L. GΓ³mez-Chova, A. Vaughan, L. Guanter, and A. Markham, Semantic segmentation of methane plumes with hyperspectral machine learning models. Scientific Reports 13, 19999 (2023). DOI: 10.1038/s41598-023-44918-6.
pip install "marshsi[mag1c]"
Step 1: download EMIT image¶
In order to download and process the EMIT image we will use the emit reader in the georeader package. See this tutorial for an example of how to load and plot the data.
EMIT product versions¶
NASA's LP DAAC ended forward processing of the EMIT v001 products on 2026-08-31 (EMIT Version 1 and Existing Version 2 Forward Processing Ended). New acquisitions are published as L1B radiance v002, with the L2A mask as a separate v003 product (EMIT Version 2 Level 1B and Level 2A, and Version 3 Level 2A Data Products Released). Older acquisitions are being reprocessed to v002.
v002 product names drop the orbit and scene numbers (EMIT_L1B_RAD_002_20220827T060753 instead of EMIT_L1B_RAD_001_20220827T060753_2223904_013), and reading them requires georeader 2.4 or later. The marshsi retrievals work on both versions. This notebook runs on v002 by default (EMIT_VERSION below), and the last section repeats the demo on v001 of the same acquisition and compares the two.
import os
import warnings
import numpy as np
import matplotlib.pyplot as plt
# mag1c and georeader call np.divide/np.log with where= and no out=; the result is masked afterwards.
warnings.filterwarnings("ignore", message="'where' used without 'out'")
EMIT_FILENAMES = {
"001": "../tests/data/EMIT_L1B_RAD_001_20220827T060753_2223904_013.nc",
"002": "../tests/data/EMIT_L1B_RAD_002_20220827T060753.nc",
}
EMIT_VERSION = "002"
EMIT_FILENAME = EMIT_FILENAMES[EMIT_VERSION]
NEED_DOWNLOAD = not os.path.exists(EMIT_FILENAME)
print(f"File '{EMIT_FILENAME}' exists: {not NEED_DOWNLOAD}")
import georeader
from georeader.readers import emit, download_utils
from georeader import plot
# NASA's data archive requires creating an account for downloading EMIT files directly.
# Create a user and a token at the NASA Earthdata portal (https://urs.earthdata.nasa.gov/profile)
# and set it as the EARTHDATA_TOKEN environment variable (or a GitHub Actions secret).
earthdata_nasa_account = True
if earthdata_nasa_account:
# Also used by EMITImage to download the OBS file next to the radiance file.
emit.AUTH_METHOD = "token"
emit.TOKEN = os.environ.get("EARTHDATA_TOKEN", "") # NASA Earthdata bearer token
def download_emit(filename):
if earthdata_nasa_account:
link = emit.get_radiance_link(filename)
headers = {"Authorization": f"Bearer {emit.TOKEN}"}
download_utils.download_product(link, headers=headers, verify=True, filename=filename)
elif filename == EMIT_FILENAMES["001"]:
# alternatively for the sake of this demo, we also uploaded the v001 scene (RAD and OBS) on our gdrive
import subprocess
out_dir = os.path.dirname(filename) + "/"
subprocess.run(["gdown", "1jyFejO80Q82qUZ5tRPaDTqsBHefuXYa0", "-O", out_dir], check=True)
subprocess.run(["gdown", "1a9k2YZDXfa4g95KlqCcKfkYcaA-QegYU", "-O", out_dir], check=True)
else:
raise ValueError("Only the v001 scene is on gdrive; download v002 with a NASA Earthdata token")
if NEED_DOWNLOAD:
download_emit(EMIT_FILENAME)
else:
print(f"Skipping download, '{EMIT_FILENAME}' already exists")
# cache_radiance=True keeps the full-spectrum windowed radiance in RAM after
# the first load_raw() call. The four downstream reads (mag1c, AT_MF_total_EMIT,
# load_rgb after to_crs, and plume_vetting.compute_emit) then all serve from
# the in-memory array instead of re-decompressing the ~1.85 GB file each time.
# The cache propagates through to_crs / read_from_bands clones by reference,
# so emit_image_utm built below shares it with rst.
rst = emit.EMITImage(EMIT_FILENAME, cache_radiance=True)
rst
Run mag1c retrieval on the EMIT product¶
Run mag1c filter retrieval based on the work of Foote et al. 2020.
from marshsi.emit import mag1c_emit
mfoutput, albedo = mag1c_emit.mag1c_emit(rst, column_step=4, georeferenced=False)
mfoutput
Plot results on satellite coordinates:
fig, ax = plt.subplots(1,2, figsize=(16,8), tight_layout=True, sharey=True)
im = ax[0].imshow(mfoutput, vmin=0,vmax=1750)
ax[0].set_title(r"$\Delta$CH$_4$ [ppm x m]")
plot.colorbar_next_to(im, ax[0])
im = ax[1].imshow(albedo, vmin=0)
ax[1].set_title("Albedo")
plot.colorbar_next_to(im, ax[1])
Run WMF retrieval on EMIT product¶
Run wide matched filter retrieval based on the work of Roger et al. 2024.
from marshsi.emit import retrieval_upv_emit
mf_classic, mf_extended, mf_combo, mf_extended_filtered, rad = retrieval_upv_emit.AT_MF_total_EMIT(rst, georeferenced=False) # georreferenced=False
Plot results on satellite coordinates:
fig, ax = plt.subplots(2,2, figsize=(20,20), tight_layout=True, sharey=True, sharex=True)
axi = ax[0,0]
im = axi.imshow(mf_classic, vmin=0,vmax=0.3)
axi.set_title(r"MF $\Delta$CH$_4$ [ppm]")
plot.colorbar_next_to(im, axi)
axi = ax[0,1]
im = axi.imshow(mf_extended, vmin=0,vmax=0.3)
axi.set_title(r"WMF $\Delta$CH$_4$ [ppm]")
plot.colorbar_next_to(im, axi)
axi = ax[1,0]
im = axi.imshow(mf_combo, vmin=0,vmax=0.3)
axi.set_title(r"MF Combo $\Delta$CH$_4$ [ppm]")
plot.colorbar_next_to(im, axi)
# rad
axi = ax[1,1]
im = axi.imshow(rad, vmin=0)
axi.set_title(f"Radiance at {retrieval_upv_emit.RAD_WAVELENGTH} in {rst.units}")
plot.colorbar_next_to(im, axi)
Georreference to UTM¶
crs_utm = georeader.get_utm_epsg(rst.footprint("EPSG:4326"))
emit_image_utm = rst.to_crs(crs_utm)
rgb = emit_image_utm.load_rgb(as_reflectance=True)
fig, ax = plt.subplots(1,3, figsize=(15,5), sharey=True)
mag1c_geo = emit_image_utm.georreference(mfoutput, fill_value_default=-1)
wmf_geo = emit_image_utm.georreference(mf_extended, fill_value_default=-1)
plot.show(rgb, ax=ax[0], title= "RGB",mask=True) # add_scalebar=True -> pip install matplotlib_scalebar
plot.show(mag1c_geo, ax=ax[1], title= r"Mag1c $\Delta X$CH$_4$ [ppm x m]",mask=True,vmin=0, vmax=1750,
add_colorbar_next_to=True)
plot.show(wmf_geo, ax=ax[2], title= r"WMF $\Delta X$CH$_4$ [ppm]", mask=True, vmin=0, vmax=1750/8_000,
add_colorbar_next_to=True)
Plume vetting¶
Plume vetting is a spectral validation step that checks whether the pixels inside each detected plume actually show the CHβ absorption signature relative to nearby background pixels. The method is described in:
> C. Xiang, D. R. Thompson, R. O. Green, J. E. Fahlen, A. K. Thorpe, P. G. Brodrick,
> R. W. Coleman, A. M. Lopez, and C. D. Elder, Identification of false methane plumes
> for orbital imaging spectrometers: A case study with EMIT,
> Remote Sensing of Environment, 2025.
>
import geopandas as gpd
import pandas as pd
import numpy as np
from georeader.geotensor import GeoTensor
from marshsi.quantification import ATMOSPHERE_HEIGHT_METHANE
from marshsi import plume_vetting
from importlib import reload
reload(plume_vetting)
PLUMES_FILE = "../tests/data/plume_vetting/plumes_emit.gpkg"
plumes_gdf = gpd.read_file(PLUMES_FILE)
polygons = [geom.geoms[0] for geom in plumes_gdf.geometry]
# Convert WMF from ppm β ppmΒ·m (multiply by effective column height in metres)
wmf_values = wmf_geo.values
if wmf_values.ndim == 3:
wmf_values = wmf_values[0]
wmf_arr = np.where(
wmf_values == wmf_geo.fill_value_default,
wmf_geo.fill_value_default,
wmf_values * ATMOSPHERE_HEIGHT_METHANE,
)
wmf_ppmxm = GeoTensor(wmf_arr, transform=wmf_geo.transform, crs=wmf_geo.crs,
fill_value_default=wmf_geo.fill_value_default)
results = plume_vetting.compute_emit(
emit_image=emit_image_utm,
cmf=wmf_ppmxm,
polygons=polygons,
min_polygon_size=20,
random_seed=42,
)
summary = (
pd.DataFrame.from_dict(results, orient="index")
.reindex(columns=["D_norm", "alpha_con_len"])
.rename_axis("polygon")
.sort_index()
)
summary
if results:
best_idx = min(results, key=lambda k: results[k]["D_norm"])
print(f"Polygon {best_idx}: D_norm={results[best_idx]['D_norm']:.3f}")
plume_vetting.plot_vetting(results[best_idx], cmf_values=wmf_ppmxm.values)
Same demo on EMIT v001¶
The cells below repeat the retrievals and the plume vetting on the other product version of the same acquisition, and compare both versions. The two versions are georeferenced to the same UTM zone but not to the same grid, because the v002 geolocation differs slightly from v001.
OTHER_VERSION = "001" if EMIT_VERSION == "002" else "002"
other_filename = EMIT_FILENAMES[OTHER_VERSION]
if not os.path.exists(other_filename):
download_emit(other_filename)
# Release the first version's radiance so only one cube is held in memory. The cache is shared by rst,
# its clones and IPython's output history, so clear it instead of only deleting names. rad is a view of
# the float64 cube used inside AT_MF_total_EMIT, so it keeps that cube alive too.
rst.clear_radiance_cache()
del rst, emit_image_utm, rad
rst_other = emit.EMITImage(other_filename, cache_radiance=True)
mfoutput_other, _ = mag1c_emit.mag1c_emit(rst_other, column_step=4, georeferenced=False, display_pbar=False)
mf_extended_other = retrieval_upv_emit.AT_MF_total_EMIT(rst_other, georeferenced=False)[1]
emit_image_utm_other = rst_other.to_crs(crs_utm)
rgb_other = emit_image_utm_other.load_rgb(as_reflectance=True)
mag1c_geo_other = emit_image_utm_other.georreference(mfoutput_other, fill_value_default=-1)
wmf_geo_other = emit_image_utm_other.georreference(mf_extended_other, fill_value_default=-1)
rst_other
fig, ax = plt.subplots(2, 3, figsize=(15, 10), sharex=True, sharey=True, tight_layout=True)
versions = {EMIT_VERSION: (rgb, mag1c_geo, wmf_geo), OTHER_VERSION: (rgb_other, mag1c_geo_other, wmf_geo_other)}
for row, (version, (rgb_v, mag1c_v, wmf_v)) in enumerate(sorted(versions.items())):
plot.show(rgb_v, ax=ax[row, 0], title=f"v{version} RGB", mask=True)
plot.show(mag1c_v, ax=ax[row, 1], title=rf"v{version} Mag1c $\Delta X$CH$_4$ [ppm x m]", mask=True,
vmin=0, vmax=1750, add_colorbar_next_to=True)
plot.show(wmf_v, ax=ax[row, 2], title=rf"v{version} WMF $\Delta X$CH$_4$ [ppm]", mask=True,
vmin=0, vmax=1750/8_000, add_colorbar_next_to=True)
wmf_values_other = wmf_geo_other.values
if wmf_values_other.ndim == 3:
wmf_values_other = wmf_values_other[0]
wmf_ppmxm_other = GeoTensor(
np.where(wmf_values_other == wmf_geo_other.fill_value_default, wmf_geo_other.fill_value_default,
wmf_values_other * ATMOSPHERE_HEIGHT_METHANE),
transform=wmf_geo_other.transform, crs=wmf_geo_other.crs, fill_value_default=wmf_geo_other.fill_value_default,
)
results_other = plume_vetting.compute_emit(
emit_image=emit_image_utm_other,
cmf=wmf_ppmxm_other,
polygons=polygons,
min_polygon_size=20,
random_seed=42,
)
summary_other = (
pd.DataFrame.from_dict(results_other, orient="index")
.reindex(columns=["D_norm", "alpha_con_len"])
.rename_axis("polygon")
.sort_index()
)
pd.concat({f"v{EMIT_VERSION}": summary, f"v{OTHER_VERSION}": summary_other}, axis=1)
The marshsi package and tutorials are released under a Creative Commons non-commercial Share-Alike licence. For using the codes in comercial pipelines written consent by the authors must be provided.
If you find this work useful please cite:
@article{ruzicka_starcop_2023,
title = {Semantic segmentation of methane plumes with hyperspectral machine learning models},
volume = {13},
issn = {2045-2322},
url = {https://www.nature.com/articles/s41598-023-44918-6},
doi = {10.1038/s41598-023-44918-6},
number = {1},
journal = {Scientific Reports},
author = {RΕ―ΕΎiΔka, VΓt and Mateo-Garcia, Gonzalo and GΓ³mez-Chova, Luis and Vaughan, Anna, and Guanter, Luis and Markham, Andrew},
month = nov,
year = {2023},
pages = {19999}
}