Run matched filters in PRISMA product¶
Author: Gonzalo Mateo-GarcĂa
See example of using PRISMA 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 PRISMA image from UNEP IMEO Cloud Storage¶
import os
SAS_TOKEN = ""
AZURE_STORAGE_ACCOUNT=""
CONTAINER = ""
folder_dest = "../tests/data/PRISMA"
tile = "PRS_L1_STD_OFFL_20250316071846_20250316071850_0001"
path_raster_dst = os.path.join(folder_dest, f"{tile}.he5")
NEED_DOWNLOAD = not os.path.exists(path_raster_dst)
print(f"File '{path_raster_dst}' exists: {not NEED_DOWNLOAD}")
if NEED_DOWNLOAD:
import adlfs
fs = adlfs.AzureBlobFileSystem(account_name=AZURE_STORAGE_ACCOUNT, sas_token=SAS_TOKEN, container=CONTAINER)
if NEED_DOWNLOAD:
os.makedirs(folder_dest, exist_ok=True)
path_raster = f"az://{CONTAINER}/PRISMA/raster/{tile}.he5"
if not fs.exists(path_raster):
raise FileNotFoundError(f"File {path_raster} does not exist")
print(f"Downloading {path_raster}...")
fs.get(path_raster, path_raster_dst)
else:
print(f"Skipping download, '{path_raster_dst}' already exists")
!tree {folder_dest}
Load image¶
from georeader.readers import prisma
rst = prisma.PRISMA(path_raster_dst)
rst
Run WW-MF retrieval on PRISMA product¶
Run wide matched filter retrieval based on the work of Roger et al. 2024.
from marshsi.prismaenmap import retrieval_upv_prisma_enmap
cube_crs = retrieval_upv_prisma_enmap.MF_sunglint_combo_prisma(rst)
cube_crs
rgb_local = rst.load_rgb(as_reflectance=True, raw=False)
rgb_local
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_local, 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, add_colorbar_next_to=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,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.plume_vetting import compute_prisma, plot_vetting
PLUMES_FILE = "../tests/data/plume_vetting/plumes_prisma.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,
wmf.fill_value_default,
wmf_values * ATMOSPHERE_HEIGHT_METHANE,
)
wmf_arr = np.where(~np.isfinite(wmf_arr), wmf.fill_value_default, wmf_arr)
wmf_ppmxm = GeoTensor(wmf_arr, transform=wmf.transform, crs=wmf.crs,
fill_value_default=wmf.fill_value_default)
results = compute_prisma(
pi=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}")
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}
}