Geoestadística 2: precipitación
Daniela Ballari
Universidad de Cuenca
daniela.ballari@ucuenca.edu.ec
Abril, 2015
Objetivo de la práctica. Realizar prácticas con los conceptos básicos de geoestadística aplicado a precipitación y evaluar distintos métodos de interpolación con validación cruzada.
Esta práctica se basa en el documento: “Methdos and Data Sources for Spatial Prediction of Rainfall” de T.Hengel, A. AghaKouchak y M. Pecec Tadie (2000) y en el scrip rainfall_mapping.R http://spatial-analyst.net/book/rainfall_mapping.R
Pasos a seguir
A- Cargar librerias y datos
B- Seleccionar ejemplo de datos diarios
B.1- Estructura espacial local: autocorrelación
B.2- Interpolación con kriging, iDW y Thiessen
B.3- Validación cruzada
C- Seleccionar ejemplo de datos acumulados anuales
C.1- Estructura espacial local: autocorrelación
C.2- Interpolación con kriging, iDW y Thiessen
C.3- Validación cruzada
A-Cargar librerias y datos
rm(list=ls())#borrar todo lo que haya en memoria
library(rgdal)## Loading required package: sp
## rgdal: version: 0.8-16, (SVN revision 498)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.10.1, released 2013/08/26
## Path to GDAL shared files: C:/Users/DOCENTE11/Documents/R/win-library/3.0/rgdal/gdal
## GDAL does not use iconv for recoding strings.
## Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
## Path to PROJ.4 shared files: C:/Users/DOCENTE11/Documents/R/win-library/3.0/rgdal/projlibrary(gstat)
#definir directorio de trabajo
setwd("C:/curso_R")
#Descargar datos y descomprimirlos en el directorio de trabajo
shell.exec("https://drive.google.com/uc?export=download&id=0B0ea6usixQHFaExTZExYREhuaUU")
#Datos
limite <- readOGR(".", "countries")## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "countries"
## with 6 features and 6 fields
## Feature type: wkbPolygon with 2 dimensionsPREC.2008 <- read.table("GSOD_2008_bilogora.csv", sep=";", header=TRUE)
names(PREC.2008)## [1] "STNID" "STN" "WBAN" "YEARMODA"
## [5] "TEMPC" "TEMP.count" "PREC" "PREC.flag"
## [9] "DATE" "STATION.NAME" "LAT" "LON"
## [13] "ELEV..1M."#STN: nombre de estación
#WBAN: identificador del weather bureau air force navy
#TEMPC: temperatura diaria acumulada (celcius)
#TEMP.count: número de observaciones usadas para calcular la temperatura media
#PREC: precipitación total diaria acumulada (mm)
#LAT y LON: coordenadas en latitud longitud, wgs84
#DATE: fecha de observación
str(PREC.2008)## 'data.frame': 38889 obs. of 13 variables:
## $ STNID : Factor w/ 111 levels "110162-99999",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ STN : int 110162 110162 110162 110162 110162 110162 110162 110162 110162 110162 ...
## $ WBAN : int 99999 99999 99999 99999 99999 99999 99999 99999 99999 99999 ...
## $ YEARMODA : int 20080802 20080804 20080806 20080808 20080805 20080810 20080801 20080803 20080809 20080624 ...
## $ TEMPC : num 22.7 27.3 24.7 23.1 21.5 22.8 26.8 25.5 18.4 25 ...
## $ TEMP.count : int 11 11 12 12 12 12 12 12 12 11 ...
## $ PREC : num NA 0 0 NA NA 0 0 0 0 NA ...
## $ PREC.flag : Factor w/ 8 levels "A","B","E","F",..: 8 7 7 8 8 7 7 7 7 8 ...
## $ DATE : Factor w/ 366 levels "2008-01-01","2008-01-02",..: 215 217 219 221 218 223 214 216 222 176 ...
## $ STATION.NAME: Factor w/ 107 levels "AGARD","AIGEN IM ENNSTAL",..: 46 46 46 46 46 46 46 46 46 46 ...
## $ LAT : num 48 48 48 48 48 ...
## $ LON : num 16.3 16.3 16.3 16.3 16.3 ...
## $ ELEV..1M. : int 233 233 233 233 233 233 233 233 233 233 ...B-DATOS DIARIOS: Seleccionar ejemplo de datos diario
# Formatear datos de fechas
PREC.2008$DATE <- as.Date(PREC.2008$DATE)
# Seleccionar datos de un día en específico (1 de mayo 2008) y cuyos valores no sean vacios
PREC.20080501 <- subset(PREC.2008, PREC.2008$DATE==as.Date("2008-05-01")&!is.na(PREC.2008$PREC))
str(PREC.20080501) # 99 estaciones## 'data.frame': 99 obs. of 13 variables:
## $ STNID : Factor w/ 111 levels "110162-99999",..: 3 4 5 6 7 8 9 10 11 12 ...
## $ STN : int 111570 111650 111700 111710 111730 111740 111750 111800 111820 111850 ...
## $ WBAN : int 99999 99999 99999 99999 99999 99999 99999 99999 99999 99999 ...
## $ YEARMODA : int 20080501 20080501 20080501 20080501 20080501 20080501 20080501 20080501 20080501 20080501 ...
## $ TEMPC : num 8.9 8.8 9.7 7.8 7 9.3 8.4 2.4 12.3 7.2 ...
## $ TEMP.count : int 20 20 8 8 7 8 8 8 8 8 ...
## $ PREC : num 1.52 2.03 0 0 7.11 ...
## $ PREC.flag : Factor w/ 8 levels "A","B","E","F",..: 4 4 4 3 3 3 3 3 4 3 ...
## $ DATE : Date, format: "2008-05-01" "2008-05-01" ...
## $ STATION.NAME: Factor w/ 107 levels "AGARD","AIGEN IM ENNSTAL",..: 2 106 47 85 22 87 12 70 97 50 ...
## $ LAT : num 47.5 47.2 47.9 47.8 47.5 ...
## $ LON : num 14.1 14.8 15.1 15.3 15.6 ...
## $ ELEV..1M. : int 649 677 615 872 1050 572 482 1554 285 994 ...head(PREC.20080501)## STNID STN WBAN YEARMODA TEMPC TEMP.count PREC PREC.flag
## 797 111570-99999 111570 99999 20080501 8.9 20 1.52 F
## 1161 111650-99999 111650 99999 20080501 8.8 20 2.03 F
## 1525 111700-99999 111700 99999 20080501 9.7 8 0.00 F
## 1966 111710-99999 111710 99999 20080501 7.8 8 0.00 E
## 2133 111730-99999 111730 99999 20080501 7.0 7 7.11 E
## 2655 111740-99999 111740 99999 20080501 9.3 8 1.02 E
## DATE STATION.NAME LAT LON ELEV..1M.
## 797 2008-05-01 AIGEN IM ENNSTAL 47.533 14.133 649
## 1161 2008-05-01 ZELTWEG 47.200 14.750 677
## 1525 2008-05-01 LUNZ 47.850 15.067 615
## 1966 2008-05-01 ST SEBASTIAN/MARIAZ 47.800 15.300 872
## 2133 2008-05-01 FISHBACH 47.450 15.650 1050
## 2655 2008-05-01 ST. MICHAEL/LEOBEN 47.333 15.000 572#Asignar coordenadas para convertir en Spatial.Points.Data.Frame
coordinates(PREC.20080501) <- ~LON+LAT # recordar que otra opción es c(LON,LAT)
#Asignar sistema de referencia
proj4string(PREC.20080501) <- CRS("+proj=longlat +datum=WGS84")
#Definir sistema de referencia UTM y convertir
utm33 <- "+proj=utm +zone=33 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
PREC.20080501.XY <- spTransform(PREC.20080501, CRS(utm33))
#En los datos hay una estación duplicada, de manera que se eliminará
table(duplicated(PREC.20080501.XY$STATION.NAME))##
## FALSE TRUE
## 98 1PREC.20080501.XY <- remove.duplicates(PREC.20080501.XY)
#Visualizar
PREC.plt1<- bubble(PREC.20080501.XY, "PREC", col="black", pch=21,
main="PREC 2008-05-01", sp.layout=list(list("sp.polygons", col="grey",
fill="transparent", limite), list("sp.points", col="black", pch="+",
cex=1.2, subset(PREC.20080501.XY, PREC.20080501.XY$PREC==0))))
PREC.plt1
dev.off()## null device
## 1B.1- Estructura espacial local: autocorrelación
#Variograma empírico
PREC.ve.d <- variogram(log1p(PREC)~1, PREC.20080501.XY) 
PREC.ve.d## np dist gamma dir.hor dir.ver id
## 1 15 9950.267 0.2559174 0 0 var1
## 2 40 22233.377 0.4779323 0 0 var1
## 3 71 34966.919 0.4814682 0 0 var1
## 4 122 49628.899 0.4875598 0 0 var1
## 5 143 63891.616 0.5043048 0 0 var1
## 6 158 77347.376 0.6710930 0 0 var1
## 7 167 92034.124 0.6796110 0 0 var1
## 8 180 105881.739 0.6618307 0 0 var1
## 9 184 120078.467 0.8772658 0 0 var1
## 10 214 133557.099 0.6997319 0 0 var1
## 11 211 148425.666 0.6473379 0 0 var1
## 12 206 161738.080 0.8461399 0 0 var1
## 13 200 175694.274 0.7072404 0 0 var1
## 14 237 189678.732 0.6737083 0 0 var1
## 15 219 204239.671 0.6712586 0 0 var1# variograma inicial
PREC.vi.d <- vgm(nugget=0, model="Exp",
range=sqrt(diff(PREC.20080501.XY@bbox["LON",])^2 +
diff(PREC.20080501.XY@bbox["LAT",])^2)/4,
psill=var(log1p(PREC.20080501.XY$PREC)))
# Regla general para evitar determinar parámetros visualmente
# Nugget = 0
# partial-sill = varianza total de los datos
# Range = 1/4 de la distancia maxima, diagonal del bounding box
PREC.vi.d## model psill range
## 1 Nug 0.0000000 0.0
## 2 Exp 0.6503578 158537.4#Ajuste del variograma (teórico)
PREC.vt.d <- fit.variogram(PREC.ve.d, model=PREC.vi.d)
PREC.vt.d## model psill range
## 1 Nug 0.1651774 0.0
## 2 Exp 0.5688303 44082.5plot(PREC.ve.d, pl = T, model = PREC.vt.d)
B.2- Interpolación con kriging, iDW y Thiessen
# Preparar gilla para predicción
PREC.20080501.XY_grid = spsample(PREC.20080501.XY, type = "regular", cellsize = c(1000,1000))
#spsample selecciona una muestra regular de puntos en la extensión geográfica de PREC.20080501.XY.
#Los puntos estarán espaciados uno de otro 1000m.
class(PREC.20080501.XY_grid)## [1] "SpatialPoints"
## attr(,"package")
## [1] "sp"gridded(PREC.20080501.XY_grid) = TRUE
class(PREC.20080501.XY_grid)## [1] "SpatialPixels"
## attr(,"package")
## [1] "sp"#Thiessen
thiessen.d = krige(PREC ~ 1, PREC.20080501.XY, PREC.20080501.XY_grid, nmax = 1)## [inverse distance weighted interpolation]pts.s <- list("sp.points", PREC.20080501.XY, col="white",pch=20)
spplot(thiessen.d, "var1.pred", asp=1, at= seq(0,16,1), zlim=c(0,16),
col.regions=gray(seq(0.9,0.1,l=30)),sp.layout = list(pts.s),main="Thiessen")
#Inverso de la distancia (idw)
idw.d = idw(PREC ~ 1, PREC.20080501.XY, PREC.20080501.XY_grid)## [inverse distance weighted interpolation]spplot(idw.d, "var1.pred", asp=1, at= seq(0,16,1), zlim=c(0,16),
col.regions=gray(seq(1,0.1,l=30)),sp.layout = list(pts.s),main="IDW")
# Ordinary kriging
ok.d <- krige(log1p(PREC) ~ 1, locations = PREC.20080501.XY,
newdata = PREC.20080501.XY_grid, model = PREC.vt.d) ## [using ordinary kriging]ok.d$PREC.pred <- expm1(ok.d$var1.pred)# Antilog para volver a valores de precipitación
par(mfrow=c(2,1))
print(spplot(ok.d, "PREC.pred", asp=1, at= seq(0,16,1), zlim=c(0,16),
col.regions=gray(seq(1,0.1,l=30)),main="OK, diaria 01/05/2008",
sp.layout = list(pts.s), zlim=c(0,16)), split=c(1,1,2,1), more=TRUE)
print(spplot(ok.d, "PREC.pred", asp=1, col.regions=gray(seq(1,0.1,l=30)),
main="Varianza OK, diaria 01/05/2008",
sp.layout = list(pts.s)), split=c(2,1,2,1), more=FALSE)
B.3-Validación cruzada
thiessen.cv.d <- krige.cv(PREC ~ 1, PREC.20080501.XY,nmax = 1)
idw.cv.d <- krige.cv(PREC ~ 1, PREC.20080501.XY)
ok.cv.d <- krige.cv(log1p(PREC) ~ 1, locations = PREC.20080501.XY, model = PREC.vt.d)
thiessen.cv.d[1:5,]## coordinates var1.pred var1.var observed residual zscore fold
## 1 (434744.4, 5264761) 0.00 NA 1.52 1.52 NA 1
## 2 (481064.7, 5227420) 1.02 NA 2.03 1.01 NA 2
## 3 (505012.4, 5299631) 0.00 NA 0.00 0.00 NA 3
## 4 (522465, 5294115) 0.00 NA 0.00 0.00 NA 4
## 5 (549000.1, 5255378) 5.08 NA 7.11 2.03 NA 5idw.cv.d[1:5,]## coordinates var1.pred var1.var observed residual zscore fold
## 1 (434744.4, 5264761) 1.738358 NA 1.52 -0.2183581 NA 1
## 2 (481064.7, 5227420) 2.658482 NA 2.03 -0.6284822 NA 2
## 3 (505012.4, 5299631) 2.032641 NA 0.00 -2.0326415 NA 3
## 4 (522465, 5294115) 2.376153 NA 0.00 -2.3761529 NA 4
## 5 (549000.1, 5255378) 4.375667 NA 7.11 2.7343332 NA 5ok.cv.d[1:5,]## coordinates var1.pred var1.var observed residual zscore
## 1 (434744.4, 5264761) 0.3287548 0.5696455 0.9242589 0.5955041 0.7890099
## 2 (481064.7, 5227420) 0.7408164 0.5236523 1.1085626 0.3677462 0.5081907
## 3 (505012.4, 5299631) 0.2126372 0.5191697 0.0000000 -0.2126372 -0.2951105
## 4 (522465, 5294115) 0.4035653 0.4818973 0.0000000 -0.4035653 -0.5813486
## 5 (549000.1, 5255378) 1.6770895 0.4539861 2.0930979 0.4160084 0.6174202
## fold
## 1 1
## 2 2
## 3 3
## 4 4
## 5 5ok.cv.d$obs <- PREC.20080501.XY$PREC
#correlación
par(mfrow=c(1,3))
print(plot(var1.pred~observed,thiessen.cv.d, main="Thiessen"), split=c(1,1,3,1), more=TRUE)## NULLprint(plot(var1.pred~observed,idw.cv.d, main="IDW"), split=c(2,1,3,1), more=TRUE)## NULLprint(plot(var1.pred~obs,ok.cv.d, main="OK"), split=c(3,1,2,1), more=FALSE)
## NULL#correlación, idealmente 1
cor(thiessen.cv.d$var1.pred,thiessen.cv.d$observed)## [1] 0.4627582cor(idw.cv.d$var1.pred,idw.cv.d$observed)## [1] 0.6617524cor(ok.cv.d$var1.pred,ok.cv.d$obs)## [1] 0.7272514#media de residuos, idealmente 0
mean(thiessen.cv.d$residual)## [1] 0.1426531mean(idw.cv.d$residual)## [1] -0.1804113mean(ok.cv.d$residual)## [1] -0.007739089#Desviaciones estandar del error de interpolación (residuos). Idealmente pequeño
sd(thiessen.cv.d$residual)## [1] 3.252934sd(idw.cv.d$residual)## [1] 2.606365sd(ok.cv.d$residual)## [1] 0.528048#Boxplots
boxplot(thiessen.cv.d$residual, idw.cv.d$residual, ok.cv.d$residual, main="Thiessen, IDW, OK")
# MSPE (mean square predictor error), idealmente pequeño
mean(thiessen.cv.d$residual^2)## [1] 10.49396mean(idw.cv.d$residual^2)## [1] 6.756368mean(ok.cv.d$residual^2)## [1] 0.2760494#Error medio cuadrático (RMSE) es una medida general. Idealmente pequeño
sqrt(sum(thiessen.cv.d$residual^2)/length(thiessen.cv.d$residual))## [1] 3.239438sqrt(sum(idw.cv.d$residual^2)/length(idw.cv.d$residual))## [1] 2.599301sqrt(sum(ok.cv.d$residual^2)/length(ok.cv.d$residual))## [1] 0.525404#Varianza de residuos, idealmente cercanos a cero
var(thiessen.cv.d$residual, na.rm=T)## [1] 10.58158var(idw.cv.d$residual, na.rm=T)## [1] 6.793137var(ok.cv.d$residual, na.rm=T)## [1] 0.2788347# Cantidad de variación explicada por el modelo en cada método, idealmente cercano a 100
(1-var(thiessen.cv.d$residual, na.rm=T)/var(PREC.20080501.XY$PREC))*100## [1] 7.4026(1-var(idw.cv.d$residual, na.rm=T)/var(PREC.20080501.XY$PREC))*100## [1] 40.55456(1-var(ok.cv.d$residual, na.rm=T)/var(PREC.20080501.XY$PREC))*100## [1] 97.55997
C-DATOS ANUALES. Seleccionar ejemplo de datos anuales acumulados
# Agregar (sumar, acumular) datos para el año 2008
PREC.2008.a <- aggregate(PREC.2008[!is.na(PREC.2008$PREC),c("PREC")],
by=list(PREC.2008$STNID[!is.na(PREC.2008$PREC)]), FUN=sum)
# Argumento 1: PREC.2008[!is.na(PREC.2008$PREC),c("PREC")]. Objeto sobre el cual se aplicará
# la agregración, se seleccionan las precipitaciones que no sean vacios y las columnas PREC.
# Deben ser valores numéricos
# Argumeno 2: by=list(PREC.2008$STNID[!is.na(PREC.2008$PREC)]). Columna sobre la cual se
# agregará: por cada estación (STNID) se agregarán los valores de PREC que no sean vacíos.
# Argumento 3: operación de agregación (media, suma, etc.). En este caso suma.
# Seleccionar coordenadas de estaciones para mantener georeferenciación.
#Unique elimina estaciones repetidas
PREC.2008.a.lon.lat <- unique(PREC.2008[!is.na(PREC.2008$PREC),c("STNID", "LON", "LAT")])
head(PREC.2008.a)## Group.1 x
## 1 110162-99999 0.00
## 2 110165-99999 0.00
## 3 111570-99999 874.31
## 4 111650-99999 745.99
## 5 111700-99999 1772.18
## 6 111710-99999 870.72head(PREC.2008.a.lon.lat)## STNID LON LAT
## 2 110162-99999 16.267 47.967
## 366 110165-99999 16.250 47.850
## 732 111570-99999 14.133 47.533
## 1097 111650-99999 14.750 47.200
## 1463 111700-99999 15.067 47.850
## 1821 111710-99999 15.300 47.800#Unir ambas variables a traves de la columna Group.1 y STNID
PREC.2008.a.unido<- merge(x = data.frame(PREC.2008.a), y = PREC.2008.a.lon.lat,
by.x = "Group.1", by.y="STNID", all.x=TRUE)
head(PREC.2008.a.unido)## Group.1 x LON LAT
## 1 110162-99999 0.00 16.267 47.967
## 2 110165-99999 0.00 16.250 47.850
## 3 111570-99999 874.31 14.133 47.533
## 4 111650-99999 745.99 14.750 47.200
## 5 111700-99999 1772.18 15.067 47.850
## 6 111710-99999 870.72 15.300 47.800#Asignar coordenadas
coordinates(PREC.2008.a.unido) <- ~LON+LAT
#Asignar sistema de referencia
proj4string(PREC.2008.a.unido) <- CRS("+proj=longlat +datum=WGS84")
#Transformar sistema de referencia a UTM
PREC.2008.a.XY <- spTransform(PREC.2008.a.unido, CRS(utm33))
#Borrar posibles valores duplicados
PREC.2008.a.XY <- remove.duplicates(PREC.2008.a.XY)
names(PREC.2008.a.XY) <- c("estacion", "PREC")
#Reemplazar valores cero por NA
PREC.2008.a.XY$PREC2 <- ifelse(PREC.2008.a.XY$PREC==0, NA, PREC.2008.a.XY$PREC)
#Seleccionar valores no vacios
PREC.2008.a.XY <- subset(PREC.2008.a.XY, !is.na(PREC.2008.a.XY$PREC2))
PREC.plt2 <- bubble(PREC.2008.a.XY, "PREC", col="black", do.sqrt=FALSE, pch=21,
main="PREC acumulada 2008", sp.layout=list("sp.polygons",
col="grey", fill="transparent", limite))
PREC.plt2
C.1- Estructura espacial local: autocorrelación
# Variograma empírico
PREC.ve.a <- variogram(log1p(PREC)~1, PREC.2008.a.XY)
PREC.ve.a## np dist gamma dir.hor dir.ver id
## 1 15 9950.267 0.03115945 0 0 var1
## 2 40 22331.174 0.04758893 0 0 var1
## 3 72 34966.627 0.03075885 0 0 var1
## 4 123 49737.590 0.07971122 0 0 var1
## 5 150 63841.723 0.06145396 0 0 var1
## 6 163 77430.847 0.05734990 0 0 var1
## 7 168 91852.213 0.06072228 0 0 var1
## 8 186 105814.799 0.06773771 0 0 var1
## 9 188 120027.715 0.06763130 0 0 var1
## 10 217 133541.982 0.07837992 0 0 var1
## 11 211 148509.610 0.07543351 0 0 var1
## 12 212 161698.148 0.09966007 0 0 var1
## 13 204 175626.769 0.09803985 0 0 var1
## 14 229 189792.835 0.09778899 0 0 var1
## 15 219 204032.112 0.10342482 0 0 var1# Variograma inicial
PREC.vi.a <- vgm(nugget=0, model="Exp",
range=sqrt(diff(PREC.2008.a.XY@bbox["LON",])^2 +
diff(PREC.2008.a.XY@bbox["LAT",])^2)/4,
psill=var(log1p(PREC.2008.a.XY$PREC)))
# Regla general
# Nugget = 0
# Sill varianza total de los datos
# Range = 1/4 de la distancia maxima, diagonal del bounding box
PREC.vi.a## model psill range
## 1 Nug 0.0000000 0.0
## 2 Exp 0.1131792 158537.4#Ajuste del variograma (teórico)
PREC.vt.a <- fit.variogram(PREC.ve.a, model=PREC.vi.a)
PREC.vt.a## model psill range
## 1 Nug 0.02639123 0.0
## 2 Exp 0.07923009 112468.2plot(PREC.ve.a, pl = T, model = PREC.vt.a)
C.3- Interpolación con kriging, iDW y Thiessen
# Preparar gilla para predicción
PREC.2008.a.XY_grid = spsample(PREC.2008.a.XY, type = "regular", cellsize = c(1000,1000))
gridded(PREC.2008.a.XY_grid) = TRUE
# Thiessen
thiessen.a = krige(PREC ~ 1, PREC.2008.a.XY, PREC.2008.a.XY_grid, nmax = 1)## [inverse distance weighted interpolation]pts.s <- list("sp.points", PREC.2008.a.XY, col="white",pch=20)
spplot(thiessen.a, "var1.pred", asp=1, at= seq(400,2100,50), zlim=c(400,2100),
col.regions=gray(seq(0.9,0.1,l=30)),sp.layout = list(pts.s),main="Thiessen anual")
#Inverso de la distancia (idw)
idw.a = idw(PREC ~ 1, PREC.2008.a.XY, PREC.2008.a.XY_grid)## [inverse distance weighted interpolation]spplot(idw.a, "var1.pred", asp=1, at= seq(400,2100,50), zlim=c(400,2100),
col.regions=gray(seq(1,0.1,l=30)),sp.layout = list(pts.s),main="IDW anual")
# Ordinary kriging
ok.a <- krige(log1p(PREC) ~ 1, locations = PREC.2008.a.XY,
newdata = PREC.2008.a.XY_grid, model = PREC.vt.a) ## [using ordinary kriging]ok.a$PREC.pred <- expm1(ok.a$var1.pred)# Antilog para volver a valores de precipitación
par(mfrow=c(2,1))
print(spplot(ok.a, "PREC.pred", asp=1, at= seq(400,2100,50), zlim=c(400,2100),
col.regions=gray(seq(1,0.1,l=35)),main="OK, anual 2008",
sp.layout = list(pts.s), zlim=c(0,16)), split=c(1,1,2,1), more=TRUE)
print(spplot(ok.a, "var1.var", asp=1, col.regions=gray(seq(1,0.1,l=30)),
main="Varianza OK, anual 2008",sp.layout = list(pts.s)),
split=c(2,1,2,1), more=FALSE)
C.3-Validación cruzada
thiessen.cv.a <- krige.cv(PREC ~ 1, PREC.2008.a.XY,nmax = 1)
idw.cv.a <- krige.cv(PREC ~ 1, PREC.2008.a.XY)
ok.cv.a <- krige.cv(log1p(PREC) ~ 1, locations = PREC.2008.a.XY, model = PREC.vt.a)
ok.cv.a$obs <- PREC.2008.a.XY$PREC
thiessen.cv.a[1:5,]## coordinates var1.pred var1.var observed residual zscore fold
## 1 (434744.4, 5264761) 1158.71 NA 874.31 -284.40 NA 1
## 2 (481064.7, 5227420) 712.76 NA 745.99 33.23 NA 2
## 3 (505012.4, 5299631) 870.72 NA 1772.18 901.46 NA 3
## 4 (522465, 5294115) 1772.18 NA 870.72 -901.46 NA 4
## 5 (549000.1, 5255378) 1005.03 NA 1005.85 0.82 NA 5idw.cv.a[1:5,]## coordinates var1.pred var1.var observed residual zscore fold
## 1 (434744.4, 5264761) 960.5473 NA 874.31 -86.23730 NA 1
## 2 (481064.7, 5227420) 837.4586 NA 745.99 -91.46858 NA 2
## 3 (505012.4, 5299631) 845.1617 NA 1772.18 927.01833 NA 3
## 4 (522465, 5294115) 1090.8605 NA 870.72 -220.14047 NA 4
## 5 (549000.1, 5255378) 835.3164 NA 1005.85 170.53358 NA 5ok.cv.a[1:5,]## coordinates var1.pred var1.var observed residual zscore
## 1 (434744.4, 5264761) 6.911983 0.05847249 6.774578 -0.13740463 -0.5682318
## 2 (481064.7, 5227420) 6.703703 0.05149363 6.616052 -0.08765163 -0.3862631
## 3 (505012.4, 5299631) 6.753763 0.05390556 7.480530 0.72676726 3.1302467
## 4 (522465, 5294115) 7.018383 0.04935535 6.770468 -0.24791499 -1.1159267
## 5 (549000.1, 5255378) 6.667441 0.04538308 6.914582 0.24714046 1.1601039
## fold obs
## 1 1 874.31
## 2 2 745.99
## 3 3 1772.18
## 4 4 870.72
## 5 5 1005.85#correlación
par(mfrow=c(1,3))
print(plot(var1.pred~observed,thiessen.cv.a, main="Thiessen"), split=c(1,1,3,1), more=TRUE)## NULLprint(plot(var1.pred~observed,idw.cv.a, main="IDW"), split=c(2,1,3,1), more=TRUE)## NULLprint(plot(var1.pred~obs,ok.cv.a, main="OK"), split=c(3,1,2,1), more=FALSE)
## NULL#correlación, idealmente 1
cor(thiessen.cv.a$var1.pred,thiessen.cv.a$observed)## [1] 0.6756613cor(idw.cv.a$var1.pred,idw.cv.a$observed)## [1] 0.7475555cor(ok.cv.a$var1.pred,ok.cv.a$obs)## [1] 0.7798622#media de residuos, idealmente 0
mean(thiessen.cv.a$residual)## [1] -17.03296mean(idw.cv.a$residual)## [1] -2.390472mean(ok.cv.a$residual)## [1] 0.0002767632#Desviaciones estandar del error de interpolación (residuos). Idealmente pequeño
sd(thiessen.cv.a$residual)## [1] 269.9929sd(idw.cv.a$residual)## [1] 218.5966sd(ok.cv.a$residual)## [1] 0.1945751#Boxplots
boxplot(thiessen.cv.a$residual, idw.cv.a$residual, ok.cv.a$residual, main="Thiessen, IDW, OK")
# MSPE (mean square predictor error), idealmente pequeño
mean(thiessen.cv.a$residual^2)## [1] 72442.47mean(idw.cv.a$residual^2)## [1] 47302.6mean(ok.cv.a$residual^2)## [1] 0.03747323#Error medio cuadrático (RMSE) es una medida general. Idealmente pequeño
sqrt(sum(thiessen.cv.a$residual^2)/length(thiessen.cv.a$residual))## [1] 269.1514sqrt(sum(idw.cv.a$residual^2)/length(idw.cv.a$residual))## [1] 217.4916sqrt(sum(ok.cv.a$residual^2)/length(ok.cv.a$residual))## [1] 0.19358#Varianza de residuos, idealmente cercanos a cero
var(thiessen.cv.a$residual, na.rm=T)## [1] 72896.19var(idw.cv.a$residual, na.rm=T)## [1] 47784.48var(ok.cv.a$residual, na.rm=T)## [1] 0.03785947# Cantidad de variación explicada por el modelo en cada método, idealmente cercano a 100
(1-var(thiessen.cv.a$residual, na.rm=T)/var(PREC.2008.a.XY$PREC))*100## [1] 26.53841(1-var(idw.cv.a$residual, na.rm=T)/var(PREC.2008.a.XY$PREC))*100## [1] 51.84489(1-var(ok.cv.a$residual, na.rm=T)/var(PREC.2008.a.XY$PREC))*100## [1] 99.99996
Ha llegado al final de la lección!