Ricerca immagini Emit (2)

Versione aggiornata che crea anche una preview 

 

 

 

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))