Rendering altezza dello zero termico

[spanezz]

📜 Documento di analisi nel wiki 📜

È un prodotto di solo contouring, sola altezza dello zero termico, per tutte le scadenze del modello.

I settaggi sul contouring sono recuperabili da: https://github.com/ARPA-SIMC/magics-maps/blob/master/json_fields/mm_fields.json

[brancomat]

sigla: alt0

parametri necessari (il livello teoricamente potrebbe essere omesso):
product:GRIB1,80,201,84;level:GRIB1,4 (gruppo COSMO)
product:GRIB1,98,228,24;level:GRIB1,1 (IFS-ECMWF)

esempio di query:

Altezza zero termico:
arki-query --data 'reftime:=today 0:00;product:GRIB1,80,201,84;level:GRIB1,4' http://arkiope.metarpa:8090/dataset/cosmo_5M_ita > hzero_cosmo.grib
arki-query --data 'reftime:=today 0:00;product:GRIB1,98,228,24;level:GRIB1,1' http://arkiope.metarpa:8090/dataset/ifs_ita010 > hzero_ecmwf.grib

contrariamente a quanto detto, è necessario un preprocessing per i dati ecmwf:
cosmo misura l'altezza dal livello del mare, ecmwf dal terreno, ai dati di ecmwf va sommata l'orografia prima dei plottaggi.

Il grib dell'orografia si può ottenere con questa query (i due dataset ifs deterministici hanno aree e risoluzioni diverse):

arki-query --data 'reftime:=today 0:00;product:GRIB1,98,128,129;level:GRIB1,1' http://arkiope.metarpa:8090/dataset/ifs_ita010 > orografia_ita.grib
arki-query --data 'reftime:=today 0:00;product:GRIB1,98,128,129;level:GRIB1,1' http://arkiope.metarpa:8090/dataset/ifs_atl025 > orografia_atl.grib

(@dcesari mi corregga se sbaglio)

L'orografia cambia molto raramente quindi si potrebbe integrare nella procedura, anche perché essendo un postprocessing non si riesce a garantire che chi si scarica lo zero termico si porti dietro anche quel campo.

[spanezz]

Ok, quindi in teoria l'orografia è un normale dato che possiamo trovare nell'input, sebbene non dipendente dallo step.

Però questo vorrebbe dire che per avere questo prodotto bisogna sempre far query, per ECMWF, sia dell'altezza dello zero termico, sia dell'orografia, cosa che una persona non è portata a fare.

arkimaps però non ha influenza sulle query che vengono fatte, e lavora coi dati che ha.

L'orografia richiesta a sua volta dipende dall'area dei dati, quindi non posso fare come con gli shapefile e tenermene una copia statica da usare se manca nell'input: se sto processando ita mi serve l'orografia ita, e se sto processando atl mi serve l'orografia atl.

Questo vuole anche dire che in caso di caching o in qualche modo riuso di orografie, serve implementare una verifica di allineamento tra l'area dei grib dell'altezza dello zero termico, e l'area dei grib dell'orografia.

Si potrebbe fare una funzione di precaching delle orografie. Questo introdurrebbe l'uso di una directory locale persistente tra invocazioni di arkimaps, la cui posizione andrebbe configurata in qualche modo. Questo introdurrebbe l'uso di un file di configurazione di arkimaps, che implica decidere in quali luoghi tenerlo e cercarlo.

Implementativamente, la cosa piú semplice che mi viene in mente è mettere in cima alla ricetta per l'altezza dello zero termico un commento che dice che per ECMWF serve piú query di quello che uno si aspetti. Posso, volendo, implementare un modo per avere commenti nello YAML che vengono mostrati anche nella documentazione della ricetta.

Può essere quest'ultima cosa un buon compromesso per il momento?

[dcesari]

Ops, vedo che sono stato invocato più di un anno fa; attenzione, ci sono due variabili di ECMWF:

• https://apps.ecmwf.int/codes/grib/param-db/?id=228024
• https://apps.ecmwf.int/codes/grib/param-db/?id=201084

tu hai indicato la prima, 228/24 che è in effetti l'altezza sopra la superficie ma qui stiamo in realtà elaborando la ricetta per la seconda. Intuisco che il motivo sia che nei nostri archivi abbiamo solo la prima per cause che non si possono spiegare in un chat pubblica, però prima verifichiamo. Se proprio dobbiamo farlo, la proposta di @spanezz mi pare accettabile, con il caveat che non è facilissimo per un utente capire se l'orografia di cui è in possesso sia allineata con i dati, non so se si può controllare il generatingProcessIdentifier, altrimenti bisogna frugare nel sito di ECMWF. Insomma decliniamo ogni responsabilità sulla veridicità del rendering sulle motagne.

[spanezz]

Davo per scontato che l'orografia fosse presente in ogni output del modello, per cui se faccio query su ifs_ita010 di:
reftime:=today 0:00;product:GRIB1,98,228,24;level:GRIB1,1 allora trovo anche reftime:=today 0:00;product:GRIB1,98,128,129;level:GRIB1,1 con l'orografia.

È giusto, o l'orografia non sempre è presente nell'output di un run di modello?

[dcesari]

Scusate, ritiro la mia polemica, il parametro 201/84 è quello della tabella DWD che figura anche nel sito di ECMWF, quindi è vero che ECMWF, per scelta, produce l'altezza dello zero termico rispetto alla superficie terrestre.

Riguardo l'ultima domanda, sarebbe ragionevole includerla sempre, aggiornata ad ogni run, ma a volte non è così nei nostri dataset ECMWF (verificato di recente con @edigiacomo)

[spanezz]

Ok, la situazione si complica. come in #24 ho creato una pagina di wiki per redigere un documento di analisi, qui: alt0.

Vi invito a metterci mano, e commentare le modifiche qui

[brancomat]

Io semplificherei di molto la questione: l'orografia di IFS viene calcolata solo se sono presenti sia altezza zero termico da terra (product:GRIB1,98,228,24;level:GRIB1,1) sia l'orografia (product:GRIB1,98,128,129;level:GRIB1,1), in caso contrario non viene plottato nulla.
Appuntiamo la cosa con un commento nella ricetta e buonanotte.

Da quel che vedo nei dataset le due variabili sono o presenti entrambe o assenti entrambe e al momento l'altezza dello zero termico viene plottata solo per ifs_ita010 (l'unico che ce l'ha), è un caso marginale per cui non andrei a complicarci la vita.

Aggiorno il documento wiki di conseguenza.

[spanezz]

La somma di altezza zero termico e orografia si può fare con vol7d o la implemento in Python?

[dcesari]

Propongo di farla in python, vol7d fa solo conversioni su dati che hanno gli stessi metadati ma coinvolgono variabili diverse, qui mi pare il caso opposto.

[brancomat]

In altre procedure lo faccio tramite eccodes, magari ci son modi migliori, comunque incollo snippet.

Nota: GRIB1,98,128,129 in realtà non è l'orografia è il geopotenziale, l'orografia va calcolata da quello dividendolo per 9.8 (vedi sotto) che è considerata approssimazione accettabile (credo, se @dcesari non mi cazzia)

geop = open("z_ita125.grb")
gid_geop = grib_new_from_file(geop)
orografia = grib_get_values(gid_geop)
orografia = orografia / 9.8
oro.close()

for gid in gid_list:
    clone_id = grib_clone(gid)
    values = grib_get_values(gid)
    # sommo i valori all'orografia
    values = values + orografia
    # setto a invalid i valori sotto l'orgrafia (zero)
    values[values <= orografia] = invalid

[spanezz]

Ok, allora seguo la stessa strategia proposta in #38

[spanezz]

Però con la differenza che il geopotenziale è lo stesso per tutti gli step

[spanezz]

Ho committato una bozza di ricetta.

Mi sono permesso di usare 9.80665 come in wreport e wikipedia.

Ho introdotto grib_set nell'input come discusso in #38:

  - model: ifs
    type: groundtomsl
    inputs: [z, alt0ground]
    grib_set:
      shortName: deg0l

La ricetta l'ho fatta copiando lo yaml da t2m, e immagino vada rivista

[spanezz]

Rispetto a questo pezzo del tuo snippet, cosa uso per invalid?

    # setto a invalid i valori sotto l'orgrafia (zero)
    values[values <= orografia] = invalid

[brancomat]

sorry, mancava un pezzo:

invalid=-999

[spanezz]

Ho fatto un prototipo di patch un po' di testa mia: 0983e4b

Siccome non mi piace hardcodare -999, l'ho messo come parametro invalid all'input: se presente, fa il capping.

Nel capping ho poi cambiato <= in <: vals[vals < z] = self.invalid, perché mi sembrava strano che "al suolo" non fosse considerato valido, però questo a sentimento di una persona profana (in questo momento sto ascoltando i Weather Report: conta?)

Propongo di partire da qui e mettere a punto con l'input degli esperti