import csv
import json
import re
from pathlib import Path
import earthaccess
import geopandas as gpd
from shapely.geometry import Polygon, box
from shapely.ops import unary_union
# ---------- parametri ----------
VERSION = "002"
CLOUD_MAX = 10
TEMPORAL = ("2022-08-01", "2026-10-18")
MIN_LAND = 0.20
MAX_OVERLAP = 0.20
OUT_DIR = Path("./emit_iberia")
DRY_RUN = False
# Regioni di ricerca (lon_min, lat_min, lon_max, lat_max).
# La penisola iberica include anche le Baleari; le isole atlantiche sono opzionali.
REGIONI = {
"penisola_baleari": (-9.6, 35.95, 4.5, 43.9), # Spagna + Portogallo continentali + Baleari
}
INCLUDI_CANARIE = False # Spagna: Canarie
INCLUDI_MADEIRA = False # Portogallo: Madeira
INCLUDI_AZZORRE = False # Portogallo: Azzorre
if INCLUDI_CANARIE:
REGIONI["canarie"] = (-18.4, 27.5, -13.3, 29.5)
if INCLUDI_MADEIRA:
REGIONI["madeira"] = (-17.4, 32.35, -16.2, 33.2)
if INCLUDI_AZZORRE:
REGIONI["azzorre"] = (-31.4, 36.8, -24.9, 39.8)
GET_RFL = True # EMIT_L2A_RFL
GET_MASK = False # EMIT_L2A_MASK
GET_UNCERT = False # EMIT_L2A_RFLUNCERT
GET_OBS = False # EMIT_L1B_OBS (dal prodotto EMITL1BRAD)
GET_RAD = False # EMIT_L1B_RAD (file grande, di solito non serve)
GET_BROWSE = True # anteprima PNG
GET_META = True # JSON UMM-G + CSV riassuntivo
# ---------- contorni Spagna e Portogallo ----------
# Il poligono Natural Earth di ciascun paese comprende anche i territori lontani (Canarie, Azzorre,
# Madeira, Ceuta, Melilla): li teniamo solo se ricadono in una delle REGIONI scelte.
URL = "https://naturalearth.s3.amazonaws.com/10m_cultural/ne_10m_admin_0_countries.zip"
world = gpd.read_file(URL)
area = unary_union([box(*b) for b in REGIONI.values()])
PAESI = {"ESP": "spagna", "PRT": "portogallo"}
paesi = {}
for a3 in PAESI:
geom = (world[world["ADM0_A3"] == a3]
.to_crs(4326)
.union_all()
.intersection(area)
.buffer(0.05)) # margine per i granuli costieri
assert not geom.is_empty, f"nessuna geometria trovata per {a3}"
paesi[a3] = geom
iberia = unary_union(list(paesi.values()))
earthaccess.login(strategy="netrc") # oppure login() / login(persist=True)
# ---------- ricerca (una per regione, poi unione senza duplicati) ----------
def search_all(short_name, **kw):
found = {}
for nome, bbox in REGIONI.items():
res = earthaccess.search_data(
short_name=short_name,
version=VERSION,
bounding_box=bbox,
temporal=TEMPORAL,
count=-1,
**kw,
)
print(f" {short_name} [{nome}]: {len(res)} granuli")
for g in res:
found[g["umm"]["GranuleUR"]] = g
return list(found.values())
l2a = search_all("EMITL2ARFL", cloud_cover=(0, CLOUD_MAX))
print(f"{len(l2a)} granuli L2A nelle regioni scelte")
def browse_links(g):
d = g.dataviz_links
d = d() if callable(d) else d # proprietà o metodo, a seconda della versione
return [u for u in (d or []) if u.endswith(".png")]
def footprint(g):
pts = (g["umm"]["SpatialExtent"]["HorizontalSpatialDomain"]
["Geometry"]["GPolygons"][0]["Boundary"]["Points"])
return Polygon([(p["Longitude"], p["Latitude"]) for p in pts])
def cloud(g):
return float(g["umm"].get("CloudCover", 100))
def stamp(g):
"""Timestamp di acquisizione YYYYMMDDTHHMMSS dal GranuleUR."""
m = re.search(r"(\d{8}T\d{6})", g["umm"]["GranuleUR"])
return m.group(1) if m else None
# ---------- filtro terra + greedy ----------
cands = []
for g in l2a:
poly = footprint(g)
land = poly.intersection(iberia)
if land.area / poly.area >= MIN_LAND:
cands.append((g, land, land.area / poly.area))
print(f"{len(cands)} granuli con almeno {MIN_LAND:.0%} su terra spagnola o portoghese")
cands.sort(key=lambda x: cloud(x[0]))
selected, covered = [], None
for g, land, frac in cands:
overlap = 0.0 if covered is None else land.intersection(covered).area / land.area
if overlap <= MAX_OVERLAP:
selected.append((g, frac))
covered = land if covered is None else unary_union([covered, land])
print(f"{len(selected)} granuli selezionati")
if covered is not None:
print(f"copertura Spagna + Portogallo: {covered.intersection(iberia).area / iberia.area:.1%}")
for a3, nome in PAESI.items():
print(f" {nome}: {covered.intersection(paesi[a3]).area / paesi[a3].area:.1%}")
# ---------- ricerca L1B (RAD + OBS) e abbinamento per timestamp ----------
l1b_by_stamp = {}
if GET_OBS or GET_RAD:
l1b = search_all("EMITL1BRAD")
print(f"{len(l1b)} granuli L1B nelle regioni scelte")
l1b_by_stamp = {stamp(g): g for g in l1b if stamp(g)}
# ---------- costruzione lista file ----------
def pick(links, key):
return [u for u in links if key in u.split("/")[-1] and u.endswith(".nc")]
links, rows, missing_l1b = [], [], []
for g, frac in selected:
ts = stamp(g)
fp = footprint(g)
l2a_links = g.data_links()
if GET_RFL: links += pick(l2a_links, "_RFL_")
if GET_MASK: links += pick(l2a_links, "_MASK_")
if GET_UNCERT: links += pick(l2a_links, "RFLUNCERT")
if GET_BROWSE: links += browse_links(g)
g1 = l1b_by_stamp.get(ts)
if g1 is None:
if GET_OBS or GET_RAD:
missing_l1b.append(g["umm"]["GranuleUR"])
else:
l1b_links = g1.data_links()
if GET_OBS: links += pick(l1b_links, "_OBS_")
if GET_RAD: links += pick(l1b_links, "_RAD_")
rows.append({
"granule_ur": g["umm"]["GranuleUR"],
"timestamp": ts,
"begin": g["umm"]["TemporalExtent"]["RangeDateTime"]["BeginningDateTime"],
"end": g["umm"]["TemporalExtent"]["RangeDateTime"].get("EndingDateTime", ""),
"cloud_pct": cloud(g),
"land_fraction": round(frac, 3),
"land_spagna": round(fp.intersection(paesi["ESP"]).area / fp.area, 3),
"land_portogallo": round(fp.intersection(paesi["PRT"]).area / fp.area, 3),
"l1b_ur": g1["umm"]["GranuleUR"] if g1 else "",
"footprint_wkt": fp.wkt,
})
links = list(dict.fromkeys(links)) # elimina duplicati mantenendo l'ordine
print(f"\n{len(links)} file da scaricare")
for r in rows:
print(r["granule_ur"], f"cloud={r['cloud_pct']:.1f}%",
f"ES={r['land_spagna']:.0%} PT={r['land_portogallo']:.0%}",
"| L1B:", r["l1b_ur"] or "NON TROVATO")
if missing_l1b:
print(f"\nATTENZIONE: nessun L1B abbinato per {len(missing_l1b)} granuli")
try:
tot = sum(g.size() for g, _ in selected)
print(f"\nDimensione L2A (indicativa, tutti i file): {tot/1024:.1f} GB")
except Exception:
pass
# ---------- metadati ----------
if GET_META and not DRY_RUN and rows:
meta_dir = OUT_DIR / "metadata"
meta_dir.mkdir(parents=True, exist_ok=True)
for g, _ in selected:
ur = g["umm"]["GranuleUR"]
(meta_dir / f"{ur}.umm.json").write_text(
json.dumps({"umm": g["umm"], "meta": g["meta"]}, indent=2))
g1 = l1b_by_stamp.get(stamp(g))
if g1 is not None:
(meta_dir / f"{g1['umm']['GranuleUR']}.umm.json").write_text(
json.dumps({"umm": g1["umm"], "meta": g1["meta"]}, indent=2))
with open(OUT_DIR / "granules_summary.csv", "w", newline="") as f:
w = csv.DictWriter(f, fieldnames=rows[0].keys())
w.writeheader()
w.writerows(rows)
# ---------- download ----------
if not DRY_RUN:
OUT_DIR.mkdir(parents=True, exist_ok=True)
earthaccess.download(links, str(OUT_DIR))
