Geoestadística con R. Parte 3. Liberia automap
Daniela Ballari
Septiembre de 2018
Maestría en Ecohidrología
Universidad de Cuenca
http://www.ucuenca.edu.ec/maestria-ecohidrologia/
Daniela Ballari (PhD)
dballari@uazuay.edu.ec
Universidad del Azuay
Publicaciones - https://scholar.google.com/citations?user=1VcAm9AAAAAJ
Objetivo de la práctica. Introducir las funciones básicas de la liberia automap para analisis geoestadistico utilizando el conjunto de datos de “Meuse”. Practica basada en https://sites.ualberta.ca/~ahamann/teaching/renr690/LabKriging2.pdf y en https://cran.r-project.org/web/packages/automap/automap.pdf.
Temario
A- Cargar librerias
B- Cargar datos Meuse
C- Funcion autofitVariogram
D- Funcion autokrige
E- Universal Kriging
F- Funcion autokrige.cv
G- Funcion compare.cv
B-Cargar datos Meuse
Utilizaremos el conjunto de datos “Meuse” de polución del suelo con muestras de concentración de metales pesados. Meuse está previamente cargado en R, de modo que solo lo llamamos con data()
data(meuse)
coordinates(meuse) = ~x+y
str(meuse) #estructura de datos meuse. Algunas variables son continuas y otras categóricas (factors)## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
## ..@ data :'data.frame': 155 obs. of 12 variables:
## .. ..$ cadmium: num [1:155] 11.7 8.6 6.5 2.6 2.8 3 3.2 2.8 2.4 1.6 ...
## .. ..$ copper : num [1:155] 85 81 68 81 48 61 31 29 37 24 ...
## .. ..$ lead : num [1:155] 299 277 199 116 117 137 132 150 133 80 ...
## .. ..$ zinc : num [1:155] 1022 1141 640 257 269 ...
## .. ..$ elev : num [1:155] 7.91 6.98 7.8 7.66 7.48 ...
## .. ..$ dist : num [1:155] 0.00136 0.01222 0.10303 0.19009 0.27709 ...
## .. ..$ om : num [1:155] 13.6 14 13 8 8.7 7.8 9.2 9.5 10.6 6.3 ...
## .. ..$ ffreq : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
## .. ..$ soil : Factor w/ 3 levels "1","2","3": 1 1 1 2 2 2 2 1 1 2 ...
## .. ..$ lime : Factor w/ 2 levels "0","1": 2 2 2 1 1 1 1 1 1 1 ...
## .. ..$ landuse: Factor w/ 15 levels "Aa","Ab","Ag",..: 4 4 4 11 4 11 4 2 2 15 ...
## .. ..$ dist.m : num [1:155] 50 30 150 270 380 470 240 120 240 420 ...
## ..@ coords.nrs : int [1:2] 1 2
## ..@ coords : num [1:155, 1:2] 181072 181025 181165 181298 181307 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:155] "1" "2" "3" "4" ...
## .. .. ..$ : chr [1:2] "x" "y"
## ..@ bbox : num [1:2, 1:2] 178605 329714 181390 333611
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:2] "x" "y"
## .. .. ..$ : chr [1:2] "min" "max"
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
## .. .. ..@ projargs: chr NAC-Funcion autofitVariogram
autoZn=autofitVariogram(log(zinc)~1, meuse)
summary(autoZn)## Experimental variogram:
## np dist gamma dir.hor dir.ver id
## 1 17 59.33470 0.1102869 0 0 var1
## 2 36 86.01449 0.1382850 0 0 var1
## 3 114 131.02870 0.1541601 0 0 var1
## 4 149 176.18845 0.2506480 0 0 var1
## 5 184 226.75652 0.2499933 0 0 var1
## 6 711 337.60359 0.3776705 0 0 var1
## 7 830 502.04773 0.4874784 0 0 var1
## 8 1349 713.21485 0.5928954 0 0 var1
## 9 1314 961.27179 0.6711660 0 0 var1
## 10 1139 1213.41157 0.6374432 0 0 var1
## 11 1355 1506.55052 0.5748673 0 0 var1
##
## Fitted variogram model:
## model psill range
## 1 Nug 0.04848089 0.0000
## 2 Sph 0.58754741 889.9084
## Sums of squares betw. var. model and sample var.[1] 1.434433e-05plot(autoZn, pch=19, col="black")
###D-Funcion autokrige #paletas de colores
bluepal=colorRampPalette(c("azure1", "steelblue4"))
brks=c(4.5,5,5.5,6,6.5,7,7.5,8)
mycol=bluepal(length(brks)-1)
data(meuse.grid)
coordinates(meuse.grid)=~x+y
proj4string(meuse.grid)=CRS("+init=epsg:28992")
meuse.grid=as(meuse.grid,"SpatialPixelsDataFrame")
ZnAKrig=autoKrige(log(zinc)~1, meuse, meuse.grid, model="Sph")## [using ordinary kriging]plot(ZnAKrig, col.regions=mycol)
automapPlot(ZnAKrig$krige_output, "var1.pred",
sp.layout = list("sp.points", meuse))
# Extracting parts from the autoKrige object
prediction_spdf = ZnAKrig$krige_output
sample_variogram = ZnAKrig$exp_var
variogram_model = ZnAKrig$var_modelE-Universal Kriging
plot(log(zinc)~sqrt(dist), meuse)
abline(lm(log(zinc)~sqrt(dist), meuse))
# Auto fit variogram
ZnUKrig=autofitVariogram(log(zinc)~sqrt(dist), meuse)
summary(ZnUKrig)## Experimental variogram:
## np dist gamma dir.hor dir.ver id
## 1 17 59.33470 0.09549033 0 0 var1
## 2 36 86.01449 0.09286878 0 0 var1
## 3 114 131.02870 0.11488929 0 0 var1
## 4 149 176.18845 0.13944591 0 0 var1
## 5 184 226.75652 0.12779452 0 0 var1
## 6 711 337.60359 0.15968833 0 0 var1
## 7 830 502.04773 0.18269393 0 0 var1
## 8 1349 713.21485 0.23735883 0 0 var1
## 9 1314 961.27179 0.24113837 0 0 var1
## 10 1139 1213.41157 0.20799233 0 0 var1
## 11 1355 1506.55052 0.18310788 0 0 var1
##
## Fitted variogram model:
## model psill range
## 1 Nug 0.07963717 0.0000
## 2 Sph 0.14652115 860.3366
## Sums of squares betw. var. model and sample var.[1] 4.655673e-06plot(ZnUKrig, pch=19, col="black")
# Universal Kriging using distance to river
ZnKrig=autoKrige(log(zinc)~sqrt(dist), meuse, meuse.grid, model="Sph")## [using universal kriging]plot(ZnKrig, col.regions=mycol)
F-Funcion autoKrige.cv
kr.cv = autoKrige.cv(log(zinc)~1, meuse, model = c("Exp"), nfold = 10)##
|
| | 0%
|
|======= | 11%
|
|============== | 22%
|
|====================== | 33%
|
|============================= | 44%
|
|==================================== | 56%
|
|=========================================== | 67%
|
|=================================================== | 78%
|
|========================================================== | 89%
|
|=================================================================| 100%kr.cv$krige.cv_output[1:5,]## coordinates var1.pred var1.var observed residual zscore
## 1 (181072, 333611) 6.825357 0.1635548 6.929517 0.10415980 0.25755407
## 2 (181025, 333558) 6.761023 0.1649729 7.039660 0.27863700 0.68601324
## 3 (181165, 333537) 6.256506 0.1879378 6.461468 0.20496169 0.47278708
## 4 (181298, 333484) 6.059823 0.2461048 5.549076 -0.51074685 -1.02954576
## 5 (181307, 333330) 5.581433 0.1735199 5.594711 0.01327875 0.03187738
## fold
## 1 4
## 2 1
## 3 1
## 4 8
## 5 1kr_dist.cv = autoKrige.cv(log(zinc)~sqrt(dist), meuse,
model = c("Exp"), nfold = 10)##
|
| | 0%
|
|======= | 11%
|
|============== | 22%
|
|====================== | 33%
|
|============================= | 44%
|
|==================================== | 56%
|
|=========================================== | 67%
|
|=================================================== | 78%
|
|========================================================== | 89%
|
|=================================================================| 100%kr_dist_ffreq.cv = autoKrige.cv(log(zinc)~sqrt(dist)+ffreq,
meuse, model = c("Exp"), nfold = 10)##
|
| | 0%
|
|======= | 11%
|
|============== | 22%
|
|====================== | 33%
|
|============================= | 44%
|
|==================================== | 56%
|
|=========================================== | 67%
|
|=================================================== | 78%
|
|========================================================== | 89%
|
|=================================================================| 100%G-Funcion compare.cv
compare.cv(kr.cv, kr_dist.cv, kr_dist_ffreq.cv)## kr.cv kr_dist.cv kr_dist_ffreq.cv
## mean_error 0.00871 -0.007814 -0.003188
## me_mean 0.00148 -0.001328 -0.0005417
## MAE 0.3086 0.2724 0.2358
## MSE 0.168 0.1457 0.1025
## MSNE 0.8815 1.075 1.062
## cor_obspred 0.8221 0.8479 0.8956
## cor_predres 0.03085 -0.02737 -0.01492
## RMSE 0.4099 0.3817 0.3201
## RMSE_sd 0.5679 0.5288 0.4434
## URMSE 0.4098 0.3817 0.3201
## iqr 0.4486 0.3798 0.3668compare.cv(kr.cv, kr_dist.cv, kr_dist_ffreq.cv,bubbleplots = TRUE, col.names = c("OK","UK1","UK2"))
## OK UK1 UK2
## mean_error 0.00871 -0.007814 -0.003188
## me_mean 0.00148 -0.001328 -0.0005417
## MAE 0.3086 0.2724 0.2358
## MSE 0.168 0.1457 0.1025
## MSNE 0.8815 1.075 1.062
## cor_obspred 0.8221 0.8479 0.8956
## cor_predres 0.03085 -0.02737 -0.01492
## RMSE 0.4099 0.3817 0.3201
## RMSE_sd 0.5679 0.5288 0.4434
## URMSE 0.4098 0.3817 0.3201
## iqr 0.4486 0.3798 0.3668compare.cv(kr.cv, kr_dist.cv, kr_dist_ffreq.cv,bubbleplots = TRUE, col.names = c("OK","UK1","UK2"),plot.diff = TRUE)
## OK UK1 UK2
## mean_error 0.00871 -0.007814 -0.003188
## me_mean 0.00148 -0.001328 -0.0005417
## MAE 0.3086 0.2724 0.2358
## MSE 0.168 0.1457 0.1025
## MSNE 0.8815 1.075 1.062
## cor_obspred 0.8221 0.8479 0.8956
## cor_predres 0.03085 -0.02737 -0.01492
## RMSE 0.4099 0.3817 0.3201
## RMSE_sd 0.5679 0.5288 0.4434
## URMSE 0.4098 0.3817 0.3201
## iqr 0.4486 0.3798 0.3668