Run matched filter on EnMAP product¶
Author: Gonzalo Mateo-García
See example of using EnMAP reader in georeader.
The Window-Wide Matched Filter is proposed in: > J. Roger, L. Guanter, J. Gorroño, and I. Irakulis-Loitxate, “Exploiting the entire near-infrared spectral range to improve the detection of methane plumes with high-resolution imaging spectrometers,” Atmospheric Measurement Techniques, vol. 17, no. 4, pp. 1333–1346, Feb. 2024, doi: 10.5194/amt-17-1333-2024.
pip install marshsi adlfs
Download EnMAP image from UNEP IMEO Cloud Storage¶
import os
SAS_TOKEN = ""
AZURE_STORAGE_ACCOUNT=""
CONTAINER = ""
folder_dest = "../tests/data/EnMAP"
tile = "ENMAP01-____L1B-DT0000149931_20250820T075156Z_002_V010502_20250827T172144Z"
folder_dest_tile = os.path.join(folder_dest, tile)
path_metadata_dst = os.path.join(folder_dest_tile, f"{tile}-METADATA.XML")
NEED_DOWNLOAD = not os.path.exists(path_metadata_dst)
print(f"File '{path_metadata_dst}' exists: {not NEED_DOWNLOAD}")
if NEED_DOWNLOAD:
import adlfs
fs = adlfs.AzureBlobFileSystem(account_name=AZURE_STORAGE_ACCOUNT, sas_token=SAS_TOKEN, container=CONTAINER)
PRODUCT_FOLDERS = {
"SPECTRAL_IMAGE_SWIR": "swir",
"QL_QUALITY_CLOUDSHADOW": "ql/cloudshadow",
"QL_PIXELMASK_SWIR": "ql/pixelmask_swir",
"QL_QUALITY_CIRRUS": "ql/cirrus",
"QL_QUALITY_SNOW": "ql/snow",
"QL_PIXELMASK_VNIR": "ql/pixelmask_vnir",
"QL_QUALITY_TESTFLAGS_SWIR": "ql/testflags_swir",
"QL_QUALITY_HAZE": "ql/haze",
"QL_QUALITY_TESTFLAGS_VNIR": "ql/testflags_vnir",
"QL_QUALITY_CLASSES": "ql/classes",
"QL_SWIR": "ql/swir",
"QL_QUALITY_CLOUD": "ql/cloud",
"QL_VNIR": "ql/vnir",
"SPECTRAL_IMAGE_VNIR": "vnir",
}
if NEED_DOWNLOAD:
os.makedirs(folder_dest_tile, exist_ok=True)
# Copy metadata
path_metadata = f"az://{CONTAINER}/EnMAP/metadata/{tile}.xml"
if not fs.exists(path_metadata):
raise FileNotFoundError(f"File {path_metadata} does not exist")
fs.get(path_metadata, path_metadata_dst)
# Download TIF files
for extension, folder in PRODUCT_FOLDERS.items():
path = f"az://{CONTAINER}/EnMAP/{folder}/{tile}.tif"
path_dst = os.path.join(folder_dest_tile, f"{tile}-{extension}.TIF")
if os.path.exists(path_dst):
continue
if not fs.exists(path):
print(f"File {path} not found skipping")
continue
print(f"Downloading file {path}...")
fs.get(path, path_dst)
else:
print(f"Skipping download, '{path_metadata_dst}' already exists")
!tree {folder_dest_tile}
Load image¶
from georeader.readers import enmap
rst = enmap.EnMAP(xml_file=path_metadata_dst,by_folder=False)
rst
Run WW-MF retrieval on EnMAP product¶
Run wide matched filter retrieval based on the work of Roger et al. 2024.
import logging
logger = logging.getLogger(__name__)
logging.basicConfig(level=logging.INFO)
from marshsi.prismaenmap import retrieval_upv_prisma_enmap
cube = retrieval_upv_prisma_enmap.MF_sunglint_combo_enmap(rst,logger=logger)
Reproject to UTM and apply RPCS¶
import georeader
from georeader import read
import numpy as np
crs = georeader.get_utm_epsg(cube.footprint(crs="EPSG:4326").centroid)
cube_crs = read.read_rpcs(
cube.values.astype(np.float32),
rpcs=rst.rpcs_swir,
dst_crs=crs,
resolution_dst_crs=30,
fill_value_default=cube.fill_value_default,
)
cube_crs
retrieval_upv_prisma_enmap.BAND_NAMES_COMBO
rgb = rst.load_rgb(as_reflectance=True,
dst_crs=crs,
apply_rpcs=True,
resolution_dst_crs=30)
rgb
Plot results¶
from georeader import plot
import matplotlib.pyplot as plt
fig, ax = plt.subplots(1,3, figsize=(15,5), sharey=True)
plot.show(rgb, ax=ax[0], mask=True, title="RGB")
mf = cube_crs.isel({"band": 1})
plot.show(mf, ax=ax[1], vmin=0, vmax=.3, title= "MF $\Delta X$CH$_4$ [ppm]", mask=True)
wmf = cube_crs.isel({"band": 2})
plot.show(wmf, ax=ax[2], vmin=0, vmax=.3, title= "WWMF $\Delta X$CH$_4$ [ppm]", mask=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_enmap.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.values
if wmf_values.ndim == 3:
wmf_values = wmf_values[0]
wmf_arr = np.where(
wmf_values == wmf.fill_value_default,
0,
wmf_values * ATMOSPHERE_HEIGHT_METHANE,
)
wmf_arr = np.where(~np.isfinite(wmf_arr), 0, wmf_arr)
wmf_ppmxm = GeoTensor(wmf_arr, transform=wmf.transform, crs=wmf.crs,
fill_value_default=0)
results = plume_vetting.compute_enmap(
enmapi=rst,
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)
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{roger_2024,
AUTHOR = {Roger, J. and Guanter, L. and Gorroño, J. and Irakulis-Loitxate, I.},
TITLE = {Exploiting the entire near-infrared spectral range to improve the detection of methane plumes with high-resolution imaging spectrometers},
JOURNAL = {Atmospheric Measurement Techniques},
VOLUME = {17},
YEAR = {2024},
NUMBER = {4},
PAGES = {1333--1346},
URL = {https://amt.copernicus.org/articles/17/1333/2024/},
DOI = {10.5194/amt-17-1333-2024}
}