CH4 enhancement con dati Emit e Enmap

Con il metodo CH4 enhancement si riesce a misurare l'eccesso di metano in una colonna d'aria rispetto al fondo mediante immagini iperspettrali utilizzando gli assorbimenti tra 2200 e 2400 nm. Non si tratta quindi di una misura assoluta ma relativa rispetto ai pixel vicini ed e' utile per il rilevamento di plumes

Una variazione di CH4 crea una variazione di radianza di  L₀·(1 − k·Δc)  dove k viene calcolato tramite un modello radiativo come MODTRAN. La variazione e' minima perche' si tratta dell'ordine di 1% sul valore della radianza

La soglia di rilevabilita' con questo metodo e' posta indicativamente a qualche centinaio di kg/h 

Come metodo di retrieval viene spesso usato Matched Filter

Ci sono un paio di progetti su GitHub che permettono il calcolo del CH4 enhancement

Hygas :  utilizza dati Enmap,Prisma e Tanager-1. Richiede installazione delle LUT da https://hive.utah.edu/records/d47hd-s1z16 e necessita' di una scheda NVidia. Ho usato solo Enmap per i problemi nello Swir di Prisma

Hypergas :  utilizza dati L1 Emit e richiede dati accessori come la direzione del vento tramite ERA5 ed il dem (scaricabile da qui)

Se non si vuole effettuare il calcolo sulla propria macchina Nasa distribuisce dati gia' processati di livello L2B  

Si puo' consultare anche https://us-ghg-center.github.io/ghgc-docs/ che riporta mappe giornaliere dei plume piu' significativi con API per il download dei dati. Interessanti sono i notebook per il trattamento dei dati

Come confronto Emit ha una migliore risoluzione spettrale ed un rapporto segnale rumore migliore di Enmap ma ha una peggiore GSD  (60 m contro 30 m) quindi Enmap vede meglio i plume dal punto di vista spaziale. Inoltre dato che Emit e' un imaging mentre Enmap un push broom i dati di Enmap devono essere processati statisticamente lungo le colonne

Considerazione a margine: i dati iperspettrali sono ancora pochi con alcune zone della Toscana che non risultano essere mai state coperte ed altre in cui e' stata fatta una sola acquisizione rendendo impossibile analisi multi-temporale. La copertura migliore e' quella effettuata da Emit 

 

L'elaborazione di CH4 enhancement crea un grafico a campagna guassiana. I valori centrali sono il background calcolato sull'immagine. I plume sono nella coda all'estrema destra della gaussiana (indicativamente si devono cercare i valori oltre i 3 sigma)

 

Discarica di Legoli Enmap 43.567755,10.806777

Il primo risultato e' stato entusiasmante in quanto HyGas ha evidenziato in modo chiaro una possibile presenza di un plume di metano sopra la discarica attiva di Legoli....ovviamente i problemi non hanno tardato a presentarsi

Il medesimo sito e' stato oggetto di osservazione da parte di una ditta privata https://carbonmapper.org/ che  usando Tanager-1 ha visto il plume con una definizione nettamente superiore. Leggendo pero' i metadati si ha che il plume e' stato registrato in uno solo di sei passaggi

 

Allargando pero' il campo si vede segnale praticamente ovunque sulla mappa anche in zone che non possono interessate da emissioni di metano come campi coltivate, corsi di fiumi etc. con una quantita' di falsi positivi difficilmente gestibili senza conoscenza a priori della verita' a terra

 

per lanciare Hygas ho usato 

python scripts/main.py --satellite enmap --mode scene --vnir ./dati/VNIR.TIF --swir ./dati/SWIR.TIF --metadata ./dati/METADATA.XML --lut ./LUT/dataset_ch4_full.hdf5 --snr-reference ./LUT/snr_reference_columnwise.npz --output ./mugello --enmap-mf-mode full-column

 

mentre per HyperGas il seguente script

#!/usr/bin/env python3

"""Recupero di CH4 o CO2 da una scena EMIT L1B con HyperGas.

Uso:

python go_emit.py "EMIT_L1B_*_001_20230814T100819*nc" --species ch4

python go_emit.py "/dati/emit/EMIT_L1B_*.nc" --species co2 --out scena_co2.nc

Servono i file EMIT L1B della stessa scena: radianza (EMIT_L1B_RAD_...nc)

e osservazione/geometria (EMIT_L1B_OBS_...nc), con lo stesso timestamp.

python go.py "EMIT_L1B_*_001_20250808T083836_2522006_003*nc" --species ch4

"""

import argparse

import glob

import sys

import numpy as np

# compatibilita' con dipendenze che usano ancora np.trapz (NumPy recente)

if not hasattr(np, "trapz"):

np.trapz = np.trapezoid

from hypergas import Hyper

def main():

p = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)

p.add_argument("pattern", help="glob dei file EMIT L1B (RAD + OBS), tra virgolette")

p.add_argument("--species", default="ch4", choices=["ch4", "co2"])

p.add_argument("--out", default=None, help="NetCDF di output (default: emit_<species>.nc)")

args = p.parse_args()

filenames = sorted(glob.glob(args.pattern))

print(f"{len(filenames)} file trovati:")

for f in filenames:

print(" ", f)

if not filenames:

sys.exit("Nessun file trovato: controlla il pattern.")

hyp = Hyper(filenames, reader="emit_l1b") # reader EMIT

hyp.load() # radianza, senza bande del vapore acqueo

print(hyp.scene)

hyp.retrieve(species=args.species) # matched filter

print(hyp.scene[args.species])

out = args.out or f"emit_{args.species}.nc"

base = out[:-3] if out.endswith(".nc") else out

# rgb ha 3 bande e il gas 1: nello stesso NetCDF xarray non riesce ad allinearli

# (AlignmentError sulla dimensione 'bands'), quindi li salvo in file separati.

available = [x["name"] for x in hyp.scene.keys()]

for name in [args.species, "rgb"]:

if name not in available:

print(f"Dataset '{name}' non disponibile, salto")

continue

fn = out if name == args.species else f"{base}_{name}.nc"

hyp.scene.save_datasets(datasets=[name], filename=fn,

header_attrs={"author": "Luca"}, writer="cf")

print("Salvato", fn)

if __name__ == "__main__":

main()

i dati Emit vengono distribuiti in formato nc non direttamente georiferito. Per fare una mappa da mettere su QGis

#!/usr/bin/env python3

"""Converte un NetCDF prodotto da HyperGas/Satpy (emit_ch4.nc, emit_ch4_rgb.nc)

in un GeoTIFF georiferito (EPSG:4326) da aprire in QGIS.

Uso:

python nc_to_geotiff.py emit_ch4.nc

python nc_to_geotiff.py emit_ch4_rgb.nc --var rgb

python nc_to_geotiff.py emit_ch4.nc --res 0.0006 # passo griglia in gradi (~60 m)

Se lon/lat sono 2D (dati in geometria sensore) il dato viene ricampionato

(nearest neighbour) su una griglia regolare lat/lon; se sono 1D viene scritto com'e'.

Richiede: xarray, netCDF4, rasterio, pyresample (gia' installati con Satpy).

"""

import argparse

import sys

import numpy as np

import rasterio

import xarray as xr

from pyresample import geometry, kd_tree

from rasterio.transform import from_origin

LON_NAMES = ("longitude", "lon", "x_lon")

LAT_NAMES = ("latitude", "lat", "y_lat")

def find_coord(ds, names):

for n in names:

if n in ds.variables:

return ds[n]

return None

def pick_var(ds, wanted):

if wanted:

if wanted not in ds.data_vars:

sys.exit(f"Variabile '{wanted}' assente. Disponibili: {list(ds.data_vars)}")

return wanted

skip = {"crs", "longitude", "latitude", "lon", "lat"}

cands = [v for v in ds.data_vars

if v not in skip and "bounds" not in v and ds[v].ndim >= 2]

if not cands:

sys.exit(f"Nessuna variabile 2D trovata: {list(ds.data_vars)}")

return cands[0]

def main():

p = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)

p.add_argument("nc")

p.add_argument("--var", default=None, help="nome variabile (default: la prima 2D)")

p.add_argument("--res", type=float, default=0.0006, help="passo griglia in gradi (default 0.0006 ~ 60 m)")

p.add_argument("--radius", type=float, default=150.0, help="raggio di ricerca in metri (default 150)")

p.add_argument("--out", default=None)

args = p.parse_args()

ds = xr.open_dataset(args.nc)

name = pick_var(ds, args.var)

da = ds[name]

print("Variabile:", name, "dims:", da.dims, "shape:", da.shape)

lon = find_coord(ds, LON_NAMES)

lat = find_coord(ds, LAT_NAMES)

if lon is None or lat is None:

sys.exit("Coordinate longitude/latitude non trovate nel file. "

"Incollami l'output di: python -c \"import xarray as xr; print(xr.open_dataset('%s'))\"" % args.nc)

# porta i dati a (bands, y, x)

spatial = [d for d in da.dims if d in lon.dims or d in lat.dims] if lon.ndim == 2 else None

if da.ndim == 2:

data = da.values[np.newaxis, ...]

elif da.ndim == 3:

band_dim = [d for d in da.dims if d not in (lon.dims if lon.ndim == 2 else ())][0]

data = da.transpose(band_dim, *[d for d in da.dims if d != band_dim]).values

else:

sys.exit(f"Dimensioni non gestite: {da.dims}")

data = data.astype("float32")

out = args.out or args.nc.rsplit(".", 1)[0] + ".tif"

if lon.ndim == 1:

# griglia regolare lon/lat

lons, lats = lon.values, lat.values

res_x, res_y = abs(lons[1] - lons[0]), abs(lats[1] - lats[0])

if lats[0] < lats[-1]:

data = data[:, ::-1, :]

lats = lats[::-1]

transform = from_origin(lons.min() - res_x / 2, lats.max() + res_y / 2, res_x, res_y)

arr = data

else:

lons2, lats2 = lon.values, lat.values

lon_min, lon_max = np.nanmin(lons2), np.nanmax(lons2)

lat_min, lat_max = np.nanmin(lats2), np.nanmax(lats2)

res = args.res

width = int(np.ceil((lon_max - lon_min) / res))

height = int(np.ceil((lat_max - lat_min) / res))

print(f"Griglia di uscita: {width} x {height} pixel, passo {res} gradi")

area = geometry.AreaDefinition(

"out", "out", "out", "EPSG:4326", width, height,

(lon_min, lat_min, lon_max, lat_max))

swath = geometry.SwathDefinition(lons=lons2, lats=lats2)

arr = np.stack([

kd_tree.resample_nearest(swath, b, area, radius_of_influence=args.radius,

fill_value=np.nan)

for b in data

]).astype("float32")

transform = from_origin(lon_min, lat_max, res, res)

with rasterio.open(out, "w", driver="GTiff", height=arr.shape[1], width=arr.shape[2],

count=arr.shape[0], dtype="float32", crs="EPSG:4326",

transform=transform, nodata=np.nan, compress="deflate") as dst:

dst.write(arr)

print("Scritto", out)

if __name__ == "__main__":

main()

Il problema e' che il picco  di misura e' presente anche nelle argille e nei carbonati quindi si possono avere falsi positivi su suoli nudi anche coperture di edifici (pannelli solari compresi)

 Per mitigare si puo' usare una maschera  usando Esa World Cover....migliore ma non di tantissimo

Per scaricare Esa World Cover per la Toscana si puo' usare questo bucket AWS
 

aws s3 cp --no-sign-request s3://esa-worldcover/v200/2021/map/ESA_WorldCover_10m_2021_v200_N42E009_Map.tif .
aws s3 cp --no-sign-request s3://esa-worldcover/v200/2021/map/ESA_WorldCover_10m_2021_v200_N42E012_Map.tif .

 

e questo script 

python bare_soil_mask.py EMIT_L2A_RFL_..._003.nc \
  --worldcover ESA_WorldCover_10m_2021_v200_N42E009_Map.tif ESA_WorldCover_10m_2021_v200_N42E012_Map.tif \
  --builtup-buffer 1 --crop-bare-only --min-pixels 3 

 

#!/usr/bin/env python3

"""Maschera di suolo nudo / roccia da riflettanza EMIT L2A (EMIT_L2A_RFL_*.nc).

Criteri (tutti da tarare sulla scena):

- vegetazione verde assente: NDVI < --ndvi-max

- acqua assente: NDWI < --ndwi-max

- superficie non scura: riflettanza media 2100-2300 nm > --min-swir

- (opzionale) vegetazione secca esclusa con CAI < --cai-max

- (opzionale) --mask-nc: esclude nuvole/cirri/acqua/nuvole dilatate dalla maschera L2A

Uso:

python bare_soil_mask.py EMIT_L2A_RFL_001_20250808T083836_2522006_003.nc

python bare_soil_mask.py RFL.nc --mask-nc EMIT_L2A_MASK_001_..._003.nc --ndvi-max 0.2

Output:

<nome>_baresoil.tif GeoTIFF EPSG:4326 (GLT di EMIT), 1 = suolo nudo/roccia, 0 = altro, 255 = no dato

<nome>_baresoil_sensor.nc stessa maschera in geometria sensore (downtrack x crosstrack)

Richiede: xarray, netCDF4, numpy, rasterio

"""

import argparse

import os

import numpy as np

import rasterio

import xarray as xr

from rasterio.transform import Affine

NODATA = -9999

def band_idx(wvl, target):

return int(np.argmin(np.abs(wvl - target)))

def load_bands(ds, wvl, centers, half_width=10.0):

"""Riflettanza media attorno a ciascun centro (nm), leggendo solo le bande necessarie."""

out = {}

for name, c in centers.items():

sel = np.where(np.abs(wvl - c) <= half_width)[0]

if sel.size == 0:

sel = np.array([band_idx(wvl, c)])

arr = ds["reflectance"].isel(bands=sel).values.astype("float32")

arr[arr == NODATA] = np.nan

out[name] = np.nanmean(arr, axis=-1)

return out

def swir_mean(ds, wvl, lo=2100.0, hi=2300.0):

sel = np.where((wvl >= lo) & (wvl <= hi))[0]

arr = ds["reflectance"].isel(bands=sel).values.astype("float32")

arr[arr == NODATA] = np.nan

return np.nanmean(arr, axis=-1)

def main():

p = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)

p.add_argument("rfl", help="EMIT_L2A_RFL_*.nc")

p.add_argument("--mask-nc", default=None, help="EMIT_L2A_MASK_*.nc (nuvole, cirri, acqua)")

p.add_argument("--ndvi-max", type=float, default=0.25)

p.add_argument("--ndwi-max", type=float, default=0.0)

p.add_argument("--min-swir", type=float, default=0.05)

p.add_argument("--cai-max", type=float, default=None,

help="se impostato, esclude vegetazione secca con CAI >= valore (es. 0.01)")

p.add_argument("--out", default=None, help="prefisso di output")

args = p.parse_args()

ds = xr.open_dataset(args.rfl)

wvl = xr.open_dataset(args.rfl, group="sensor_band_parameters")["wavelengths"].values.astype(float)

loc = xr.open_dataset(args.rfl, group="location")

gt = ds.attrs["geotransform"]

centers = {"green": 560, "red": 650, "nir": 860}

if args.cai_max is not None:

centers.update({"s2030": 2030, "s2100": 2100, "s2210": 2210})

b = load_bands(ds, wvl, centers)

swir = swir_mean(ds, wvl)

eps = 1e-6

ndvi = (b["nir"] - b["red"]) / (b["nir"] + b["red"] + eps)

ndwi = (b["green"] - b["nir"]) / (b["green"] + b["nir"] + eps)

valid = np.isfinite(ndvi) & np.isfinite(ndwi) & np.isfinite(swir)

bare = valid & (ndvi < args.ndvi_max) & (ndwi < args.ndwi_max) & (swir > args.min_swir)

if args.cai_max is not None:

cai = 0.5 * (b["s2030"] + b["s2210"]) - b["s2100"]

bare &= np.isfinite(cai) & (cai < args.cai_max)

if args.mask_nc:

m = xr.open_dataset(args.mask_nc)["mask"].values # (down, cross, flags)

# indici dei flag EMIT L2A MASK (da verificare nel file): 0 nuvola, 1 cirri, 2 acqua, 4 nuvola dilatata

bad = np.zeros(bare.shape, dtype=bool)

for k in (0, 1, 2, 4):

if k < m.shape[-1]:

bad |= m[..., k] > 0

bare &= ~bad

pct = 100.0 * bare.sum() / max(valid.sum(), 1)

print(f"Suolo nudo/roccia: {pct:.1f}% dei pixel validi")

mask = np.where(valid, bare.astype("uint8"), 255).astype("uint8")

prefix = args.out or os.path.splitext(args.rfl)[0]

# maschera in geometria sensore

xr.Dataset({"baresoil": (("downtrack", "crosstrack"), mask)}).to_netcdf(prefix + "_baresoil_sensor.nc")

# ortorettifica con la GLT di EMIT (indici 1-based, 0 = nessun dato)

glt = np.nan_to_num(np.stack([loc["glt_x"].values, loc["glt_y"].values], axis=-1), nan=0).astype(int)

ok = np.all(glt != 0, axis=-1)

glt[ok] -= 1

out = np.full(glt.shape[:2], 255, dtype="uint8")

out[ok] = mask[glt[ok, 1], glt[ok, 0]]

with rasterio.open(prefix + "_baresoil.tif", "w", driver="GTiff", height=out.shape[0],

width=out.shape[1], count=1, dtype="uint8", crs="EPSG:4326",

transform=Affine.from_gdal(*gt), nodata=255, compress="deflate") as dst:

dst.write(out, 1)

print("Scritti:", prefix + "_baresoil.tif", "e", prefix + "_baresoil_sensor.nc")

if __name__ == "__main__":

main()