Ricerca immagini Enmap

La ricerca di EoWeb non mi hai soddisfatto moltissimo

Ho provato a chiedere a Claude un tool per avere una selezione di immagini per posizione geografia, data e copertura nuvolosa

 

 e non mi ero mai accorto di quanto fosse cosi' poco diffusa la copertura di Enmap (che non ha una acquisizione continua ma solo on demand)

 

"""

Download EnMAP L2A per l'Italia in agosto, con filtro nuvolosità lato client,

calcolo della copertura per anno e cumulata, e export dei buchi.

Installazione:

    pip install pystac-client geopandas shapely requests pyogrio

Credenziali EOWEB (PowerShell, poi riapri il terminale):

    setx EOWEB_USER "tuo_utente"

    setx EOWEB_PASS "tua_password"

"""

import os

import sys

import warnings

import requests

import geopandas as gpd

from pystac_client import Client

from shapely.geometry import shape

from shapely.ops import unary_union

warnings.filterwarnings("ignore", message=".*does not conform to.*")

# ---------------- CONFIGURAZIONE ----------------

STAC_URL   = "https://geoservice.dlr.de/eoc/ogc/stac/v1/"

COLLECTION = "ENMAP_HSI_L2A"          # verificata a runtime (vedi sotto)

BBOX       = [6.6, 35.4, 18.6, 47.2]  # Italia con isole

YEARS      = [2023, 2024, 2025]       # il 2022 era in commissioning (poche scene)

MONTH_FROM = "01-01"

MONTH_TO   = "12-31"

MAX_CLOUD  = 5                         # % nuvolosità massima

OUT_DIR    = "enmap"

DRY_RUN    = True                      # True = solo analisi, nessun download

ASSET_KEYS = None                      # es. ["spectral_image"]; None = li stampa e basta

# -------------------------------------------------

NE_URL = ("https://naturalearth.s3.amazonaws.com/50m_cultural/"

          "ne_50m_admin_0_countries.zip")

def get_italy():

    world = gpd.read_file(NE_URL)

    col = "ADMIN" if "ADMIN" in world.columns else "admin"

    italy = world[world[col] == "Italy"].to_crs(3035)

    geom = unary_union(italy.geometry.values)

    return geom

def cloud_of(item):

    return item.properties.get("eo:cloud_cover")

def cloud_ok(item):

    cc = cloud_of(item)

    return cc is not None and cc <= MAX_CLOUD

def footprint_3035(items):

    g = gpd.GeoSeries([shape(i.geometry) for i in items], crs=4326).to_crs(3035)

    return unary_union(g.values)

def main():

    cat = Client.open(STAC_URL)

    # controllo collezione

    ids = [c.id for c in cat.get_collections()]

    enmap_ids = [c for c in ids if "ENMAP" in c.upper()]

    print("Collezioni EnMAP disponibili:", enmap_ids)

    coll = COLLECTION if COLLECTION in ids else (enmap_ids[0] if enmap_ids else None)

    if coll is None:

        sys.exit("Nessuna collezione EnMAP trovata nel catalogo.")

    if coll != COLLECTION:

        print(f"ATTENZIONE: uso '{coll}' invece di '{COLLECTION}'")

    italy = get_italy()

    italy_area = italy.area

    selected = {}      # anno -> lista item filtrati

    rows = []          # per export footprint

    for y in YEARS:

        search = cat.search(

            collections=[coll],

            bbox=BBOX,

            datetime=f"{y}-{MONTH_FROM}/{y}-{MONTH_TO}",

            max_items=5000,

        )

        items = list(search.items())

        n_tot = len(items)

        n_nocc = sum(cloud_of(i) is None for i in items)

        good = [i for i in items if cloud_ok(i)]

        selected[y] = good

        print(f"\n{y}: {n_tot} scene totali | "

              f"{len(good)} con nuvole <= {MAX_CLOUD}% | "

              f"{n_nocc} senza valore nuvolosità")

        if n_tot and n_nocc == n_tot:

            print("  -> il campo 'eo:cloud_cover' non esiste. Proprietà disponibili:")

            print("  ", sorted(items[0].properties.keys()))

        if good:

            cov = footprint_3035(good).intersection(italy)

            print(f"  copertura Italia: {cov.area / italy_area:.1%}")

        else:

            print("  copertura Italia: 0%")

        for i in good:

            rows.append({

                "id": i.id, "anno": y,

                "data": str(i.datetime), "nuvole": cloud_of(i),

                "geometry": shape(i.geometry),

            })

    all_items = [i for its in selected.values() for i in its]

    if not all_items:

        sys.exit("\nNessuna scena selezionata: alza MAX_CLOUD o allarga le date.")

    # copertura cumulata e buchi

    cov_all = footprint_3035(all_items).intersection(italy)

    print(f"\nCopertura cumulata {YEARS[0]}-{YEARS[-1]}: "

          f"{cov_all.area / italy_area:.1%}  ({len(all_items)} scene)")

    gaps = italy.difference(cov_all)

    gpd.GeoDataFrame({"tipo": ["buco"]}, geometry=[gaps], crs=3035)\

       .to_file("buchi_italia.gpkg", driver="GPKG")

    gpd.GeoDataFrame(rows, crs=4326).to_file("footprint_scene.gpkg", driver="GPKG")

    print("Scritti: buchi_italia.gpkg, footprint_scene.gpkg (aprili in QGIS)")

    # asset disponibili

    print("\nAsset della prima scena:")

    for k, a in all_items[0].assets.items():

        print(f"  {k:25s} {a.href}")

    if DRY_RUN:

        print("\nDRY_RUN attivo: nessun download. Imposta DRY_RUN = False "

              "e ASSET_KEYS dopo aver guardato gli asset sopra.")

        return

    if not ASSET_KEYS:

        sys.exit("Imposta ASSET_KEYS con le chiavi degli asset da scaricare.")

    # download

    user, pwd = os.environ.get("EOWEB_USER"), os.environ.get("EOWEB_PASS")

    if not user or not pwd:

        sys.exit("Imposta EOWEB_USER e EOWEB_PASS.")

    os.makedirs(OUT_DIR, exist_ok=True)

    sess = requests.Session()

    sess.auth = (user, pwd)

    for it in all_items:

        for key in ASSET_KEYS:

            asset = it.assets.get(key)

            if asset is None:

                continue

            ext = os.path.splitext(asset.href.split("?")[0])[1] or ".bin"

            out = os.path.join(OUT_DIR, f"{it.id}_{key}{ext}")

            if os.path.exists(out):

                print("già presente:", out)

                continue

            print("scarico:", out)

            tmp = out + ".part"

            with sess.get(asset.href, stream=True, timeout=120) as r:

                r.raise_for_status()

                with open(tmp, "wb") as f:

                    for chunk in r.iter_content(1 << 20):

                        f.write(chunk)

            os.replace(tmp, out)

if __name__ == "__main__":

    main()