Tecnica geostatistica d'interpolazione dati: il Kriging...

**aleceppi** · 10/06/2008, 16:11

Chi di voi sa usare questa tecnica o ha già fatto dei lavori, magari usando dati riguardanti le precipitazioni?

Ordinary Kriging, Universal Kriging, Block Kriging, ect...?

Per saperne di più vi lascio due articoli: uno sulla distribuzione e ricostruzione del campo pluviometrico in Calabria e uno in Slovenia.

Che cos'è il Kriging...

Da Wikipedia (copia/incollo)...

Il kriging è un metodo di regressione usato nell'ambito dell'analisi spaziale (geostatistica) che permette di interpolare una grandezza nello spazio, minimizzando l’errore quadratico medio. Nell'ambito della statistica è meglio noto come processo gaussiano.
Deve il suo nome a Danie Krige, ingegnere minerario sudafricano che sviluppò negli anni ’50 alcuni metodi empirici per la previsione della distribuzione di minerale nel sottosuolo a partire da campionamenti del terreno. Diversi altri autori si dedicarono alla materia, tra i quali Wold, Kolmogorov, Norbert Wiener, Georges Matheron e Gandin.
Conoscendo il valore di una grandezza in alcuni punti nello spazio (per esempio la temperatura misurata in ogni città di una regione), possiamo determinare il valore della grandezza in altri punti per i quali non esistono misure, per esempio una località di campagna sprovvista di termometri.
Nel kriging, questa interpolazione spaziale si basa sull'autocorrelazione della grandezza, cioè l’assunto che la grandezza in oggetto vari nello spazio con continuità; detto in parole più semplici le cose più vicine sono più simili rispetto alle cose più lontane (Legge di Tobler).
Il valore incognito in un punto viene calcolato con una media pesata dei valori noti.
I pesi che vengono dati alle misure note (cioè alle temperature misurate nelle città) dipendono dalla relazione spaziale tra i valori misurati nell'intorno del punto incognito (cioè il punto in campagna). Per calcolare i pesi si usa il semivariogramma, un grafico che mette in relazione la distanza tra due punti e il valore di semivarianza tra le misure effettuate in questi due punti. Il semivariogramma espone, sia in maniera qualitativa che quantitativa, il grado di dipendenza spaziale, che altro non è che l’autocorrelazione vista prima.

Per saperne di più...
http://mox.polimi.it/~bonavent/tutorials.html

ac

**Ricky nibi** · 10/06/2008, 18:29

grandissimo Ale

Grazie per i pdf, dopo li leggo per bene.

Io ho usato il kriging per fare la mappa degli isospessori della neve al Ghiacciaio del Lupo. Una volta creato il modello mi sono calcolato il volume di acqua presente sul ghiacciaio ad inizio stagione. Ho fatto un po' di test con altri interpolatori. Non ho trovato grandissime differenze nei risultati fra kriging e inverso della distanza. Il kriging però mi permette di inserire una anisotropia e questo mi ha aiutato

se hai altro materiale didattico in merito posta perchè la cosa mi interessa parecchio e vorrei capirci meglio

questo week end saliamo al Lupo con il tuo collega Davide B.

rifacciamo le analisi degli spessori e poi, auspicabilmente, posteremo una mappa di confronto degli spessori 2007 - 2008

**aleceppi** · 10/06/2008, 20:41

Ciao Riky,

sto cercando anch'io del materiale per una sana e robusta bibliografia, poichè devo fare un esame, dove devo ricostruire il campo di precipitazione su tutto il Piemonte e su un'area di pianura piemontese con i dati ARPA che abbiamo al PoliMi.
L'idea è valutare 2 eventi uno convettivo e uno frontale.
Da qui con Matlab e Surfer valutare quale sia la funzione migliore che interpola il semivariogramma, valutare oltre al trend spaziale anche quello temporale (i dati ARPA Piemonte sono ogni 10 minuti), valutare quante stazioni sono necessarie per ricostruire il campo di pioggia, verifica con la cross-validation, controllo dell'anisotropia e infine generare una simulazione Monte Carlo con 50 campi di pioggia (sulla base della varianza dell'errore che ti dà il Kriging) e far girare il modello idrologico avendo in uscita 50 portate.

Ora noi dottorandi avevamo pensato a una cosa più semplice, il professore ci ha dato una legnata mica da ridere. Tant'è siamo in quattro e ci sbatteremo, anche perchè a me questo argomento interessa.

Ti dico solo una cosa, questa mappa in Grads è interpolata con l'algoritmo di Cressman del 1959...

Hai già capito cosa voglio fare....

Tu che programma hai usato per l'interpolazione?

Bavera si dottora martedì, oggi ho visto la sua tesi. C'era anche il tuo nome nei ringraziamenti.
Invece oggi ho fatto il viaggio con Carlo D'Agata in treno.

**Ricky nibi** · 10/06/2008, 21:51

Ho usato surfer

il mio progetto è molto meno ambizioso

questo è un output usando il kriging

gli ho dato un filo di anisotropia per assecondare meglio la morfologia del ghiacciaio

come detto sabato e domenica rifaremo i sondaggi per fare la mappa di confronto. La geostatistica è una figata, soprattutto applicata ai dati climatologici

ringrazierò a mia volta Davide

**Tormenta** · 11/06/2008, 07:54

Si lo conosco, ci ho lavorato spesso.
Molto utile ma anche molto delicato....
La cosa buona è che è un interpolatore esatto, quindi laddove hai il valore ritroverai dopo l'interpolazione il valore che hai messo.
Con Cressman ad esempio puoi avere 25 in un punto e ritrovarti 23 dopo l'interpolazione.

Se è per farci le mappe giornaliere o orarie ti dico che il kriging funziona bene solo se hai tante osservazioni...quindi nei punti della penisola dove ci sono pochi valori, non ha senso usarlo. Servirebbe comunque, anche nel caso tu abbia una buona copertura dati, un autofit del semivariogramma, cosa che lo rende meno preciso. Non puoi infatti metterti li per ogni mappa che vuoi fare a scegliere i parametri (sill,nugget e range).
Ti suggerisco questo bellissimo articolo che mi ha dato molti spunti.
http://ams.allenpress.com/perlserv/?...%2FJAM2295%2E1
Se ti serve il PDF contattami nella chat dello staff.

Ottimo anche il cokriging che ti consente di mantenere il legame con l'orografia, anche se per me è più ostico da usare. Un modello digitale del terreno ad un km di risoluzione lo puoi trovare qui:
http://www.ngdc.noaa.gov/mgg/topo/

Spero di esserti stato di aiuto.
Paolo

**Borat** · 11/06/2008, 09:44

Qualcuno sa se c'è modo di farlo con Stata 9 o con R?

Vorrei elaborare questa cosa qui (www.castionmeteo.it/analisi/analisi.htm) aggiugendo delle mappe. Ho le coordinate delle stazioni, ho le temperature, ma non trovo un cavolo di programma interfacciabile con Stata o R per fare le interpolazioni e creare i grafici.

Ho provato con la funzione IPOLATE di R ma non cavo un ragno dal buco.

Grazie in anticipo.

**Borat** · 11/06/2008, 09:51

PS: avrei bisogno anche che le operazioni fossero totalmente automatizzabili attraverso comandi batch, eventualmente anche con macro o altro.

Basta che si possa automatizzare totalmente, senza dover aprire manualmente il file di dati, definire il formato del grafico, cliccare "esegui" per farglielo plottare, e poi salvarmelo manualmente.

Se Surfer può fare questo, potrei anche prenderlo in considerazione.

Qualcuno ne sa qualcosa?

**Tormenta** · 11/06/2008, 09:56

Originariamente Scritto da Borat

Qualcuno sa se c'è modo di farlo con Stata 9 o con R?

Vorrei elaborare questa cosa qui (www.castionmeteo.it/analisi/analisi.htm) aggiugendo delle mappe. Ho le coordinate delle stazioni, ho le temperature, ma non trovo un cavolo di programma interfacciabile con Stata o R per fare le interpolazioni e creare i grafici.

Ho provato con la funzione IPOLATE di R ma non cavo un ragno dal buco.

Grazie in anticipo.

R fa già tutto, perchè devi trovare una interfaccia?
Fai uno script, lo lanci e hai le mappe.

**atlantic** · 11/06/2008, 10:00

ma che lingua è, aramaico?

**Borat** · 11/06/2008, 10:18

Originariamente Scritto da Tormenta

R fa già tutto, perchè devi trovare una interfaccia?
Fai uno script, lo lanci e hai le mappe.

Uau, tutto in R?

OK, allora la domanda di riserva è (ammetto che andrebbe meglio altrove): che procedura devo usare?

Arrivo a creare un oggetto lista con nomi stazioni, altitudini, coordinate, temperature, ma da lì non so come procedere. A suo tempo ho provato a invocare ipolate ma mi esce qualche errore che non capisco (credo che voglia un dataframe e non una lista e che lo voglia ordinato per coordinate, cose che non mi riescono ancora).

Mi scuso per il contenuto tecnico e i chiari sintomi di dilettantismo totale in R...