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
Le condizioni migliori di rilevabilita' sono costituite da un background molto riflettente (quindi no mare dove peraltro i riflessi potrebbero ingannare il retrival)
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)
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()
Visto che l'algoritmo e' basata sul background ho provato a vedere cosa accade se,invece di usare tutta l'immagine, si usa solo un taglio rispetto al bersaglio
Per tagliare una immagine Emit si puo' usare
python emit_crop.py EMIT_L1B_RAD_001_20250808T083836_2522006_003.nc --point 43.56726 10.80698 --point-km 1 --margin-rows 50 --margin-cols 50 --out-dir crop/
#!/usr/bin/env python3
"""
emit_crop.py - Ritaglia un granulo EMIT L1B (RAD, ed eventualmente OBS) sull'area di interesse.
Di default il taglio avviene SOLO lungo la traccia (downtrack) e mantiene tutta la
larghezza (crosstrack): si conservano tutte le colonne del sensore (utile per MF
column-wise). Con --margin-cols (o --cols) si taglia anche in crosstrack attorno alla
bbox/punto. La struttura dei gruppi NetCDF resta quella originale, quindi il file
ridotto puo' essere letto come un normale granulo EMIT.
Esempi
------
# bbox in gradi: lon_min lat_min lon_max lat_max
python emit_crop.py EMIT_L1B_RAD_001_20240801T101010_xxxx.nc \
--bbox 10.5 43.6 10.9 43.9 --margin-rows 150 --out-dir crop/
# punto di interesse (lat lon) con bbox di +/- 1 km e 150 righe di margine
python emit_crop.py EMIT_L1B_RAD_....nc --point 43.56726 10.80698 --margin-rows 150
# taglio anche in crosstrack: 150 righe e 200 colonne di margine attorno al punto
python emit_crop.py EMIT_L1B_RAD_....nc --point 43.56726 10.80698 \
--margin-rows 150 --margin-cols 200
# se conosci gia' righe (e colonne)
python emit_crop.py EMIT_L1B_RAD_....nc --rows 400 1200 --cols 300 900 --out-dir crop/
# OBS indicato a mano (altrimenti viene cercato nella stessa cartella)
python emit_crop.py RAD.nc --obs OBS.nc --bbox ...
Gestione della GLT (glt_x / glt_y, nel gruppo "location"):
--glt adjust (default) gli indici vengono riallineati alle nuove righe; i pixel
fuori dal ritaglio diventano "no data". La griglia
ortorettificata resta quella originale.
--glt drop glt_x e glt_y vengono eliminate (utile se il reader usa
solo lat/lon per pixel).
--glt keep copiate cosi' come sono (SCONSIGLIATO: gli indici non
corrispondono piu' alle righe).
Requisiti: numpy, netCDF4
"""
import argparse
import sys
from datetime import datetime, timezone
from pathlib import Path
import numpy as np
import netCDF4 as nc
# ----------------------------------------------------------------------------
# Ricerca delle righe
# ----------------------------------------------------------------------------
def find_window(path, bbox, margin_rows, margin_cols):
"""Restituisce (r0, r1, c0, c1), con r1 e c1 esclusivi.
margin_cols=None -> tutte le colonne (nessun taglio in crosstrack).
"""
lon_min, lat_min, lon_max, lat_max = bbox
with nc.Dataset(path) as ds:
loc = ds["location"]
lat = np.ma.filled(loc["lat"][:].astype("float64"), np.nan)
lon = np.ma.filled(loc["lon"][:].astype("float64"), np.nan)
nrows, ncols = lat.shape
inside = (lat >= lat_min) & (lat <= lat_max) & (lon >= lon_min) & (lon <= lon_max)
rows = np.flatnonzero(inside.any(axis=1))
cols = np.flatnonzero(inside.any(axis=0))
if rows.size == 0:
sys.exit("Nessun pixel del granulo cade dentro la bbox indicata.")
r0 = max(int(rows[0]) - margin_rows, 0)
r1 = min(int(rows[-1]) + 1 + margin_rows, nrows)
if margin_cols is None:
c0, c1 = 0, ncols
else:
c0 = max(int(cols[0]) - margin_cols, 0)
c1 = min(int(cols[-1]) + 1 + margin_cols, ncols)
return r0, r1, c0, c1
# ----------------------------------------------------------------------------
# Copia dei gruppi
# ----------------------------------------------------------------------------
def copy_attrs(src, dst, skip=()):
for a in src.ncattrs():
if a not in skip:
dst.setncattr(a, src.getncattr(a))
def adjusted_glt(group, r0, r1, c0, c1):
"""Riallinea glt_x / glt_y (indici 1-based su righe/colonne originali) al ritaglio."""
gx = group["glt_x"]
gy = group["glt_y"]
gx.set_auto_maskandscale(False)
gy.set_auto_maskandscale(False)
x = gx[:].astype("int64")
y = gy[:].astype("int64")
new_y = y - r0
new_x = x - c0
invalid = (y < 1) | (x < 1) | (new_y < 1) | (new_y > (r1 - r0)) | (new_x < 1) | (new_x > (c1 - c0))
fill_x = gx.getncattr("_FillValue") if "_FillValue" in gx.ncattrs() else 0
fill_y = gy.getncattr("_FillValue") if "_FillValue" in gy.ncattrs() else 0
out_y = np.where(invalid, fill_y, new_y).astype(gy.dtype)
out_x = np.where(invalid, fill_x, new_x).astype(gx.dtype)
return {"glt_x": out_x, "glt_y": out_y}
def copy_group(src, dst, r0, r1, c0, c1, glt_mode):
# dimensioni
for name, dim in src.dimensions.items():
if name == "downtrack":
size = r1 - r0
elif name == "crosstrack":
size = c1 - c0
else:
size = None if dim.isunlimited() else len(dim)
dst.createDimension(name, size)
# eventuale GLT riallineata
overrides = {}
if glt_mode == "adjust" and "glt_x" in src.variables and "glt_y" in src.variables:
overrides = adjusted_glt(src, r0, r1, c0, c1)
# variabili
for name, var in src.variables.items():
if glt_mode == "drop" and name in ("glt_x", "glt_y"):
continue
var.set_auto_maskandscale(False)
fill = var.getncattr("_FillValue") if "_FillValue" in var.ncattrs() else None
kwargs = {}
if var.dtype != str and var.dtype.kind != "O":
kwargs.update(zlib=True, complevel=4)
out = dst.createVariable(name, var.dtype, var.dimensions, fill_value=fill, **kwargs)
out.set_auto_maskandscale(False)
copy_attrs(var, out, skip=("_FillValue",))
if name in overrides:
out[:] = overrides[name]
continue
sl = {"downtrack": slice(r0, r1), "crosstrack": slice(c0, c1)}
idx = tuple(sl.get(d, slice(None)) for d in var.dimensions)
out[...] = var[idx]
# sottogruppi
for gname, g in src.groups.items():
copy_group(g, dst.createGroup(gname), r0, r1, c0, c1, glt_mode)
def update_geo_attrs(dst):
"""Aggiorna gli attributi globali di estensione geografica in base al ritaglio."""
loc = dst["location"]
lat = np.ma.masked_less_equal(loc["lat"][:].astype("float64"), -9990)
lon = np.ma.masked_less_equal(loc["lon"][:].astype("float64"), -9990)
mapping = {
"northernmost_latitude": lat.max(),
"southernmost_latitude": lat.min(),
"easternmost_longitude": lon.max(),
"westernmost_longitude": lon.min(),
}
for k, v in mapping.items():
if k in dst.ncattrs():
dst.setncattr(k, float(v))
def crop_file(path, out_path, r0, r1, c0, c1, glt_mode):
with nc.Dataset(path) as src:
nrows = src.dimensions["downtrack"].size
ncols = src.dimensions["crosstrack"].size
r1 = min(r1, nrows)
c1 = min(c1, ncols)
if r0 >= r1:
sys.exit(f"Intervallo righe non valido: {r0}-{r1} (downtrack={nrows})")
if c0 >= c1:
sys.exit(f"Intervallo colonne non valido: {c0}-{c1} (crosstrack={ncols})")
with nc.Dataset(out_path, "w", format="NETCDF4") as dst:
copy_attrs(src, dst)
copy_group(src, dst, r0, r1, c0, c1, glt_mode)
update_geo_attrs(dst)
stamp = datetime.now(timezone.utc).strftime("%Y-%m-%dT%H:%M:%SZ")
hist = (dst.getncattr("history") + "\n") if "history" in dst.ncattrs() else ""
dst.setncattr(
"history",
f"{hist}{stamp} emit_crop.py: downtrack rows {r0}:{r1} of {nrows}, "
f"crosstrack cols {c0}:{c1} of {ncols} (glt={glt_mode})",
)
print(f"OK {Path(path).name} -> {out_path} (righe {r0}:{r1}, {r1 - r0} di {nrows}; colonne {c0}:{c1}, {c1 - c0} di {ncols})")
# ----------------------------------------------------------------------------
# OBS
# ----------------------------------------------------------------------------
def guess_obs(rad_path):
p = Path(rad_path)
if "_RAD_" not in p.name:
return None
cand = p.with_name(p.name.replace("_RAD_", "_OBS_"))
return cand if cand.exists() else None
# ----------------------------------------------------------------------------
# main
# ----------------------------------------------------------------------------
def main():
ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
ap.add_argument("rad", help="file EMIT L1B RAD (.nc)")
ap.add_argument("--obs", help="file EMIT L1B OBS (.nc); se omesso viene cercato accanto al RAD")
ap.add_argument("--no-obs", action="store_true", help="non ritagliare il file OBS")
g = ap.add_mutually_exclusive_group(required=True)
g.add_argument("--bbox", nargs=4, type=float, metavar=("LON_MIN", "LAT_MIN", "LON_MAX", "LAT_MAX"))
g.add_argument("--point", nargs=2, type=float, metavar=("LAT", "LON"),
help="punto di interesse (lat lon): usa una bbox di +/- --point-km attorno ad esso")
g.add_argument("--rows", nargs=2, type=int, metavar=("R0", "R1"), help="righe downtrack [R0, R1)")
ap.add_argument("--cols", nargs=2, type=int, metavar=("C0", "C1"),
help="colonne crosstrack [C0, C1) (con --rows; default tutte)")
ap.add_argument("--point-km", type=float, default=1.0,
help="semi-lato della bbox attorno a --point, in km (default 1)")
ap.add_argument("--margin-rows", type=int, default=150,
help="righe extra sopra e sotto la bbox (default 150, ~9 km); solo con --bbox")
ap.add_argument("--margin-cols", type=int, default=None,
help="colonne extra a sinistra e destra della bbox/punto; se omesso si tengono "
"tutte le colonne (nessun taglio in crosstrack)")
ap.add_argument("--glt", choices=["adjust", "drop", "keep"], default="adjust")
ap.add_argument("--out-dir", default=".", help="cartella di output")
args = ap.parse_args()
out_dir = Path(args.out_dir)
out_dir.mkdir(parents=True, exist_ok=True)
if args.point:
lat_c, lon_c = args.point
dlat = args.point_km / 111.32
dlon = args.point_km / (111.32 * np.cos(np.radians(lat_c)))
args.bbox = [lon_c - dlon, lat_c - dlat, lon_c + dlon, lat_c + dlat]
if args.bbox:
r0, r1, c0, c1 = find_window(args.rad, args.bbox, args.margin_rows, args.margin_cols)
else:
r0, r1 = args.rows
c0, c1 = args.cols if args.cols else (0, 10**9)
print(f"Righe selezionate: {r0}:{r1} colonne: {c0}:{c1 if c1 < 10**9 else 'tutte'}")
rad = Path(args.rad)
crop_file(rad, out_dir / f"{rad.stem}_crop.nc", r0, r1, c0, c1, args.glt)
if not args.no_obs:
obs = Path(args.obs) if args.obs else guess_obs(rad)
if obs is None:
print("File OBS non trovato: salto (usa --obs per indicarlo).")
else:
with nc.Dataset(obs) as o, nc.Dataset(rad) as r:
if (o.dimensions["downtrack"].size != r.dimensions["downtrack"].size
or o.dimensions["crosstrack"].size != r.dimensions["crosstrack"].size):
sys.exit("RAD e OBS hanno dimensioni diverse: file non corrispondenti.")
crop_file(obs, out_dir / f"{obs.stem}_crop.nc", r0, r1, c0, c1, args.glt)
if __name__ == "__main__":
main()
Ma anche in questo caso non ho avuto risultati significativi
Utilizzando un crop su dati Enmap
python enmap_crop.py ./dati_crop/ENMAP01-____L1B-DT0000147167_20250731T104357Z_010_V010506_20261009T201416Z/ --point 43.56726 10.80698 --margin-rows 100 --margin-cols 100 --out-dir ./discarica_crop/
#!/usr/bin/env python3
"""
emit_crop.py - Ritaglia un granulo EMIT L1B (RAD, ed eventualmente OBS) sull'area di interesse.
Il taglio avviene SOLO lungo la traccia (downtrack) e mantiene tutta la larghezza
(crosstrack): si conservano tutte le colonne del sensore (utile per MF column-wise)
e la struttura dei gruppi NetCDF resta quella originale, quindi il file ridotto
puo' essere letto come un normale granulo EMIT.
Esempi
------
# bbox in gradi: lon_min lat_min lon_max lat_max
python emit_crop.py EMIT_L1B_RAD_001_20240801T101010_xxxx.nc \
--bbox 10.5 43.6 10.9 43.9 --margin-rows 150 --out-dir crop/
# punto di interesse (lat lon) con bbox di +/- 1 km e 150 righe di margine
python emit_crop.py EMIT_L1B_RAD_....nc --point 43.56726 10.80698 --margin-rows 150
# se conosci gia' le righe
python emit_crop.py EMIT_L1B_RAD_....nc --rows 400 1200 --out-dir crop/
# OBS indicato a mano (altrimenti viene cercato nella stessa cartella)
python emit_crop.py RAD.nc --obs OBS.nc --bbox ...
Gestione della GLT (glt_x / glt_y, nel gruppo "location"):
--glt adjust (default) gli indici vengono riallineati alle nuove righe; i pixel
fuori dal ritaglio diventano "no data". La griglia
ortorettificata resta quella originale.
--glt drop glt_x e glt_y vengono eliminate (utile se il reader usa
solo lat/lon per pixel).
--glt keep copiate cosi' come sono (SCONSIGLIATO: gli indici non
corrispondono piu' alle righe).
Requisiti: numpy, netCDF4
"""
import argparse
import sys
from datetime import datetime, timezone
from pathlib import Path
import numpy as np
import netCDF4 as nc
# ----------------------------------------------------------------------------
# Ricerca delle righe
# ----------------------------------------------------------------------------
def find_rows(path, bbox, margin):
"""Restituisce (r0, r1) con r1 esclusivo, usando lat/lon per pixel del gruppo location."""
lon_min, lat_min, lon_max, lat_max = bbox
with nc.Dataset(path) as ds:
loc = ds["location"]
lat = np.ma.filled(loc["lat"][:].astype("float64"), np.nan)
lon = np.ma.filled(loc["lon"][:].astype("float64"), np.nan)
nrows = lat.shape[0]
inside = (lat >= lat_min) & (lat <= lat_max) & (lon >= lon_min) & (lon <= lon_max)
rows = np.flatnonzero(inside.any(axis=1))
if rows.size == 0:
sys.exit("Nessun pixel del granulo cade dentro la bbox indicata.")
r0 = max(int(rows[0]) - margin, 0)
r1 = min(int(rows[-1]) + 1 + margin, nrows)
return r0, r1
# ----------------------------------------------------------------------------
# Copia dei gruppi
# ----------------------------------------------------------------------------
def copy_attrs(src, dst, skip=()):
for a in src.ncattrs():
if a not in skip:
dst.setncattr(a, src.getncattr(a))
def adjusted_glt(group, r0, r1):
"""Riallinea glt_x / glt_y (indici 1-based sulle righe/colonne originali) al ritaglio."""
gx = group["glt_x"]
gy = group["glt_y"]
gx.set_auto_maskandscale(False)
gy.set_auto_maskandscale(False)
x = gx[:].astype("int64")
y = gy[:].astype("int64")
new_y = y - r0
invalid = (y < 1) | (x < 1) | (new_y < 1) | (new_y > (r1 - r0))
fill_x = gx.getncattr("_FillValue") if "_FillValue" in gx.ncattrs() else 0
fill_y = gy.getncattr("_FillValue") if "_FillValue" in gy.ncattrs() else 0
out_y = np.where(invalid, fill_y, new_y).astype(gy.dtype)
out_x = np.where(invalid, fill_x, x).astype(gx.dtype)
return {"glt_x": out_x, "glt_y": out_y}
def copy_group(src, dst, r0, r1, glt_mode):
# dimensioni
for name, dim in src.dimensions.items():
if name == "downtrack":
size = r1 - r0
else:
size = None if dim.isunlimited() else len(dim)
dst.createDimension(name, size)
# eventuale GLT riallineata
overrides = {}
if glt_mode == "adjust" and "glt_x" in src.variables and "glt_y" in src.variables:
overrides = adjusted_glt(src, r0, r1)
# variabili
for name, var in src.variables.items():
if glt_mode == "drop" and name in ("glt_x", "glt_y"):
continue
var.set_auto_maskandscale(False)
fill = var.getncattr("_FillValue") if "_FillValue" in var.ncattrs() else None
kwargs = {}
if var.dtype != str and var.dtype.kind != "O":
kwargs.update(zlib=True, complevel=4)
out = dst.createVariable(name, var.dtype, var.dimensions, fill_value=fill, **kwargs)
out.set_auto_maskandscale(False)
copy_attrs(var, out, skip=("_FillValue",))
if name in overrides:
out[:] = overrides[name]
continue
idx = tuple(slice(r0, r1) if d == "downtrack" else slice(None) for d in var.dimensions)
out[...] = var[idx]
# sottogruppi
for gname, g in src.groups.items():
copy_group(g, dst.createGroup(gname), r0, r1, glt_mode)
def update_geo_attrs(dst):
"""Aggiorna gli attributi globali di estensione geografica in base al ritaglio."""
loc = dst["location"]
lat = np.ma.masked_less_equal(loc["lat"][:].astype("float64"), -9990)
lon = np.ma.masked_less_equal(loc["lon"][:].astype("float64"), -9990)
mapping = {
"northernmost_latitude": lat.max(),
"southernmost_latitude": lat.min(),
"easternmost_longitude": lon.max(),
"westernmost_longitude": lon.min(),
}
for k, v in mapping.items():
if k in dst.ncattrs():
dst.setncattr(k, float(v))
def crop_file(path, out_path, r0, r1, glt_mode):
with nc.Dataset(path) as src:
nrows = src.dimensions["downtrack"].size
r1 = min(r1, nrows)
if r0 >= r1:
sys.exit(f"Intervallo righe non valido: {r0}-{r1} (downtrack={nrows})")
with nc.Dataset(out_path, "w", format="NETCDF4") as dst:
copy_attrs(src, dst)
copy_group(src, dst, r0, r1, glt_mode)
update_geo_attrs(dst)
stamp = datetime.now(timezone.utc).strftime("%Y-%m-%dT%H:%M:%SZ")
hist = (dst.getncattr("history") + "\n") if "history" in dst.ncattrs() else ""
dst.setncattr(
"history",
f"{hist}{stamp} emit_crop.py: downtrack rows {r0}:{r1} of {nrows} (glt={glt_mode})",
)
print(f"OK {Path(path).name} -> {out_path} (righe {r0}:{r1}, {r1 - r0} di {nrows})")
# ----------------------------------------------------------------------------
# OBS
# ----------------------------------------------------------------------------
def guess_obs(rad_path):
p = Path(rad_path)
if "_RAD_" not in p.name:
return None
cand = p.with_name(p.name.replace("_RAD_", "_OBS_"))
return cand if cand.exists() else None
# ----------------------------------------------------------------------------
# main
# ----------------------------------------------------------------------------
def main():
ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
ap.add_argument("rad", help="file EMIT L1B RAD (.nc)")
ap.add_argument("--obs", help="file EMIT L1B OBS (.nc); se omesso viene cercato accanto al RAD")
ap.add_argument("--no-obs", action="store_true", help="non ritagliare il file OBS")
g = ap.add_mutually_exclusive_group(required=True)
g.add_argument("--bbox", nargs=4, type=float, metavar=("LON_MIN", "LAT_MIN", "LON_MAX", "LAT_MAX"))
g.add_argument("--point", nargs=2, type=float, metavar=("LAT", "LON"),
help="punto di interesse (lat lon): usa una bbox di +/- --point-km attorno ad esso")
g.add_argument("--rows", nargs=2, type=int, metavar=("R0", "R1"), help="righe downtrack [R0, R1)")
ap.add_argument("--point-km", type=float, default=1.0,
help="semi-lato della bbox attorno a --point, in km (default 1)")
ap.add_argument("--margin-rows", type=int, default=150,
help="righe extra sopra e sotto la bbox (default 150, ~9 km); solo con --bbox")
ap.add_argument("--glt", choices=["adjust", "drop", "keep"], default="adjust")
ap.add_argument("--out-dir", default=".", help="cartella di output")
args = ap.parse_args()
out_dir = Path(args.out_dir)
out_dir.mkdir(parents=True, exist_ok=True)
if args.point:
lat_c, lon_c = args.point
dlat = args.point_km / 111.32
dlon = args.point_km / (111.32 * np.cos(np.radians(lat_c)))
args.bbox = [lon_c - dlon, lat_c - dlat, lon_c + dlon, lat_c + dlat]
if args.bbox:
r0, r1 = find_rows(args.rad, args.bbox, args.margin_rows)
else:
r0, r1 = args.rows
print(f"Righe selezionate: {r0}:{r1}")
rad = Path(args.rad)
crop_file(rad, out_dir / f"{rad.stem}_crop.nc", r0, r1, args.glt)
if not args.no_obs:
obs = Path(args.obs) if args.obs else guess_obs(rad)
if obs is None:
print("File OBS non trovato: salto (usa --obs per indicarlo).")
else:
with nc.Dataset(obs) as o, nc.Dataset(rad) as r:
if o.dimensions["downtrack"].size != r.dimensions["downtrack"].size:
sys.exit("RAD e OBS hanno un numero diverso di righe: file non corrispondenti.")
crop_file(obs, out_dir / f"{obs.stem}_crop.nc", r0, r1, args.glt)
if __name__ == "__main__":
main()
L'elaborazione del solo crop dell'immagine completa sui dati Enmap ha migliorato il risultato in quanto ci sono molto meno falsi positivi a fronte di una migliore definizione del plume
