Geoestadística: variables secundarias o explicatorias

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. Aplicar geoestadística para predecir la distribución espacial de la precipitación anual, utilizando variables secundarias o explicatorias.

Las variables secundarias a utilizar son imágenes de radar (ground-based) y la topografía a través de un DEM (digital elevation model).

Esta práctica se basa parcialmente 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- Seleccionar ejemplo de datos anuales acumulados
C- Kriging ordinario
D- Kriging Universal con variables secundarias
E- Validación cruzada

A-Cargar librerias y datos

library(rgdal)
library(gstat)
library(RSAGA)
library(raster)

#definir directorio de trabajo
setwd("C:/R_ecohidrologia/Geoestadistica")

#Definir sistema de referencia UTM
utm33 <- "+proj=utm +zone=33 +ellps=WGS84 +datum=WGS84 +units=m +no_defs" 

#Datos
limite <- readOGR(".", "countries")

## OGR data source with driver: ESRI Shapefile 
## Source: ".", layer: "countries"
## with 6 features
## It has 6 fields

  #Convertir paises a capa de lineas
  rsaga.geoprocessor(lib="shapes_lines", module=0, param=list(LINES="borders.shp",POLYGONS="countries.shp"))

## 
## 
## library path:    C:\Progra~1\SAGA-GIS\modules\shapes_lines.dll
## library name:    Shapes - Lines
## module name :    Convert Polygons to Lines
## author      :    O.Conrad (c) 2005
## 
## Load shapes: countries.shp...
## 
## 
## Parameters
## 
## 
## Lines: Lines
## Polygons: countries
## 
## Save shapes: borders.shp...

  borders <- readOGR(".", "borders")#leer archivo de lineas

## OGR data source with driver: ESRI Shapefile 
## Source: ".", layer: "borders"
## with 6 features
## It has 6 fields

PREC.2008 <- read.table("GSOD_2008_bilogora.csv", sep=";", header=TRUE)

# Leer variables auxiliares
grids1km <- readGDAL("dem.asc")#modelo digital de elevaciones

## dem.asc has GDAL driver AAIGrid 
## and has 408 rows and 408 columns

  names(grids1km) <- "dem"
  proj4string(grids1km) <- CRS(utm33)

## Warning in `proj4string<-`(`*tmp*`, value = <S4 object of class structure("CRS", package = "sp")>): A new CRS was assigned to an object with an existing CRS:
## +proj=utm +zone=33 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## without reprojecting.
## For reprojection, use function spTransform in package rgdal

  #crear raster con valores de LAT, LON
  grids.ll <- spTransform(grids1km, CRS("+init=epsg:4326"))

## Warning in spTransform(grids1km, CRS("+init=epsg:4326")): Grid warping not
## available, coercing to points

## Warning in spTransform(as(x, "SpatialPixelsDataFrame"), CRSobj, ...): Grid
## warping not available, coercing to points

  grids1km$LAT <- grids.ll@coords[,2]
  grids1km$LON <- grids.ll@coords[,1]

# Derivar Wetness Index a partir del DEM. http://en.wikipedia.org/wiki/Topographic_Wetness_Index
writeGDAL(grids1km["dem"], "dem.sdat", "SAGA")
  rsaga.wetness.index("dem.sgrd","swi.sgrd")

## 
## 
## library path:    C:\Progra~1\SAGA-GIS\modules\ta_hydrology.dll
## library name:    Terrain Analysis - Hydrology
## module name :    SAGA Wetness Index
## author      :    (c) 2001 by J.Boehner, O.Conrad
## 
## Load grid: dem.sgrd...
## 
## 
## Parameters
## 
## 
## Grid system: 1000; 408x 408y; 467000x 4879000y
## Elevation: dem
## Catchment area: Catchment area
## Catchment slope: Catchment slope
## Modified Catchment Area: Modified Catchment Area
## Topographic Wetness Index: Topographic Wetness Index
## Suction: 10.000000
## Type of Area: square root of catchment area
## Type of Slope: catchment slope
## Minimum Slope: 0.000000
## Offset Slope: 0.100000
## Slope Weighting: 1.000000
## 
## catchment area and slope...
## Create index: dem
## pass 1 (74495 > 0)
## pass 2 (48243 > 0)
## pass 3 (45858 > 0)
## pass 4 (37650 > 0)
## pass 5 (39993 > 0)
## pass 6 (34970 > 0)
## pass 7 (32188 > 0)
## pass 8 (31622 > 0)
## pass 9 (31500 > 0)
## pass 10 (30575 > 0)
## pass 11 (29629 > 0)
## pass 12 (29181 > 0)
## pass 13 (28386 > 0)
## pass 14 (28288 > 0)
## pass 15 (27698 > 0)
## pass 16 (27571 > 0)
## pass 17 (27099 > 0)
## pass 18 (26711 > 0)
## pass 19 (25180 > 0)
## pass 20 (24559 > 0)
## pass 21 (24304 > 0)
## pass 22 (23629 > 0)
## pass 23 (23278 > 0)
## pass 24 (21349 > 0)
## pass 25 (20561 > 0)
## pass 26 (19582 > 0)
## pass 27 (19254 > 0)
## pass 28 (18642 > 0)
## pass 29 (17140 > 0)
## pass 30 (15633 > 0)
## pass 31 (14560 > 0)
## pass 32 (13233 > 0)
## pass 33 (13034 > 0)
## pass 34 (11211 > 0)
## pass 35 (10562 > 0)
## pass 36 (9614 > 0)
## pass 37 (9357 > 0)
## pass 38 (8608 > 0)
## pass 39 (7924 > 0)
## pass 40 (7703 > 0)
## pass 41 (7970 > 0)
## pass 42 (7785 > 0)
## pass 43 (7321 > 0)
## pass 44 (7258 > 0)
## pass 45 (7006 > 0)
## pass 46 (7030 > 0)
## pass 47 (6782 > 0)
## pass 48 (6824 > 0)
## pass 49 (6552 > 0)
## pass 50 (6506 > 0)
## pass 51 (6442 > 0)
## pass 52 (6060 > 0)
## pass 53 (5831 > 0)
## pass 54 (5572 > 0)
## pass 55 (5443 > 0)
## pass 56 (5076 > 0)
## pass 57 (5234 > 0)
## pass 58 (4816 > 0)
## pass 59 (4799 > 0)
## pass 60 (4470 > 0)
## pass 61 (4391 > 0)
## pass 62 (4221 > 0)
## pass 63 (4368 > 0)
## pass 64 (4157 > 0)
## pass 65 (4078 > 0)
## pass 66 (3901 > 0)
## pass 67 (3711 > 0)
## pass 68 (3716 > 0)
## pass 69 (3669 > 0)
## pass 70 (4109 > 0)
## pass 71 (4154 > 0)
## pass 72 (4053 > 0)
## pass 73 (4151 > 0)
## pass 74 (4257 > 0)
## pass 75 (4082 > 0)
## pass 76 (3872 > 0)
## pass 77 (4128 > 0)
## pass 78 (4214 > 0)
## pass 79 (4174 > 0)
## pass 80 (3934 > 0)
## pass 81 (4296 > 0)
## pass 82 (4375 > 0)
## pass 83 (4485 > 0)
## pass 84 (4455 > 0)
## pass 85 (4071 > 0)
## pass 86 (4224 > 0)
## pass 87 (4454 > 0)
## pass 88 (4501 > 0)
## pass 89 (4178 > 0)
## pass 90 (4149 > 0)
## pass 91 (4141 > 0)
## pass 92 (4261 > 0)
## pass 93 (4143 > 0)
## pass 94 (3600 > 0)
## pass 95 (3754 > 0)
## pass 96 (4090 > 0)
## pass 97 (3860 > 0)
## pass 98 (3583 > 0)
## pass 99 (3628 > 0)
## pass 100 (3616 > 0)
## pass 101 (3605 > 0)
## pass 102 (3409 > 0)
## pass 103 (2731 > 0)
## pass 104 (2495 > 0)
## pass 105 (2762 > 0)
## pass 106 (2644 > 0)
## pass 107 (2316 > 0)
## pass 108 (2054 > 0)
## pass 109 (1772 > 0)
## pass 110 (1565 > 0)
## pass 111 (1412 > 0)
## pass 112 (1342 > 0)
## pass 113 (1221 > 0)
## pass 114 (1126 > 0)
## pass 115 (1095 > 0)
## pass 116 (1072 > 0)
## pass 117 (889 > 0)
## pass 118 (826 > 0)
## pass 119 (839 > 0)
## pass 120 (749 > 0)
## pass 121 (655 > 0)
## pass 122 (593 > 0)
## pass 123 (575 > 0)
## pass 124 (485 > 0)
## pass 125 (428 > 0)
## pass 126 (385 > 0)
## pass 127 (342 > 0)
## pass 128 (283 > 0)
## pass 129 (236 > 0)
## pass 130 (209 > 0)
## pass 131 (184 > 0)
## pass 132 (147 > 0)
## pass 133 (116 > 0)
## pass 134 (87 > 0)
## pass 135 (67 > 0)
## pass 136 (55 > 0)
## pass 137 (44 > 0)
## pass 138 (28 > 0)
## pass 139 (17 > 0)
## pass 140 (10 > 0)
## pass 141 (3 > 0)
## pass 142 (0 > 0)
## post-processing...
## topographic wetness index...
## Save grid: C:\Users\dani\AppData\Local\Temp\RtmpiGiIQQ\file18385e377850...
## Save grid: C:\Users\dani\AppData\Local\Temp\RtmpiGiIQQ\file18382dee2b0...
## Save grid: C:\Users\dani\AppData\Local\Temp\RtmpiGiIQQ\file18386fc650a9...
## Save grid: swi.sgrd...

  grids1km$swi <- readGDAL("swi.sdat")$band1

## swi.sdat has GDAL driver SAGA 
## and has 408 rows and 408 columns

names(grids1km)

## [1] "dem" "LAT" "LON" "swi"

#Leer imagenes de radar
#crear un listado de imagenes disponibles. Deben estar dentro de un directorio llamado RN
listado <- dir(path=paste(getwd(), "/RN", sep=""), 
               pattern=glob2rx("RN1_2008_*.tif"), recursive=FALSE, full.names=TRUE) 
listado

##  [1] "C:/R_ecohidrologia/Geoestadistica/RN/RN1_2008_01.tif"
##  [2] "C:/R_ecohidrologia/Geoestadistica/RN/RN1_2008_02.tif"
##  [3] "C:/R_ecohidrologia/Geoestadistica/RN/RN1_2008_03.tif"
##  [4] "C:/R_ecohidrologia/Geoestadistica/RN/RN1_2008_04.tif"
##  [5] "C:/R_ecohidrologia/Geoestadistica/RN/RN1_2008_05.tif"
##  [6] "C:/R_ecohidrologia/Geoestadistica/RN/RN1_2008_06.tif"
##  [7] "C:/R_ecohidrologia/Geoestadistica/RN/RN1_2008_07.tif"
##  [8] "C:/R_ecohidrologia/Geoestadistica/RN/RN1_2008_08.tif"
##  [9] "C:/R_ecohidrologia/Geoestadistica/RN/RN1_2008_09.tif"
## [10] "C:/R_ecohidrologia/Geoestadistica/RN/RN1_2008_10.tif"
## [11] "C:/R_ecohidrologia/Geoestadistica/RN/RN1_2008_11.tif"
## [12] "C:/R_ecohidrologia/Geoestadistica/RN/RN1_2008_12.tif"

#leer primera imagen
radar.mensual <- readGDAL(listado[1])

## C:/R_ecohidrologia/Geoestadistica/RN/RN1_2008_01.tif has GDAL driver GTiff 
## and has 408 rows and 408 columns

#Asignar como nombre de la imagen el mes al que corresponde. 
#La función "Format" y "%b" extrae el mes (abreviado). 
#Otras opciones son %Y para extraer el año y %d para el dia.
names(radar.mensual) <- format(as.Date(paste("2008-", 1, "-01", sep="")), "%b")
names(radar.mensual) 

## [1] "Jan"

#Leer desde la imágen 2 hasta la 12 y acumular en el SpatialGridDataFrame "radar.mensual".
#El script lee la imagen en la posición j y asigna el nombre del mes respectivo.
#Silent=TRUE es para que no muestre resultados intermedios en la consola. 
for(j in 2:length(listado)) 
  { 
radar.mensual@data[,format(as.Date(paste("2008-", j, "-01", sep="")), "%b")] <- 
  readGDAL(listado[j], silent=TRUE)$band1 
  }
names(radar.mensual)

##  [1] "Jan" "Feb" "Mar" "Apr" "May" "Jun" "Jul" "Aug" "Sep" "Oct" "Nov"
## [12] "Dec"

#visualizar las 12 imagenes. Esto tardará un momento ya que son 12 imágenes.
windows(width=14, height=6)
spplot(radar.mensual, at=seq(0,9,9/40)^3, col.regions=grey(rev(seq(0,1,0.025))),
       sp.layout=list("sp.lines", col="black", borders))

dev.off() 

## png 
##   2

# Calcular la precipitación acumulada anual medida con el radar
grids1km$radar.anual <- rowSums(radar.mensual@data[,1:12], na.rm=FALSE, dims=1) 

names(grids1km)

## [1] "dem"         "LAT"         "LON"         "swi"         "radar.anual"

#Visualizar
radarplot <- spplot(grids1km["radar.anual"], at=seq(5,15,(15-5)/40)^3, 
                    col.regions=grey(rev(seq(0,1,0.025))), 
                    sp.layout=list("sp.lines", col="black", borders), 
                    main="Precipitación radar acumulada anual")
radarplot

B-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)   
# Seleccionar coordenadas de estaciones para mantener georeferenciación.
PREC.2008.a.lon.lat <- unique(PREC.2008[!is.na(PREC.2008$PREC),c("STNID", "LON", "LAT")])

#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)
#Asignar coordenadas, sistema de referencia y transformar
coordinates(PREC.2008.a.unido) <- ~LON+LAT
proj4string(PREC.2008.a.unido) <- CRS("+proj=longlat +datum=WGS84")
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 y seleccionar valores no vacios
PREC.2008.a.XY$PREC2 <- ifelse(PREC.2008.a.XY$PREC==0, NA, PREC.2008.a.XY$PREC)
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- Kriging Ordinario

# Variograma empírico
PREC.ve.a <- variogram(log1p(PREC)~1, PREC.2008.a.XY)
# 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)))
#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)

# 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

#Inverso de la distancia (idw)
idw.a = idw(PREC ~ 1, PREC.2008.a.XY, PREC.2008.a.XY_grid)

## [inverse distance weighted interpolation]

# 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

pts.s <- list("sp.points", PREC.2008.a.XY, col="white",pch=19)
spplot(idw.a, "var1.pred", at= seq(400,2100,50), zlim=c(400,2100), asp=1, 
             col.regions=gray(seq(1,0.1,l=35)),sp.layout = list(pts.s),main="IDW")

spplot(ok.a, "PREC.pred", at= seq(400,2100,50), zlim=c(400,2100), asp=1, 
             col.regions=gray(seq(1,0.1,l=35)),main="OK anual PREC mm",
             sp.layout = list(pts.s))

spplot(ok.a, "var1.var", asp=1, col.regions=gray(seq(1,0.1,l=30)),
             main="Varianza OK, PREC mm",sp.layout = list(pts.s))

graphics.off()

D- Kriging Universal con variables secundarias

#Extraer datos de radar para cada estación superponiendo (extract) ambas variables
PREC.2008.ov <- extract(brick(grids1km), PREC.2008.a.XY) 
PREC.2008.ov <- data.frame(PREC.2008.a.XY@data, data.frame(PREC.2008.ov))#unir a datos existentes

#Se observan valores NA, ello es porque no existen datos de grilla en estas zonas. 
head(PREC.2008.ov)

##       estacion    PREC   PREC2 dem      LAT      LON       swi radar.anual
## 3 111570-99999  874.31  874.31  NA       NA       NA        NA          NA
## 4 111650-99999  745.99  745.99 661 47.19622 14.74916 10.207204      1998.4
## 5 111700-99999 1772.18 1772.18  NA       NA       NA        NA          NA
## 6 111710-99999  870.72  870.72  NA       NA       NA        NA          NA
## 7 111730-99999 1005.85 1005.85 976 47.44660 15.64996  9.054108      2344.7
## 8 111740-99999  712.76  712.76 740 47.33147 15.00000  8.729557       365.0

PREC.2008.ov <- subset(PREC.2008.ov, !is.na(PREC.2008.ov$radar.anual)) 

Exploración de relaciones entre variables

## [1] 0.5118432

## [1] 0.1177819

## [1] -0.3599588

## [1] -0.3590517

## [1] -0.6055487

cor(PREC.2008.ov[,c("PREC", "dem" ,"LAT" , "LON","swi" , "radar.anual")])

##                   PREC        dem        LAT        LON        swi
## PREC         1.0000000  0.5118432 -0.3590517 -0.6055487 -0.3599588
## dem          0.5118432  1.0000000  0.2674194 -0.4718539 -0.5912024
## LAT         -0.3590517  0.2674194  1.0000000 -0.1148591 -0.1674120
## LON         -0.6055487 -0.4718539 -0.1148591  1.0000000  0.6306198
## swi         -0.3599588 -0.5912024 -0.1674120  0.6306198  1.0000000
## radar.anual  0.1177819  0.5222903  0.3584995 -0.2513354 -0.4189161
##             radar.anual
## PREC          0.1177819
## dem           0.5222903
## LAT           0.3584995
## LON          -0.2513354
## swi          -0.4189161
## radar.anual   1.0000000

plot(PREC.2008.ov[,c("PREC", "dem" ,"LAT" , "LON","swi" , "radar.anual")])

#Regresión lineal (https://www0.gsb.columbia.edu/premba/analytical/s7/s7_6.cfm)
lm1 <- lm(log1p(PREC) ~ LAT + LON, PREC.2008.ov)
summary(lm1)

## 
## Call:
## lm(formula = log1p(PREC) ~ LAT + LON, data = PREC.2008.ov)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.44456 -0.12014  0.00621  0.09828  0.49563 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 15.42396    1.29282  11.930  < 2e-16 ***
## LAT         -0.13849    0.02663  -5.201 2.22e-06 ***
## LON         -0.14511    0.01660  -8.741 1.60e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1965 on 64 degrees of freedom
## Multiple R-squared:  0.5956, Adjusted R-squared:  0.583 
## F-statistic: 47.13 on 2 and 64 DF,  p-value: 2.62e-13

lm2 <- lm(log1p(PREC) ~ dem+radar.anual+swi, PREC.2008.ov)
summary(lm2)

## 
## Call:
## lm(formula = log1p(PREC) ~ dem + radar.anual + swi, data = PREC.2008.ov)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.56340 -0.19686  0.02745  0.12359  0.66687 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.938e+00  3.231e-01  21.471  < 2e-16 ***
## dem          4.223e-04  1.455e-04   2.903  0.00509 ** 
## radar.anual -6.378e-05  5.647e-05  -1.129  0.26301    
## swi         -3.847e-02  2.619e-02  -1.469  0.14691    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2693 on 63 degrees of freedom
## Multiple R-squared:  0.2521, Adjusted R-squared:  0.2165 
## F-statistic: 7.078 on 3 and 63 DF,  p-value: 0.0003569

lm3 <- lm(log1p(PREC) ~ LAT+LON+dem+radar.anual+swi, PREC.2008.ov)
summary(lm3)

## 
## Call:
## lm(formula = log1p(PREC) ~ LAT + LON + dem + radar.anual + swi, 
##     data = PREC.2008.ov)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.40777 -0.08999 -0.00010  0.11859  0.37845 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.596e+01  1.200e+00  13.298  < 2e-16 ***
## LAT         -1.662e-01  2.547e-02  -6.525 1.51e-08 ***
## LON         -1.289e-01  1.931e-02  -6.675 8.36e-09 ***
## dem          3.723e-04  9.693e-05   3.841 0.000294 ***
## radar.anual  2.395e-05  3.836e-05   0.624 0.534754    
## swi          3.096e-02  1.973e-02   1.569 0.121868    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1757 on 61 degrees of freedom
## Multiple R-squared:  0.6919, Adjusted R-squared:  0.6667 
## F-statistic:  27.4 on 5 and 61 DF,  p-value: 2.008e-14

Variograma y Kriging universal

coordinates(PREC.2008.ov) <- ~LON+LAT
proj4string(PREC.2008.ov) <- CRS("+proj=longlat +datum=WGS84")
PREC.2008.ov <- spTransform(PREC.2008.ov, CRS(utm33))
#Variograma empírico, nube de variograma 
PREC.ve.a.u <- variogram(log1p(PREC) ~ LAT+LON+dem+radar.anual+swi, PREC.2008.ov)
# variograma inicial
PREC.vi.a.u <- vgm(nugget=0, model="Exp", 
                   range=sqrt(diff(PREC.2008.ov@bbox["LON",])^2 + 
                   diff(PREC.2008.ov@bbox["LAT",])^2)/4, 
                   psill=var(log1p(PREC.2008.ov$PREC)))
#Ajuste del variograma (teórico)
PREC.vt.a.u <- fit.variogram(PREC.ve.a.u, model=PREC.vi.a.u)

#Kriging universal
PREC.uk.a <- krige(log1p(PREC) ~ LAT+LON+dem+radar.anual+swi, PREC.2008.ov, newdata = grids1km, model = PREC.vt.a.u) 

## [using universal kriging]

#new data: debe ser el raster de predicción que contenga las variables secundarias a utilizar.

PREC.uk.a$PREC.pred <- expm1(PREC.uk.a$var1.pred) #volver a valores de precipitación

par(mfrow=c(2,1))
pts.s <- list("sp.points", PREC.2008.ov, col="black",pch=1)
print(spplot(PREC.uk.a, "PREC.pred", asp=1,col.regions=gray(seq(1,0.1,l=35)),
              main="Predicción UK, PREC mm",sp.layout = list(pts.s)), split=c(1,1,2,1), more=TRUE)
print(spplot(PREC.uk.a, "var1.var", asp=1, col.regions=gray(seq(1,0.1,l=30)),
             main="Varianza UK, PREC mm",sp.layout = list(pts.s)), 
             split=c(2,1,2,1), more=FALSE)

E- Validación cruzada

idw.cv <- krige.cv(PREC ~ 1, PREC.2008.a.XY)
ok.cv <- krige.cv(log1p(PREC) ~ 1, locations = PREC.2008.a.XY, model = PREC.vt.a)
uk.cv <- krige.cv(log1p(PREC) ~ LAT+LON+dem+radar.anual+swi, locations = PREC.2008.ov,model = PREC.vt.a.u)

idw.cv$obs<- PREC.2008.a.XY$PREC
ok.cv$obs <- PREC.2008.a.XY$PREC
uk.cv$obs <-PREC.2008.ov$PREC

summary(idw.cv)

## Object of class SpatialPointsDataFrame
## Coordinates:
##           min       max
## LON  434744.4  892122.9
## LAT 4876238.0 5315499.4
## Is projected: NA 
## proj4string : [NA]
## Number of points: 98
## Data attributes:
##    var1.pred         var1.var      observed         residual      
##  Min.   : 569.6   Min.   : NA   Min.   : 432.1   Min.   :-378.61  
##  1st Qu.: 679.3   1st Qu.: NA   1st Qu.: 591.6   1st Qu.:-111.34  
##  Median : 792.2   Median : NA   Median : 746.8   Median : -53.40  
##  Mean   : 820.9   Mean   :NaN   Mean   : 818.5   Mean   :  -2.39  
##  3rd Qu.: 906.8   3rd Qu.: NA   3rd Qu.: 949.7   3rd Qu.:  49.70  
##  Max.   :1348.9   Max.   : NA   Max.   :2056.1   Max.   : 928.18  
##                   NA's   :98                                      
##      zscore         fold            obs        
##  Min.   : NA   Min.   : 1.00   Min.   : 432.1  
##  1st Qu.: NA   1st Qu.:25.25   1st Qu.: 591.6  
##  Median : NA   Median :49.50   Median : 746.8  
##  Mean   :NaN   Mean   :49.50   Mean   : 818.5  
##  3rd Qu.: NA   3rd Qu.:73.75   3rd Qu.: 949.7  
##  Max.   : NA   Max.   :98.00   Max.   :2056.1  
##  NA's   :98

summary(ok.cv)

## Object of class SpatialPointsDataFrame
## Coordinates:
##           min       max
## LON  434744.4  892122.9
## LAT 4876238.0 5315499.4
## Is projected: NA 
## proj4string : [NA]
## Number of points: 98
## Data attributes:
##    var1.pred        var1.var          observed        residual         
##  Min.   :6.268   Min.   :0.03671   Min.   :6.071   Min.   :-0.5375807  
##  1st Qu.:6.411   1st Qu.:0.04744   1st Qu.:6.384   1st Qu.:-0.0951852  
##  Median :6.626   Median :0.05226   Median :6.617   Median : 0.0006323  
##  Mean   :6.648   Mean   :0.05342   Mean   :6.648   Mean   : 0.0002768  
##  3rd Qu.:6.828   3rd Qu.:0.05737   3rd Qu.:6.857   3rd Qu.: 0.1076166  
##  Max.   :7.254   Max.   :0.08367   Max.   :7.629   Max.   : 0.7267673  
##      zscore                fold            obs        
##  Min.   :-2.5092450   Min.   : 1.00   Min.   : 432.1  
##  1st Qu.:-0.3950501   1st Qu.:25.25   1st Qu.: 591.6  
##  Median : 0.0028543   Median :49.50   Median : 746.8  
##  Mean   : 0.0007853   Mean   :49.50   Mean   : 818.5  
##  3rd Qu.: 0.4420903   3rd Qu.:73.75   3rd Qu.: 949.7  
##  Max.   : 3.1302467   Max.   :98.00   Max.   :2056.1

summary(uk.cv)

## Object of class SpatialPointsDataFrame
## Coordinates:
##         min     max
## LON  481000  868000
## LAT 4900000 5283000
## Is projected: NA 
## proj4string : [NA]
## Number of points: 67
## Data attributes:
##    var1.pred        var1.var           observed        residual         
##  Min.   :6.056   Min.   :0.008305   Min.   :6.071   Min.   :-0.5874215  
##  1st Qu.:6.399   1st Qu.:0.023091   1st Qu.:6.379   1st Qu.:-0.0623772  
##  Median :6.608   Median :0.024363   Median :6.582   Median : 0.0162590  
##  Mean   :6.596   Mean   :0.024541   Mean   :6.595   Mean   :-0.0009843  
##  3rd Qu.:6.778   3rd Qu.:0.025920   3rd Qu.:6.781   3rd Qu.: 0.1098597  
##  Max.   :7.623   Max.   :0.063625   Max.   :7.538   Max.   : 0.3568743  
##      zscore               fold           obs        
##  Min.   :-3.165656   Min.   : 1.0   Min.   : 432.1  
##  1st Qu.:-0.410090   1st Qu.:17.5   1st Qu.: 588.6  
##  Median : 0.097654   Median :34.0   Median : 721.1  
##  Mean   : 0.004916   Mean   :34.0   Mean   : 766.7  
##  3rd Qu.: 0.710471   3rd Qu.:50.5   3rd Qu.: 879.5  
##  Max.   : 2.236873   Max.   :67.0   Max.   :1876.5

#Correlación
par(mfrow=c(1,3))
print(plot(var1.pred~obs,idw.cv, main="idw"), split=c(1,1,3,1), more=TRUE)

## NULL

print(plot(var1.pred~obs,ok.cv, main="ok"), split=c(2,1,3,1), more=TRUE)

## NULL

print(plot(var1.pred~obs,uk.cv, main="uK"), split=c(3,1,2,1), more=FALSE)

## NULL

#correlación, idealmente 1
cor(idw.cv$var1.pred,idw.cv$obs)

## [1] 0.7475555

cor(ok.cv$var1.pred,ok.cv$obs)

## [1] 0.7798622

cor(uk.cv$var1.pred,uk.cv$obs)

## [1] 0.8207795

#media de residuos, idealmente 0
mean(idw.cv$residual)

## [1] -2.390472

mean(ok.cv$residual)

## [1] 0.0002767632

mean(uk.cv$residual)

## [1] -0.0009843329

#Desviaciones estandar del error de interpolación (residuos), idealmente pequeño
sd(idw.cv$residual)

## [1] 218.5966

sd(ok.cv$residual)

## [1] 0.1945751

sd(uk.cv$residual)

## [1] 0.1797338

# MSPE (mean square predictor error), idealmente pequeño
mean(idw.cv$residual^2)

## [1] 47302.6

mean(ok.cv$residual^2)

## [1] 0.03747323

mean(uk.cv$residual^2)

## [1] 0.03182306

#Error medio cuadrático (RMSE) es una medida general. Idealmente pequeño
sqrt(sum(idw.cv$residual^2)/length(idw.cv$residual))

## [1] 217.4916

sqrt(sum(ok.cv$residual^2)/length(ok.cv$residual))

## [1] 0.19358

sqrt(sum(uk.cv$residual^2)/length(uk.cv$residual))

## [1] 0.1783902

#Varianza de residuos, idealmente cercanos a cero
var(idw.cv$residual, na.rm=T)

## [1] 47784.48

var(ok.cv$residual, na.rm=T)

## [1] 0.03785947

var(uk.cv$residual, na.rm=T)

## [1] 0.03230424

#Cantidad de variación explicada por el modelo en cada método, idealmente cercano a 100
(1-var(idw.cv$residual, na.rm=T)/var(PREC.2008.a.XY$PREC))*100

## [1] 51.84489

(1-var(ok.cv$residual, na.rm=T)/var(PREC.2008.a.XY$PREC))*100

## [1] 99.99996

(1-var(uk.cv$residual, na.rm=T)/var(PREC.2008.a.XY$PREC))*100

## [1] 99.99997

Has llegado al final del material!