Wrong grib2 level metadata when two surfaces are present

[dcesari]

Nella versione attuale dello scanner, quando il livello verticale in grib2 è costituito da due superfici (layer) in cui la prima ha solo il tipo di livello con valore mancante (es. superficie terrestre), il metadato Level: è sbagliato e induce una serie di errori, ad es in arkiweb/meteohub.

Un goffo tentativo di modificare a mano /etc/arkimet/scan/grib2.py:

@@ -76,7 +76,7 @@
     if "typeOfSecondFixedSurface" in grib and grib["typeOfSecondFixedSurface"] != 255:
         level = {
             "style": "GRIB2D",
-            "l1": ltype1, "scale1": lscale1, "value1": lvalue1,
+            "l1": ltype1, "scale": lscale1, "value": lvalue1,
             "l2": grib["typeOfSecondFixedSurface"],
         }
         if grib.get_long("scaleFactorOfSecondFixedSurface") != -1:

sistema il caso di cui sopra, ma rovina il caso in cui la prima superficie ha un livello con valori presenti.

Allego un file level_defekt.zip grib con le varie casistiche incontrate.

arki-scan --yaml grib:level_defekt.grib |grep Level:
# risultato attuale:
Level: GRIB2D(150, 000, 0000000066, 101, 002, 0000000000)
Level: GRIB2D(001,   -, 2147483647, 162, 000, 0000000000)
# risultato con la modifica di cui sopra:
Level: GRIB2D(150,   -,          -,101, 002, 0000000000)
Level: GRIB2D(001,   -,          -,162, 000, 0000000000)
# risultato atteso:
Level: GRIB2D(150, 000, 0000000066, 101, 002, 0000000000)
Level: GRIB2D(001,   -,          -,162, 000, 0000000000)

[dcesari]

Preciso che questi messaggi non sono sempre codificati come si deve, quando il tipo di livello non prevede un valore, i valori dei livelli sono a volte mancanti, a volte presenti (tipicamente =0), ma non è compito di arkimet speculare su questi casi, basta che rispetti i valori scritti nel messaggio, anche se inconsistenti.

[dcesari]

Ciao, sollecito un feedback su questo punto perché è importante per il futuro a breve termine.

[edigiacomo]

Mi viene da dire che il problema non sia nello scanner ma nella definizione del missing nel caso del valore del livello GRIB2, che è definito come GRIB2_MISSING_VALUE = 0xffffffff (vedi https://github.com/ARPA-SIMC/arkimet/blob/v1.47-1/arki/types/level.cc#L34) cioè il massimo di uint32_t (che equivale a 4294967295) ma, nel caso del secondo messaggio GRIB allegato da Davide, ha un valore di 2147483647 = 0x7fffffff.

Un grib_dump del file conferma che viene interpretato come missing:

$ grib_dump level_defekt.grib | grep scaledValueOfFirstFixedSurface
  scaledValueOfFirstFixedSurface = 66;
  scaledValueOfFirstFixedSurface = MISSING;

Guando i sorgenti di eccodes (2.27.1, la versione per Fedora 36), vedo che GRIB_MISSING_LONG ha una storia travagliata e che adesso coincide con 2147483647 (da https://github.com/ecmwf/eccodes/blob/2.27.1/src/grib_api.h#L95-L96):

/* #define GRIB_MISSING_LONG   0x80000001*/
/* #define GRIB_MISSING_LONG 0xffffffff */
#define GRIB_MISSING_LONG 2147483647

E pare che il valore sia passato da 0xffffffff a 2147483647 circa 8 anni fa (https://github.com/ecmwf/eccodes/blame/2.27.1/src/grib_api.h#L95-L96)

Possibile che la definizione del GRIB2_MISSING_VALUE sia stata fatta ai tempi in cui era 0xffffffff? E soprattutto, se definissimo GRIB2_MISSING_VALUE = GRIB_MISSING_LONG avremmo problemi di retrocompatibilità?

[edigiacomo]

Con la seguente modifica:

diff --git a/arki/types/level.cc b/arki/types/level.cc
index 4e99acee..c69762f3 100644
--- a/arki/types/level.cc
+++ b/arki/types/level.cc
@@ -31,7 +31,7 @@ const size_t traits<Level>::type_sersize_bytes = SERSIZELEN;
 
 const uint8_t Level::GRIB2_MISSING_TYPE = 0xff;
 const uint8_t Level::GRIB2_MISSING_SCALE = 0xff;
-const uint32_t Level::GRIB2_MISSING_VALUE = 0xffffffff;
+const uint32_t Level::GRIB2_MISSING_VALUE = GRIB_MISSING_LONG;

Il pacchetto compila, non ci sono errori nei test e il comportamento sembra quello desiderato:

# arki-scan --yaml grib:level_defekt.grib |grep Level:
Level: GRIB2D(150, 000, 0000000066, 101, 002, 0000000000)
Level: GRIB2D(001,   -,          -,162, 000, 0000000000)

È solo una prova per verificare la mia ipotesi, lascio a @spanezz (che ringrazio in anticipo) l'implementazione di una patch più appropriata che tenga eventualmente conto della variazione del valore missing.

[dcesari]

Aggiungo qualche ragionamento, che spero aiuti a fare chiarezza.
Il WMO purtroppo su questo punto fa orecchie da mercante, inizialmente, nel famoso manuale 306, scrive genericamente:

FM 92 GRIB
I.2 – GRIB Reg — 3
REGULATIONS:
...
92.1.4 All bits set to “1” for any value indicates that value is missing. This rule shall not apply to
packed data.
92.1.5 If applicable, negative values shall be indicated by setting the most significant bit to “1”.

Poi però, a quanto ho potuto capire, non specifica quale valore di quale code table può avere segno negativo e quale no (cioè quando si materializza quell'"if applicable").
D'altra parte tutti i bit a 1 per un signed significa -1 che tipicamente è un valore utile, per cui la regola del missing è incoerente.
Nel caso specifico del livello verticale, tabella 4.5, si intuisce uno sforzo per far sì che i valori utili siano >= 0 (se sono sopra il suolo mi propongono il livello altezza sul suolo, se sono nel suolo o nel mare c'è il livello profondità...); alla fine dei conti l'unico caso che mi viene in mente è il definire un livello in atmosfera ad un'altezza rispetto alla superficie del mare che tenga anche conto delle depressioni geografiche, che comunque non è un caso terribilmente esotico, per cui secondo me il livello dovrebbe essere "signed". Evidentemente ad ECMWF hanno fatto lo stesso ragionamento e ad un certo punto hanno cambiato il missing da 0xffffffff a 0x7fffffff intepretando la regola del missing in maniera ragionevole ma difforme da quanto scritto dal WMO ("il massimo numero rappresentabile" e non "tutti i bit a 1").

[spanezz]

È solo una prova per verificare la mia ipotesi, lascio a @spanezz (che ringrazio in anticipo) l'implementazione di una patch più appropriata che tenga eventualmente conto della variazione del valore missing.

Io la ufficializzerei cosí. Mi piace l'idea di usare la costante fornita da eccodes, nella speranza che voglia dire delegare ad eccodes l'interpretazione del valore missing presente nel grib

[spanezz]

Uhm, però, visto che questa roba finisce anche nei metadati arkimet, può non essere cosí semplice. Abbiamo in archivio dei file con dei livelli mancanti?

[dcesari]

Abbiamo ancora pochi dataset grib2 permanenti, erg5 mi pare abbia tutti i livelli in ordine, cosmo 2I ensemble no, da arkiweb vedo questo:

GRIB2D,1,null,2147483647,101,0,0 | sfc Ground or water surface - 2147483647 sfc Mean sea level 0 0
-- | --
GRIB2D,1,null,2147483647,162,0,0 | Formatting failed: /etc/arkimet/format/10-level:51: attempt to concatenate a nil value
GRIB2D,150,0,1,101,2,2200000 | Formatting failed: /etc/arkimet/format/10-level:51: attempt to concatenate a nil value
GRIB2D,150,0,2,101,2,2104000 | Formatting failed: /etc/arkimet/format/10-level:51: attempt to concatenate a nil value

anche se l'errore Formatting failed:... forse è più legato ad arkiweb stesso. Quindi ci sarebbe un archivio da sanare, arkiweb diche che sono, ehm, 19 TB, @brancomat che ne dici? D'altra parte non escludo che i campi affetti dall'errore siano campi che nessuno andrà mai ad estrarre e possiamo lasciarli indicizzati così.

[edigiacomo]

D'altra parte non escludo che i campi affetti dall'errore siano campi che nessuno andrà mai ad estrarre e possiamo lasciarli indicizzati così.

Dimmi che un archivio si può eliminare senza dirmi che un archivio si può eliminare 😄

[spanezz]

che archivio è, cosí vado a pescare qualche file da aggiungere ai test?

[dcesari]

Il dataset è cosmo_2I_fcens (_arc). @edigiacomo intendevo dire che i campi con i livelli bacati potrebbero essere campi di poco interesse (l'altezza geometrica dei livelli del modello). Spero che gli altri campi a qualcuno servano, altrimenti anch'io qua servo a poco, poi sono d'accordo che gli accessi storici sono pochi, ma sto andando fuori tema github.

[spanezz]

Ok, ci sto ragionando un po' piú approfonditamente. level.cc è incentrato sulla codifica binaria dei metadati arkimet, ed è meglio che non cambi perché altrimenti rischiamo di rompere la compatibilità coi dataset già esistenti. Lato metadati arkimet, i livelli usano uint32_t e hanno tutti i bit a 1 per dire "mancante".

I valori dei livelli nei metadati arkimet hanno un fattore moltiplicativo e un valore, ma non un fattore additivo. In GRIB vedo l'incoerenza tra il pezzo che dice "tutti i bit a 1 vuol dire mancante" e il pezzo che dice "il primo bit a 1 vuol dire negativo".

Finora ho lavorato con l'assunto che tutti i valori dei livelli presenti nel GRIB siano positivi, con un bias che viene da BUFR che codifica sempre valori positivi, poi nel caso nella tabella B fornisce un valore di base da sottrarre.

Per il caso in questione, posso modificare lo scanner dei GRIB2 per gestire sia 0xFFFFFFFF che 0x7FFFFFFF come valori mancanti: cosí gestiamo in input entrambe le casistiche, e il formato su disco dei metadati arkimet rimane invariato e coerente. Questo richiede rifare la scansione dei GRIB2 coi livelli mancanti codificati nel modo che arkimet fino ad adesso non si aspetta, ma poi siamo coerenti.

Domanda: abbiamo dei GRIB con dei valori negativi nei livelli?

Al momento i metadati di arkimet non possono rappresentare valori di livelli negativi. C'è da capire se è un caso che non abbiamo ancora visto, o se è un caso che stiamo trattando senza accorgercene. In tutti i modi, se vogliamo gestirlo posso e devo implementare un nuovo stile di codifica livelli che possa salvare anche valori negativi, e farebbe comodo avere dei GRIB di riferimento per guidare lo sviluppo

[spanezz]

Ho fatto un prototipo di fix nella PR #320: i test passano con le aspettative documentate da @dcesari nel primo post di questa issue.

Sono tentato di implementare un arki-check che faccia il rescan dei grib in un dataset e controlli se i metadati risultanti sono uguali a quelli che ci sono nel dataset, in modo da poter validare sia questa modifica nello scanning, che modifiche future

[spanezz]

Ho aggiunto alla PR l'implementazione di arki-maint scantest.

Funziona simile ad arki-query (cioè: arki-maint scantest "query" dataset1 dataset2...) ma quello che fa è querare i dati, riscansionarli, e dire se è cambiato qualcosa.

Con questo dovrebbe essere possibile testare abbastanza ragionevolmente questa e altre possibili future modifiche proposte per la scansione dei dati

[brancomat]

Sulla bontà dell'implementazione proposta in PR passo la palla a @dcesari , per l'archivio da sistemare (i 19Tb di cosmo_2I_fcens_arc), posso lanciare un (lungo) batch quando va in produzione la nuova versione di arkimet, chiedo solo conferma sulle operazioni da fare:

• unarchive dei dati
• arki-check -f
• repack per riarchiviazione
  (la domanda in realtà è: devo prima cancellare gli indici o se ne accorge da solo? ipotizzo la seconda ma non ne sono certo)

[spanezz]

Sulla bontà dell'implementazione proposta in PR passo la palla a @dcesari , per l'archivio da sistemare (i 19Tb di cosmo_2I_fcens_arc), posso lanciare un (lungo) batch quando va in produzione la nuova versione di arkimet, chiedo solo conferma sulle operazioni da fare:

Ci sarebbero da cancellare gli indici, e forse si fa prima a tirar fuori a mano i segmenti grib dagli archivi e reimportarli nel dataset.

Non ho idea di quanto ci possa volere a fare questo discorso su 19Tb (non è forse tanto il numero di tera quanto il numero di GRIB all'interno di quei tera): se si vede che è una cosa eterna per riscansionare una piccola percentuale dei file contenuti, posso provare a implementare in arki-check qualcosa di piú furbo

[dcesari]

Grazie per tutte le discussioni, stavo per scrivere poi ho ulteriormente approfondito qua e là. Scusate se scopro solo adesso queste cose. Ho scoperto con delusione che eccodes considera comunque i valori dei livelli come unsigned /usr/share/eccodes/definitions/grib2/template.4.horizontal.def:

#  Scale factor of first fixed surface
signed[1] scaleFactorOfFirstFixedSurface = missing()  : can_be_missing,dump,no_copy,edition_specific;

#  Scaled value of first fixed surface
unsigned[4] scaledValueOfFirstFixedSurface = missing()  : can_be_missing,dump,no_copy,edition_specific;

quindi non si capisce perché abbiano cambiato il valore missing nel tempo, anche se continuo ad essere convinto che avrebbe meteorologicamente senso, ad es. sul Mar Caspio, definire una superficie atmosferica a -10m sul livello del mare.

Ad ogni modo questo ha come conseguenza che con eccodes attuale non riusciremo mai a codificare dei livelli negativi (provato con grib_set -s e si rifiuta), per cui possiamo evitare di gestirli.

L'altra cosa che ho scoperto è che (1.non bisogna mai dare nulla per scontato) 2.il WMO non usa il complemento a 2 per gl'interi, in effetti non è scritto da nessuna parte, per cui un eventuale signed con tutti i bit a 1 va intepretato come il minimo intero rappresentabile e non -1. Ho verificato controllando il dump esadecimale di un grib in corrispondenza di una latitudine negativa ed è così (per eccodes). Quindi non c'è problema ad intepretare entrambi i valori 0xffffffff e 0x7fffffff come mancanti, nel senso che int32_t_WMO 0xffffffff non è -1 ma un valore di nessun interesse, non so se questo può aiutare la gestione della transizione.

Riguardo la patch, non ho capito molto, ma mi fido.

Coinsiderando il volume di dati da sanare e il fatto che ci sarebbero da sanare anche dei dataset su meteohub di Cineca, sarebbe utile avere uno strumento ad hoc per correggere i metadati, se necessario, senza rileggere e decodificare tutti i dati, non so se è possibile.

[dcesari]

Possiamo procedere con la patch e, una volta aggiornato arkimet, provare a riscannerizzare qualche segmento per capire se è fattibile o se serva qualcosa di più furbo.

[brancomat]

Pull request inclusa in arkimet v1.48-1 (attualmente in build in CI)

[brancomat]

arkimet v1.48-1 installato su server interno "rocky8"

[dcesari]

Con la nuova versione funziona; resta da decidere cosa fare con gli archivi indicizzati con le vecchie versioni.