Analisi dati pozzi
Indice
Ricostruzione tramite GIS della superficie piezometrica della falda alluvionale nell’Alta valle del Tevere
6. Analisi dei livelli di falda con approccio stocastico
6.3 Analisi dati pozzi ’91
Tutti i dati inseriti in Grass sono stati analizzati tramite il programma di statistica[RW. N. Venables, D. M. Smith, 2006].
Per prima cosa si è avviato R da dentro Grass. Sono poi stati importati, attraverso la libreria spgrass6, da Grass in R, i punti e i relativi attributi della mappa vettoriale pozzi91_min50m_pianura tramite i comandi sotto specificati (con # si indicano i commenti delle stringhe di codice R):
> library(spgrass)
# carica la libreria spgrass;
> pozzi91_meno
# trasforma i punti del file vettoriale “pozzi91_min50m_pianura” in uno SpatialPointDataframe, cioè una matrice avente come colonne le coordinate dei punti e tutti i gli attributi e come righe il numero dei punti;
> region
# esporta la regione di Grass in R e la assegna alla variabile region;
> grd<-
GridTopology(cellcentre.offset=c(region$west+(region$ewres/2),
region$south+(region$nsres/2)),cellsize=
c(region$ewres, region$nsres),cells.dim=c(region$cols, region$rows))
# crea una griglia utilizzando i dati della regione di Grass. Le celle che andranno a formare la griglia dovranno essere quadrate. Questo si può ottenere o settando all’inizio la regione di Grass, oppure immettendo direttamente da questo punto i valori che impostano le dimensioni della cella, ad esempio:
“cell.dim=c(50,50)”;
> mask_SG (grd,data=list(k=rep(1,region$cols*region$rows))
proj4string=CRS(region$proj4))
# crea uno SpatialGridDataframe, cioè riferisce nello spazio la griglia precedentemente costruita trasformandola in un DataFrame avente un preciso sistema di riferimento. Si prepara così uno spazio atto ad accogliere i dati della predizione;
Si è poi passato all’analisi statistica vera e propria tramite i comandi di R:
> library(gstat)
# carica la libreria gstat
6.3.1 Analisi statistica univariata dei dati relativi ai pozzi 91
La necessità di questa analisi è verificare il comportamento statistico del dataset, verificare ad esempio se esistono valori di outlier (legati ad esempio un qualche errore grossolano che potrebbe essere stato commesso) o verificare se i dati appartengono ad una stessa popolazione o meno.
> hist(pozzi91_meno$LEV_MAG91,nclass=20)
# plotta l’istogramma dei dati del Maggio ’91 dividendo il range in 20 classi (Figura 6.10), LEV_MAG91 indica il nome della colonna che contiene i dati di livello piezometrico;
> plot(density(pozzi91_meno$LEV_MAG91))
# plotta la curva densità di frequenza del Maggio ’91 (Figura 6.9), secondo un kernel gaussiano
Questi due grafici dicono che il dataset presenta una sola popolazione, poiché presenta un unico massimo, si può quindi affermare quasi con certezza che i dati si riferisco ad un’unica falda.
Osservato che esiste un’unica popolazione andiamo a verificare la presenza di eventuali outliers, cioè quei valori molto lontano dalla media che ricadono prima del 5° percentile o dopo 95° percentile, tramite la costruzione del boxplot del dataset.
> boxplot(pozzi91_meno$LEV_MAG91)
Dal boxplot (Figura 6.11)si può vedere che non esistono outlier nel nostro set di dati in quanto non si ha nessun punto che sta oltre il 5° (linea orizzontale più in basso) o il 95° (linea orizzontale più in alto) percentile e che la media dei valori è a circa 300 metri.
La fase successiva è quella relativa alla creazione dei variogrammi tramite i comandi di R:
>
plot(variogram(LEV_MAG91~1,loc=pozzi91_meno,cutoff=8000,width=800)
)
# plotta il variogramma del livello piezometrico del Maggio ’91(Figura 6.12), LEV_MAG91 è la variabile da stimare, ~1 indica il kriging ordinario (cioè si assume il dato grezzo senza eseguire una detrendizzazione), loc localizza i dati da trattare, cutoff è la distanza massima entro la quale si prendono i dati da interpolare, width è il passo con qui vengono presi i dati (lag)
Dal variogramma si vede chiaramente che la serie di dati presenta un drift (aumenta indefinitamente), ciò è dovuto essenzialmente al fatto che le quote della superficie piezometrica seguono un trend e non risultano casuali, diminuendo da NW verso SE. Bisogna quindi detrendizzare il dato. Per fare questo, prendiamo come dato da studiare la profondità della falda dal piano campagna, in questo modo eliminiamo la parte deterministica del dato ed andiamo a studiare solo la parte stocastica (residuo
Y(x))[vedi 6.2.2].
>
plot(variogram(LEV_MAG91~DEM,loc=pozzi91_meno,cutoff=8000,width=80
0))
# sostituendo ~1 con ~DEM passiamo da kriging ordinario al kriging universale, la colonna DEM contiene le quote del piano campagna nei punti dove sono ubicati i pozzi. Il kriging universale mette in relazione i livello piezometrico con il livello del piano campagna andando ad estrapolare la parte aleatoria del dato (si studiano così solo i residui);
Il variogramma ottenuto mostra ora chiaramente (Figura 6.13) che la variabile residuo, si comporta in modo stazionario, assumendo il tipico andamento dei variogrammi, cioè si stabilizza ad una certa distanza alla quale non si ha più correlazione spaziale.
Una volta isolato il residuo si sono cercate eventuali anisotropie del dato facendo stampare a video la mappa del variogramma tramite il comando:
>
plot(variogram(LEV_MAG91~DEM,loc=pozzi91_meno,cutoff=8000,
width=800,map=T))
# plotta la mappa del variogramma per il dataset del Maggio ’91 aggiungendo il comando map=T (Figura 6.14);
Da questa mappa risulta quasi immediato affermare che esiste un’anisotropia in direzione 135° azimutali. Osservando infatti la mappa possiamo notare una maggiore continuità di colore rosa scuro ad indicare un basso valore di varianza per questa direzione piuttosto che nella direzione ortogonale.
Individuata l’anisotropia questa viene inserita nel calcolo del variogramma direzionale:
>plot(variogram(LEV_MAG91~DEM,loc=pozzi91_meno,cutoff=8000,
width=80 0,alpha=c(45,135)))
# plotta il variogramma direzionale secondo le direzioni 135° e 45° con il comando alpha=c(45,135), la tolleranza angolare è data di default ed pari a 45° cioè la massima apertura per cui i due angoli non si sovrappongono (Figura 6.15);
Si osserva come la forma del variogramma abbia meno oscillazioni ad elevati lag per 135° piuttosto che per 45°, questo perché nella direzione 135° si trovano più punti a distanze elevate piuttosto che nella direzione ortogonale.
Dall’osservazione del variogramma possiamo osservare come esista tra le due direzioni una differenza di range (anisotropia geometrica). Per la direzione 135° il range è ~3000, mentre per la direzione ortogonale vale ~2000. Il rapporto tra questi range è detto rapporto di anisotropia.
6.3.3 Modellazione variogramma
Il passo successivo è quello di modellare il variogramma attraverso delle curve che interpolino i punti discreti (fitting). In base all’andamento generale del variogramma i modelli che intuitivamente meglio lo possono approssimare sono quello sferico o quello esponenziale.
Il modello che deve essere scelto è quello che fornisce gli errori più contenuti nel processo di cross validation del kriging. Nel nostro caso il processo di fitting è stato svolto sia ad occhio che dal calcolatore ed in base ai valori ottenuti si è poi scelto il modello migliore.
>
plot(variogram(LEV_MAG91~DEM,loc=pozzi91_meno,cutoff=8000,width=80
0,alpha=c(135,45)),vgm(model=”Exp”,psill=22,nugget=4,range=3000,an
is=c(135,0.7)))
#plotta il variogramma (variogram()) e il modello fittato manualmente tramite il comando vgm(), “Exp” indica il modello esponenziale, psill il partial sill, nugget il nugget, range il range, anis=c(135,0.7) indica l’anisotropia con la direzione di massima continuità 135 e il rapporto di anisotropia 0.7 tra le due direzioni;
>
plot(variogram(LEV_MAG91~DEM,loc=pozzi91_meno,cutoff=8000,width=80
0,alpha=c(135,45)),fit.variogram(variogram(LEV_MAG91~DEM,loc=pozzi
91_meno,cutoff=8000,width=800,alpha=c(135,45)),vgm(model=”Exp”,psi
ll=36,nugget=5,range=3000,anis=c(135,0.7))))
#plotta il variogramma e il modello fittato dal calcolatore, fit.variogram() è il comando che effettua il fitting automatico;
Lo stesso procedimento è stato svolto con il modello sferico. In questo caso si è cambiato model=”Exp” con model=”Sph” e i parametri di psill, range e nugget nel fitting manuale.
>
plot(variogram(LEV_MAG91~DEM,loc=pozzi91_meno,cutoff=8000,width=80
0,alpha=c(135,45)),vgm(model=”Sph”,psill=18,nugget=4,range=6000,an
is=c(135,0.7)))
>
plot(variogram(LEV_MAG91~DEM,loc=pozzi91_meno,cutoff=8000,width=80
0,alpha=c(135,45)),fit.variogram(variogram(LEV_MAG91~DEM,loc=pozzi
91_meno,cutoff=8000,width=800,alpha=c(135,45)),vgm(model=”Sph”,psi
ll=36,nugget=5,range=3000,anis=c(135,0.7))))
Dopo aver determinato i 4 modelli di variogramma si è passato allo studio sugli errori per i quattro modelli tramite la cross-validation. Per fare la cross-validation salviamo prima i modelli su appositi files.
>
fit_variogramma_mag91_exp=fit.variogram(variogram(LEV_MAG91~DEM,lo
c=pozzi91_meno,cutoff=8000,width=800,alpha=c(135,45)),vgm(model=”E
xp”,psill=1,range=1000,nugget=1,anis=c(135,0.6)))
#salva nel file fit_variogramma_mag91_exp il fitting automatico del variogramma del Maggio ’91 con modello esponenziale;
>
variogramma_mag91_exp=vgm(model=”Exp”,psill=22,nugget=4
range=3000,anis=c(135,0.7))
#salva nel file variogramma_mag91_exp il fitting manuale del variogramma del Maggio ’91 con modello esponenziale;
>
fit_variogramma_mag91_sph=fit.variogram(variogram(LEV_MAG91~DEM,lo
c=pozzi91_meno,cutoff=8000,width=800,alpha=c(135,45)),vgm(model=”S
ph”,psill=36,nugget=5,range=3000,anis=c(135,0.7)))
#salva nel file fit_variogramma_mag91_sph il fitting automatico del variogramma del Maggio ’91 con modello sferico;
>
variogramma_mag91_sph=vgm(model=”Sph”,psill=18,nugget=4,
range=6000,anis=c(135,0.7))
#salva nel file variogramma_mag91_sph il fitting manuale del variogramma del Maggio ’91 con modello sferico.
6.3.4 Cross-validation dei dati
La cross-validation è una parte molto importante del lavoro in quanto permette di determinare quale sia il modello di variogramma migliore per il nostro set di dati. La cross-validation si esegue in R tramite il comando krige.cv.
>
CV_fit_variogramma_mag91_exp=krige.cv(LEV_MAG91~DEM,loc=pozzi91_me
no,model=fit_variogramma_mag91_exp,nfold=nrow(pozzi91_meno))
#salva nel file CV_fit_variogramma_mag91_exp la cross-validation per il modello fit_variogramma_mag91_exp, nfold=nrow(pozzi91_meno) indica che si devono estrarre dal calcolo un pozzo alla volta (LOOCV = Leaveoneout cross validation)
>
CV_variogramma_mag91_exp=krige.cv(LEV_MAG91~DEM,loc=pozzi91_meno
model=variogramma_mag91_exp,nfold=nrow(pozzi91_meno))
#salva nel file CV_variogramma_mag91_exp la cross-validation per il modello variogramma_mag91_exp;
>
CV_fit_variogramma_mag91_sph=
krige.cv(LEV_MAG91~DEM,loc=pozzi91_meno
model=fit_variogramma_mag91_sph,nfold=nrow(pozzi91_meno))
#salva nel file CV_fit_variogramma_mag91_sph la cross-validation per il modello fit_variogramma_mag91_sph;
>
CV_variogramma_mag91_sph=krige.cv(LEV_MAG91~DEM,loc=pozzi91_meno,
model=variogramma_mag91_sph,nfold=nrow(pozzi91_meno))
#salva nel file CV_variogramma_mag91_sph la crossvalidation per il modello variogramma_mag91_sph;
Ognuno dei quattro file così ottenuto è formato da una matrice avente come colonne le coordinate dei punti, la variabile osservata, la variabile stimata, la varianza di stima e i residui (differenza tra variabile osservata e stimata).
Per poter decidere quale sia il modello migliore bisogna ora fare un analisi statistica univariata dei dati stimati dalla cross-validation, più precisamente sui residui.
Per prima cosa si costruiscono gli istogrammi di frequenza dei residui dei quattro modelli (due esponenziali e due sferici) tramite il comando hist() (Figura 6.20 e Figura 6.21).
Le distribuzioni degli errori sono tutte strette attorno allo zero ad indicare una buona stima da parte dell’interpolatore. Presentano anche tutti lo stesso range compreso tra 15 e 10 metri. Essendo gli istogrammi tutti pressoché uguali, si passa ad analizzare i valori degli indici statistici quali media, mediana, primo quartile, terzo quartile, minimo e massimo per i residui, lanciando il comando summary in R e calcolando quindi anche i valori di:
e
> summary(CV_variogramma_mag91_exp$residual)
> summary(CV_fit_variogramma_mag91_exp$residual)
> summary(CV_variogramma_mag91_sph$residual)
> summary(CV_fit_variogramma_mag91_sph$residual)
Dai risultati riassunti in Tabella 2 possiamo vedere che il modello che presenta la maggiore accuratezza, cioè con media dei residui più vicina a zero sia CV_fit_variogramma_mag91_exp, anche se tutti sono molto prossimi a zero.
La precisione è invece determinata dall’errore minimo quadrato (RMSE), quindi il modello più preciso è quello che presenta RMSE più basso, nel nostro caso, CV_fit_variogramma_mag91_sph.
poiché questo modello presenta anche il range più piccolo (minimo 13,74000 e massimo 8,53900), e
la deviazione standard media quadrata (MSDR) più prossima ad 1 è stato possibile affermare che il modello che meglio interpola il variogramma sperimentale del dataset pozzi91 è CV_fit_variogramma_mag91_sph, cioè il modello sferico fittato automaticamente.
Dopo aver accuratamente analizzato i dati della crossvalidation e scelto il modello più appropriato per l’interpolazione, si passa alla costruzione della superficie attraverso il comando krige di R. In particolare si esegue l’universal kriging:
>
UK=krige(LEV_MAG91~DEM,loc=pozzi91_meno,newdata=DEM,
model=fit_variogramma_mag91_sph)
#crea la mappa delle altezze piezometriche nel file UK, newdata indica il data frame che contiene i dati stimati dal kriging, il quale deve contenere la stessa colonna DEM (variabile indipendente) presente nel file pozzi91_meno, con cui viene calcolato il drift.
Completato l’universal kriging la mappa può essere stampata a video attraverso il comando spplot.
> spplot(UK,zcol=1,main=”Altezze piezometriche Maggio ’91”)
ottenendo la mappa in Figura 6.22.
Questa è stata infine importata in Grass attraverso il comando:
> writeRast6sp(UK,”UK”,zcol=”var1.pred”)
La mappa UK importata in Grass (Figura 6.23) è stata in seguito sottratta dalla mappa dtm25 del DEM a 25 metri ottenendo così la mappa delle profondità della falda dal piano campagna (Figura 6.24) attraverso r.mapcalc.
> r.mapcalc “prof_falda_mag91=dtm25UK”
La mappa prodotta evidenzia delle soggiacenze variabili tra 1 e 7 metri. Più precisamente si osserva che la falda ha una profondità di circa 6 metri nella zona sud orientale, vicino Città di Castello e nella zona nord occidentale, poco sotto lo sbarramento di Montedoglio. Profondità maggiori, fino a 50 metri, si creano nella parte orientale (zona azzurra), questi sono, tuttavia, errori legati al fatto che in quella zona il territorio comincia ad essere collinare e quindi di scarso interesse per lo studio relativo alla falda alluvionale.
0 commenti