Maestría en Ecohidrología
Universidad de Cuenca
http://www.ucuenca.edu.ec/maestria-ecohidrologia/

Daniela Ballari (PhD)
daniela.ballari@ucuenca.edu.ec
http://www.researchgate.net/profile/Daniela_Ballari/publications

Curso completo en: http://rpubs.com/daniballari/analisis_espacial_indice



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

Temario

A- Cargar librerias y datos
B- 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- Ejemplo de datos 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 

library(rgdal)
library(gstat)

#definir directorio de trabajo
setwd("C:/R_ecohidrologia/Geoestadistica")
## OGR data source with driver: ESRI Shapefile 
## Source: ".", layer: "countries"
## with 6 features
## It has 6 fields
# Cargar los datos de countries.shp a una variable llamada limite
# Cargar los datos de GSOD_2008_bilogora.csv a una variable llamada PREC.2008

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 

# 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     1
PREC.20080501.XY <- remove.duplicates(PREC.20080501.XY)
hist(PREC.20080501.XY$PREC)

# Se aplicará log para mejorar la normalización de los datos. 
#log1p = log(x+1) se utiliza para evitar tener valores negativos resultado de log
hist(log1p(PREC.20080501.XY$PREC)) 

#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

B.1- Estructura espacial local: autocorrelación

#Variograma empírico
PREC.ve.d <- variogram(log1p(PREC)~1, PREC.20080501.XY) 
plot(PREC.ve.d, pl = T)

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.5
plot(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
## 797  (434744.4, 5264761)      0.00       NA     1.52     1.52     NA    1
## 1161 (481064.7, 5227420)      1.02       NA     2.03     1.01     NA    2
## 1525 (505012.4, 5299631)      0.00       NA     0.00     0.00     NA    3
## 1966   (522465, 5294115)      0.00       NA     0.00     0.00     NA    4
## 2133 (549000.1, 5255378)      5.08       NA     7.11     2.03     NA    5
idw.cv.d[1:5,]
##              coordinates var1.pred var1.var observed   residual zscore
## 797  (434744.4, 5264761)  1.738358       NA     1.52 -0.2183581     NA
## 1161 (481064.7, 5227420)  2.658482       NA     2.03 -0.6284822     NA
## 1525 (505012.4, 5299631)  2.032641       NA     0.00 -2.0326415     NA
## 1966   (522465, 5294115)  2.376153       NA     0.00 -2.3761529     NA
## 2133 (549000.1, 5255378)  4.375667       NA     7.11  2.7343332     NA
##      fold
## 797     1
## 1161    2
## 1525    3
## 1966    4
## 2133    5
ok.cv.d[1:5,]
##              coordinates var1.pred  var1.var  observed   residual
## 797  (434744.4, 5264761) 0.3287548 0.5696455 0.9242589  0.5955041
## 1161 (481064.7, 5227420) 0.7408164 0.5236523 1.1085626  0.3677462
## 1525 (505012.4, 5299631) 0.2126372 0.5191697 0.0000000 -0.2126372
## 1966   (522465, 5294115) 0.4035653 0.4818973 0.0000000 -0.4035653
## 2133 (549000.1, 5255378) 1.6770895 0.4539861 2.0930979  0.4160084
##          zscore fold
## 797   0.7890099    1
## 1161  0.5081907    2
## 1525 -0.2951105    3
## 1966 -0.5813486    4
## 2133  0.6174202    5
ok.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)
## NULL
print(plot(var1.pred~observed,idw.cv.d, main="IDW"), split=c(2,1,3,1), more=TRUE)
## NULL
print(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.4627582
cor(idw.cv.d$var1.pred,idw.cv.d$observed)
## [1] 0.6617524
cor(ok.cv.d$var1.pred,ok.cv.d$obs)
## [1] 0.7272514
#media de residuos, idealmente 0
mean(thiessen.cv.d$residual)
## [1] 0.1426531
mean(idw.cv.d$residual)
## [1] -0.1804113
mean(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.252934
sd(idw.cv.d$residual)
## [1] 2.606365
sd(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.49396
mean(idw.cv.d$residual^2)
## [1] 6.756368
mean(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.239438
sqrt(sum(idw.cv.d$residual^2)/length(idw.cv.d$residual))
## [1] 2.599301
sqrt(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.58158
var(idw.cv.d$residual, na.rm=T)
## [1] 6.793137
var(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. 

# 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.72
head(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

hist(PREC.2008.a.XY$PREC)

hist(log1p(PREC.2008.a.XY$PREC))

C.1- Estructura espacial local: autocorrelación

# Variograma empírico
PREC.ve.a <- variogram(log1p(PREC)~1, PREC.2008.a.XY)
plot(PREC.ve.a,pl = T)

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.2
plot(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
## 3 (434744.4, 5264761)   1158.71       NA   874.31  -284.40     NA    1
## 4 (481064.7, 5227420)    712.76       NA   745.99    33.23     NA    2
## 5 (505012.4, 5299631)    870.72       NA  1772.18   901.46     NA    3
## 6   (522465, 5294115)   1772.18       NA   870.72  -901.46     NA    4
## 7 (549000.1, 5255378)   1005.03       NA  1005.85     0.82     NA    5
idw.cv.a[1:5,]
##           coordinates var1.pred var1.var observed   residual zscore fold
## 3 (434744.4, 5264761)  960.5473       NA   874.31  -86.23730     NA    1
## 4 (481064.7, 5227420)  837.4586       NA   745.99  -91.46858     NA    2
## 5 (505012.4, 5299631)  845.1617       NA  1772.18  927.01833     NA    3
## 6   (522465, 5294115) 1090.8605       NA   870.72 -220.14047     NA    4
## 7 (549000.1, 5255378)  835.3164       NA  1005.85  170.53358     NA    5
ok.cv.a[1:5,]
##           coordinates var1.pred   var1.var observed    residual     zscore
## 3 (434744.4, 5264761)  6.911983 0.05847249 6.774578 -0.13740463 -0.5682318
## 4 (481064.7, 5227420)  6.703703 0.05149363 6.616052 -0.08765163 -0.3862631
## 5 (505012.4, 5299631)  6.753763 0.05390556 7.480530  0.72676726  3.1302467
## 6   (522465, 5294115)  7.018383 0.04935535 6.770468 -0.24791499 -1.1159267
## 7 (549000.1, 5255378)  6.667441 0.04538308 6.914582  0.24714046  1.1601039
##   fold     obs
## 3    1  874.31
## 4    2  745.99
## 5    3 1772.18
## 6    4  870.72
## 7    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)
## NULL
print(plot(var1.pred~observed,idw.cv.a, main="IDW"), split=c(2,1,3,1), more=TRUE)
## NULL
print(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.6756613
cor(idw.cv.a$var1.pred,idw.cv.a$observed)
## [1] 0.7475555
cor(ok.cv.a$var1.pred,ok.cv.a$obs)
## [1] 0.7798622
#media de residuos, idealmente 0
mean(thiessen.cv.a$residual)
## [1] -17.03296
mean(idw.cv.a$residual)
## [1] -2.390472
mean(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.9929
sd(idw.cv.a$residual)
## [1] 218.5966
sd(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.47
mean(idw.cv.a$residual^2)
## [1] 47302.6
mean(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.1514
sqrt(sum(idw.cv.a$residual^2)/length(idw.cv.a$residual))
## [1] 217.4916
sqrt(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.19
var(idw.cv.a$residual, na.rm=T)
## [1] 47784.48
var(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 este material!