Spazializzazione dei dati di temperatura e precipitazione su griglia

La creazione dei dataset di precipitazioni e temperature massime e minime giornaliere su gliglia regolare tridimensionale si fonda sulla tecnica dell' "optimal interpolation". 

 Questa tecnica soddisfa la necessità di effettuare un processo di interpolazione spaziale, su griglia regolare, dei dati a disposizione, derivanti da stazioni di rilevamento appartenenti a differenti reti di misurazione, sia con diversa copertura temporale sia con una distribuzione spaziale arbitraria e variabile sul territorio regionale, sfruttando il principio che solo le osservazioni più vicine possono influenzare il valore finale sul punto griglia.

Attraverso tale approccio si rinuncia ad avere un'informazione dettagliata sul singolo punto-stazione ottenendo un campo di analisi su griglia, completo del maggior numero possibile di informazioni del passato, omogeneo ed integrabile con altri dati eventualmente disponibili in futuro e che può essere utilizzato per analisi statistiche sulla variazione relativa nel tempo dei parametri meteo-climatici.

 Il campo di analisi, ricavato con l'optimal interpolation x^asui nodi di una griglia regolare, è ottenuto ponendo in relazione lineare l'incremento dell'analisi (x^a-x^b, dove x^b è il valore del campo di background in corrispondenza degli stessi nodi della griglia) ed il vettore innovazione (y^o-y^b) (Kalnay 1  e Uboldi 2 )

x^a = x^b + K (y^o - y^b)

dove y^o rappresenta il vettore dei valori della variabile meteorologica misurati dalle stazioni e y^b è il valore che il campo di background assume sui punti corrispondenti alle stazioni stesse. Minimizzando la varianza dell'errore dell'analisi si ottiene la definizione della matrice guadagno K:

K = G (S + O) ^-1

dove G è la matrice delle covarianze tra l'errore del campo di background sul punto griglia e l'errore del campo di background sul punto stazione, S è la matrice delle covarianze dell'errore del campo di background su coppie di punti stazione e O è la matrice delle covarianze dell'errore sull'osservazione.

I dati osservati e la scelta del passo griglia

I dati osservati di temperatura e precipitazione, implementati nel dataset, derivano da due reti di misurazione che differiscono non solo per il numero di stazioni e per il periodo di funzionamento, ma presentano anche caratteristiche differenti nei sensori e nella loro gestione:

Rete di stazioni meccaniche e manuali SIMN (disponibile dal 1913 al 2002)Rete di telemisura di Arpa Piemonte (disponibile a partire dal 1988).

Sono stati realizzati controlli di qualità per verificare la fruibilità di tali dati, sottoponendoli preventivamente a verifiche di tipo logico (che garantissero la coerenza del dato ed escludessero valori registrati fisicamente impossibili), di tipo climatologico (tramite i quali sono stati confrontati i singoli valori delle serie storiche con la climatologia delle stazioni stesse ed eventualmente con quelle vicine) e di tipo temporale (onde evitare persistenze anomale di un singolo valore).

Una prima analisi esplorativa è stata effettuata per valutare l'andamento della quantità di dati disponibili, consultando le anagrafiche dei pluviometri e dei termometri e ricavando, anno per anno, il numero di stazioni attive. Dalla figura 1 si evince infatti come il numero di sensori termometrici sia aumentato di oltre dieci volte nell'arco di 50 anni, partendo da circa 25 unità negli anni '50 fino ad arrivare a 371 sensori nel 2009, con un ripido incremento a partire dagli anni '90 in seguito all'introduzione della rete di telemisura. Per quanto riguarda le stazioni pluviometriche, si registra un aumento dei sensori da 119 unità negli anni '50 fino a 386 unità nel 2009.

• 1 Kalnay, E. - Atmospheric modeling, data assimilation and predictability. Cambridge Univ. Press, 341 pp.(2003)
• 2Uboldi, F., Lussana, C. and Salvati, M. - Three-dimensional spatial interpolation of surface meteorological observations from high-resolution local networks. Meteorological Applications (2008)

L'analisi del numero e della distribuzione spaziale delle stazioni è fondamentale anche per stimare un passo adeguato della griglia sulla quale rappresentare il campo di temperatura e di precipitazione ottenuto con l'optimal interpolation. A tal fine è stata effettuata un'analisi in ambiente GIS per ottenere una stima della media del raggio di influenza di ciascuna stazione. E' stata quindi realizzata una tassellatura di Voronoi, attuando una decomposizione del territorio sulla base della distanza di ciascun punto dello spazio bidimensionale rispetto all'insieme delle localizzazioni delle stazioni di rilevamento. L'insieme dei punti più vicini a una determinata stazione rispetto a tutte le altre è definito come la "cella di Voronoi" di quella stazione e costituisce l'unità di base della tassellatura. Per ciascun anno compreso nel periodo 1957-2009 è stata effettuata una tassellatura sulla base delle stazioni attive e sono state calcolate le aree di ciascuna cella. La media di queste aree, insieme alla media del loro raggio (ottenuto approssimando le celle a delle superfici circolari), è indicativa della distanza media tra le stazioni. Con questo metodo è stato possibile ottenere una stima finale non influenzata da valori estremi a causa dell'aumento delle stazioni presenti, come avverrebbe invece considerando la media delle distanze di ciascuna stazione rispetto a tutte le altre. Non è inoltre necessario individuare a priori un numero soggettivo di k vicini, rispetto ai quali calcolare la distanza di ciascuna stazione, come avverrebbe in analisi basate su algoritmi nearest neighbour delle semplici distanze tra le stazioni. Inoltre questo metodo ha il vantaggio di accompagnarsi ad una rappresentazione cartografica che può costituire uno strumento di analisi e di supporto per la localizzazione di criticità locali nel numero di stazioni presenti. I risultati ottenuti sono mostrati in figura di seguito. 

Al termine di questa analisi preliminare, la risoluzione di 15 km è stata individuata come meglio rappresentativa dell'intero periodo preso in esame, 1957-2009.

L'altezza topografica di ciascuna cella della griglia è stata ottenuta con un up-scaling alla risoluzione di 15 km a partire da un DEM di risoluzione 100 m, applicando una tecnica di interpolazione bilineare. Successivamente, ove possibile, i valori della quota di ciascun punto griglia sono stati sostituiti dal valor medio delle quote delle stazioni reali presenti nel momento di massima densità della rete.

Il campo di background

Come precedentemente citato, nel processo di interpolazione l'analisi è calcolata sulla base di un campo di appoggio denominato campo di background. Nel caso delle temperature è stato possibile utilizzare un campo modellistico, ottenuto mediante un downscaling del set di reanalisi ERA-40 elaborato dall'ECMWF – European Centre for Medium range Weather Forecast (disponibile a risoluzione 1.125°) per il periodo 1957-2001 e del set di analisi dell'ECMWF (disponibile a risoluzione 0.500° fino al 2006 e successivamente a risoluzione 0.250°) per il periodo 2002-2009. Preventivamente si è verificato che i segnali climatologici principali, forniti da ERA-40 per il territorio regionale sulle temperature, fossero congruenti con quelli derivanti da un sottoinsieme omogeneo, in termini di tipologia di sensore e, laddove possibile, di densità spaziale, di stazioni presenti sul territorio nel periodo tra il 1950 e il 2000 (Ciccarelli). Nel caso della precipitazione non è stato possibile usufruire di un campo esterno ed il background è stato quindi calcolato a partire dalle osservazioni stesse, mediante una procedura di detrending (Uboldi  2 ).

Il campo IDI e la definizione dei parametri di scala

Uno dei principali obiettivi del lavoro svolto è stato quello di mantenere la matrice guadagno K dell'algoritmo dell'optimal interpolation il più costante possibile, al fine di compensare la grande variabilità nella densità delle stazioni di rilevamento sul territorio nell'intero periodo in esame. Mantenere la matrice guadagno costante permette, infatti, di conservare il più possibile l'omogeneità temporale del segnale e garantisce che il dataset finale sia affetto solo marginalmente dalla non uniformità della distribuzione dei dati.

E' stata quindi svolta un'analisi approfondita del campo IDI, Integral Data Influence, cioè l'analisi dei risultati che si ottengono ponendo i valori del campo di background pari a zero ed il valore delle osservazioni pari a uno. L'equazione si riduce a:

x^a = K

K= G(S + O) ^-1

Assumendo che le correlazioni dell'errore di background nelle matrici G e S possano essere separate nel prodotto tra la correlazione orizzontale e quella verticale (Kalnay) e utilizzando funzioni gaussiane per definire le correlazioni in funzione della sola distanza tra le stazioni ed i punti griglia, la matrice guadagno dipende fondamentalmente dai parametri di scala σ_oriz e σ_vert, oltre che da ε², che è legato alla matrice di covarianza degli errori delle osservazioni e che stima la confidenza nelle osservazioni rispetto al campo di background:

G(ii, jj) = exp (- d²_ijoriz / 2 σ²_oriz) * exp ( - d²_ijvert/2σ²_vert)

dove dij rappresenta la distanza tra il punto griglia i-esimo e la stazione j-esima

S(ii, jj) = exp (- d²_ijoriz / 2 σ²_oriz) * exp ( - d²_ijvert/2σ²_vert)

dove dij rappresenta la distanza tra le coppie di stazioni i-esima e j-esima

O=  ε² I

dove I è la matrice identità.

Ne consegue che, se le funzioni di correlazione si approssimano a zero al crescere della distanza dal punto stazione, il campo IDI sarà approssimativamente pari a 1 nelle porzioni di territorio prossime alle stazioni mentre tenderà a zero nelle zone prive di sensori, e l'area di "influenza" di ciascuna stazione sarà proporzionale al valore dei parametri di scala σ_oriz e σ_vert.

• 2Uboldi, F., Lussana, C. and Salvati, M. - Three-dimensional spatial interpolation of surface meteorological observations from high-resolution local networks. Meteorological Applications (2008)

La scelta dei parametri di scala rappresenta dunque un importante elemento per mascherare o amplificare il contributo totale delle stazioni presenti sul territorio e filtrare la disomogeneità del segnale dovuta al numero variabile di sensori.

Per ogni anno sono stati quindi modulati opportunamente i parametri di scala in funzione del numero di stazioni, per mantenere il campo IDI prossimo ad un valore di riferimento individuato pari a 0.8. Si tenga presente che incrementare notevolmente ed arbitrariamente i parametri di scala non porta conseguentemente ad un valor medio del campo IDI sul dominio pari a 1, poiché l'IDI tende ad un valore asintotico che è strettamente legato al numero delle stazioni. Per questo motivo il valore di riferimento 0.8 è stato definito come il minimo tra i valori asintotici calcolati per ogni anno (vedi figura 5).

Sono stati inizialmente fissati i parametri σvert e ε² pari a valori suggeriti dalla letteratura (Uboldi  2  ) ed è stato variato σoriz in un intervallo compreso tra 10 e 60 km, in modo tale che il valore del corrispondente campo IDI, mediato sull'intero territorio regionale, fosse prossimo al valore di riferimento 0.8. Successivamente è stato ripetuto lo stesso metodo utilizzando il set di σ_oriz precedentemente trovato e variando ε² in un intervallo tra 0.2 e 0.8. Il parametro σ_vert è stato mantenuto fisso a 500 m, in quanto si è verificato operativamente che la sua variazione non influisce significativamente sul campo risultante (figura 6).

• 2Uboldi, F., Lussana, C. and Salvati, M. - Three-dimensional spatial interpolation of surface meteorological observations from high-resolution local networks. Meteorological Applications (2008)

Per quanto riguarda la precipitazione, si è ripetuto il procedimento per la determinazione dei parametri di scala del campo IDI in maniera analoga alla temperatura, variando il termine σoriz da 10 a 30 km, a causa del maggior numero di sensori pluviometrici disponibili nel periodo. Inoltre l'analisi del campo IDI è stata effettuata impostando un campo quota costante, con valore pari a 2m, per annullare nell'algoritmo dell'optimal interpolation la relazione tra il segnale e la quota.

La distribuzione delle precipitazioni in Piemonte evidenzia indubbiamente una relazione con l’altezza topografica del territorio. I flussi umidi provenienti dal Mediterraneo, ai quali è esposto il territorio regionale, interagiscono con la catena montuosa e sono costretti ad un moto ascensionale, determinando una precipitazione di tipo orografico più intensa in corrispondenza della zona di interfaccia con la zona pianeggiante. Tuttavia, dal momento che le diverse aree montane e pedemontane non sono egualmente esposte ai flussi umidi e presentano, a parità di quota, diversa distribuzione spaziale delle precipitazioni anche a causa della differente esposizione dei versanti alpini e della conformazione delle vallate interne, si è preferito eliminare direttamente nell’interpolazione ogni relazione tra precipitazione e quota.

Si è quindi ottenuto, sia per la precipitazione sia per la temperatura, un set di parametri tali da mantenere la matrice guadagno pressoché costante nell'arco temporale studiato e, utilizzando tali parametri nel processo di interpolazione, si è ottenuta la miglior stima possibile del campo di analisi in quelle porzioni di territorio in cui non sono presenti sensori di rilevazione e, al contempo, non è stato sommato arbitrariamente un falso segnale laddove la densità di stazioni varia molto nel tempo.

Una volta ottenuti i campi di temperatura e precipitazione spazializzati sulla griglia a 15 km di risoluzione, sono stati effettuati alcuni controlli per valutare la bontà del processo di interpolazione.

Sono stati quindi confrontate statisticamente alcune serie storiche di stazioni synop della rete GTS, non incluse nel processo di assimilazione, con il corrispondente punto griglia dell'analisi ottenuta con l'optimal interpolation. Nel confronto è stato tenuto conto che, trattandosi di una tecnica di spazializzazione, il punto griglia non riproduce fedelmente i dati della stazione corrispondente in quanto non può tenere conto delle caratteristiche della stazione stessa quali quota, esposizione, ecc. L'obiettivo è stato quindi quello di verificare che nelle due serie storiche a confronto si mantenessero le stesse caratteristiche nelle distribuzioni dei dati.

Come esempio, si mostra nelle figure 8 e 9 l'andamento e la dispersione della temperatura massima e minima della stazione SYNOP di Caselle (quota pari a 287 m) e della temperatura del corrispondente punto griglia (quota pari a 273 m), dopo che le serie storiche giornaliere sono state de-stagionalizzate attraverso il filtro di loess.

Da entrambe le figure si evidenzia come l’andamento della temperatura massima assegnata al punto griglia in corrispondenza della stazione di Caselle segua meglio lo sviluppo del SYNOP corrispondente, rispetto alla temperatura minima. In particolare, il segnale interpolato evidenzia un andamento mediamente superiore all’osservato nel periodo che arriva fino a fine anni ’70, oltre che una dispersione visivamente maggiore nello scatter-plot (figura 9).

Va sottolineato tuttavia che il dato interpolato sulla cella della griglia rappresenta una media elaborata di tutti i punti presenti nella cella stessa, e quindi è ragionevole aspettarsi discrepanze contenute quando la corrispondente serie storica viene confrontata con quella di un punto-stazione singolo, seppur facente parte della cella che lo rappresenta.3

• 3N.Ciccarelli, J. von Hardenberg, A. Provenzale, C. Ronchi, A. Vargiu, R. Pelosini , Climate variability in north-western Italy during the second half of the 20th century. Global and Planetary Change, (2008)