Geoestadística 3: variables secundarias o explicatorias
Daniela Ballari
Universidad de Cuenca
daniela.ballari@ucuenca.edu.ec
Abril, 2015
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
Pasos a seguir
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
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)
library(RSAGA)## Warning: package 'RSAGA' was built under R version 3.0.3## Loading required package: shapefiles
## Loading required package: foreign
##
## Attaching package: 'shapefiles'
##
## The following objects are masked from 'package:foreign':
##
## read.dbf, write.dbf
##
## Loading required package: plyrlibrary(raster)
#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=0B0ea6usixQHFakhhamtDR1lJeWs")
shell.exec("https://drive.google.com/uc?export=download&id=0B0ea6usixQHFaExTZExYREhuaUU")
#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 and 6 fields
## Feature type: wkbPolygon with 2 dimensions #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 and 6 fields
## Feature type: wkbLineString with 2 dimensionsPREC.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("+proj=longlat +datum=WGS84"))## Warning in spTransform(grids1km, CRS("+proj=longlat +datum=WGS84")): 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 (79333 > 0)
## pass 2 (64535 > 0)
## pass 3 (55417 > 0)
## pass 4 (44157 > 0)
## pass 5 (41795 > 0)
## pass 6 (40685 > 0)
## pass 7 (40516 > 0)
## pass 8 (36919 > 0)
## pass 9 (38435 > 0)
## pass 10 (35282 > 0)
## pass 11 (34408 > 0)
## pass 12 (35048 > 0)
## pass 13 (33480 > 0)
## pass 14 (31538 > 0)
## pass 15 (32266 > 0)
## pass 16 (29792 > 0)
## pass 17 (31493 > 0)
## pass 18 (30593 > 0)
## pass 19 (30305 > 0)
## pass 20 (26255 > 0)
## pass 21 (26954 > 0)
## pass 22 (24546 > 0)
## pass 23 (25476 > 0)
## pass 24 (24393 > 0)
## pass 25 (22093 > 0)
## pass 26 (20155 > 0)
## pass 27 (20547 > 0)
## pass 28 (19893 > 0)
## pass 29 (17634 > 0)
## pass 30 (16787 > 0)
## pass 31 (15191 > 0)
## pass 32 (14389 > 0)
## pass 33 (12682 > 0)
## pass 34 (11539 > 0)
## pass 35 (10347 > 0)
## pass 36 (9668 > 0)
## pass 37 (9398 > 0)
## pass 38 (8873 > 0)
## pass 39 (8505 > 0)
## pass 40 (8041 > 0)
## pass 41 (8046 > 0)
## pass 42 (7899 > 0)
## pass 43 (7447 > 0)
## pass 44 (7398 > 0)
## pass 45 (7522 > 0)
## pass 46 (7409 > 0)
## pass 47 (7103 > 0)
## pass 48 (7092 > 0)
## pass 49 (6934 > 0)
## pass 50 (6582 > 0)
## pass 51 (6279 > 0)
## pass 52 (5964 > 0)
## pass 53 (5849 > 0)
## pass 54 (5781 > 0)
## pass 55 (5544 > 0)
## pass 56 (5108 > 0)
## pass 57 (5032 > 0)
## pass 58 (5004 > 0)
## pass 59 (4801 > 0)
## pass 60 (4589 > 0)
## pass 61 (4518 > 0)
## pass 62 (4565 > 0)
## pass 63 (4309 > 0)
## pass 64 (4161 > 0)
## pass 65 (4107 > 0)
## pass 66 (4081 > 0)
## pass 67 (3881 > 0)
## pass 68 (4171 > 0)
## pass 69 (4377 > 0)
## pass 70 (4366 > 0)
## pass 71 (4018 > 0)
## pass 72 (4250 > 0)
## pass 73 (4450 > 0)
## pass 74 (4478 > 0)
## pass 75 (4490 > 0)
## pass 76 (3788 > 0)
## pass 77 (4386 > 0)
## pass 78 (4488 > 0)
## pass 79 (3808 > 0)
## pass 80 (4330 > 0)
## pass 81 (4533 > 0)
## pass 82 (3853 > 0)
## pass 83 (4505 > 0)
## pass 84 (3934 > 0)
## pass 85 (4404 > 0)
## pass 86 (3800 > 0)
## pass 87 (4443 > 0)
## pass 88 (3729 > 0)
## pass 89 (4291 > 0)
## pass 90 (3656 > 0)
## pass 91 (4268 > 0)
## pass 92 (4270 > 0)
## pass 93 (4210 > 0)
## pass 94 (3449 > 0)
## pass 95 (3978 > 0)
## pass 96 (3963 > 0)
## pass 97 (3090 > 0)
## pass 98 (3563 > 0)
## pass 99 (3442 > 0)
## pass 100 (2820 > 0)
## pass 101 (3211 > 0)
## pass 102 (3136 > 0)
## pass 103 (3320 > 0)
## pass 104 (2554 > 0)
## pass 105 (2785 > 0)
## pass 106 (2274 > 0)
## pass 107 (1971 > 0)
## pass 108 (1737 > 0)
## pass 109 (1517 > 0)
## pass 110 (1409 > 0)
## pass 111 (1374 > 0)
## pass 112 (1252 > 0)
## pass 113 (1116 > 0)
## pass 114 (1170 > 0)
## pass 115 (1078 > 0)
## pass 116 (906 > 0)
## pass 117 (898 > 0)
## pass 118 (827 > 0)
## pass 119 (757 > 0)
## pass 120 (658 > 0)
## pass 121 (631 > 0)
## pass 122 (527 > 0)
## pass 123 (468 > 0)
## pass 124 (407 > 0)
## pass 125 (352 > 0)
## pass 126 (312 > 0)
## pass 127 (261 > 0)
## pass 128 (220 > 0)
## pass 129 (178 > 0)
## pass 130 (163 > 0)
## pass 131 (130 > 0)
## pass 132 (108 > 0)
## pass 133 (91 > 0)
## pass 134 (70 > 0)
## pass 135 (54 > 0)
## pass 136 (44 > 0)
## pass 137 (24 > 0)
## pass 138 (17 > 0)
## pass 139 (7 > 0)
## pass 140 (2 > 0)
## pass 141 (0 > 0)
## post-processing...
## topographic wetness index...
## Save grid: C:\Users\DOCENT~1\AppData\Local\Temp\RtmpETuqxz\filef245cc37c0a...
## Save grid: C:\Users\DOCENT~1\AppData\Local\Temp\RtmpETuqxz\filef2478b5571c...
## Save grid: C:\Users\DOCENT~1\AppData\Local\Temp\RtmpETuqxz\filef24106f6f6d...
## Save grid: swi.sgrd... grids1km$swi <- readGDAL("swi.sdat")$band1## swi.sdat has GDAL driver SAGA
## and has 408 rows and 408 columns#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:/curso_R/RN/RN1_2008_01.tif" "C:/curso_R/RN/RN1_2008_02.tif"
## [3] "C:/curso_R/RN/RN1_2008_03.tif" "C:/curso_R/RN/RN1_2008_04.tif"
## [5] "C:/curso_R/RN/RN1_2008_05.tif" "C:/curso_R/RN/RN1_2008_06.tif"
## [7] "C:/curso_R/RN/RN1_2008_07.tif" "C:/curso_R/RN/RN1_2008_08.tif"
## [9] "C:/curso_R/RN/RN1_2008_09.tif" "C:/curso_R/RN/RN1_2008_10.tif"
## [11] "C:/curso_R/RN/RN1_2008_11.tif" "C:/curso_R/RN/RN1_2008_12.tif"#leer primera imagen
radar.mensual <- readGDAL(listado[1])## C:/curso_R/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] "ene"#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] "ene" "feb" "mar" "abr" "may" "jun" "jul" "ago" "sep" "oct" "nov"
## [12] "dic"#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.2plot(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 (overlay) ambas variables
PREC.2008.ov <- overlay(grids1km, PREC.2008.a.XY) ## Warning: '.local' is deprecated.
## Use 'over' instead.
## See help("Deprecated") and help("sp-deprecated").#aparecerá un mensaje de advertencia ya que es una función en antigua
PREC.2008.ov@data <- cbind(PREC.2008.a.XY@data, PREC.2008.ov@data)#unir a datos existentes
#Se observan valores NA, ello es porque no existen datos de grilla en estas zonas.
#Esto puede observarse en la siguiente imagen
head(PREC.2008.ov@data)## 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.0PREC.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.3767924
## [1] -0.6015258
cor(PREC.2008.ov@data[,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.0000000plot(PREC.2008.ov@data[,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-13lm2 <- 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.0003569lm3 <- 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-14Variograma y Kriging universal
#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, at= seq(400,2100,50), zlim=c(400,2100),
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 :98summary(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.1summary(uk.cv)## Object of class SpatialPointsDataFrame
## Coordinates:
## min max
## LON 481064.7 868438.2
## LAT 4900065.3 5282731.9
## Is projected: NA
## proj4string : [NA]
## Number of points: 67
## Data attributes:
## var1.pred var1.var observed residual
## Min. :6.056 Min. :0.007874 Min. :6.071 Min. :-0.584657
## 1st Qu.:6.400 1st Qu.:0.022906 1st Qu.:6.379 1st Qu.:-0.059341
## Median :6.600 Median :0.024260 Median :6.582 Median : 0.012103
## Mean :6.596 Mean :0.024471 Mean :6.595 Mean :-0.000951
## 3rd Qu.:6.778 3rd Qu.:0.026047 3rd Qu.:6.781 3rd Qu.: 0.111181
## Max. :7.628 Max. :0.063066 Max. :7.538 Max. : 0.356432
## zscore fold obs
## Min. :-3.197934 Min. : 1.0 Min. : 432.1
## 1st Qu.:-0.412339 1st Qu.:17.5 1st Qu.: 588.6
## Median : 0.072216 Median :34.0 Median : 721.1
## Mean : 0.005082 Mean :34.0 Mean : 766.7
## 3rd Qu.: 0.705131 3rd Qu.:50.5 3rd Qu.: 879.5
## Max. : 2.203468 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)## NULLprint(plot(var1.pred~obs,ok.cv, main="ok"), split=c(2,1,3,1), more=TRUE)## NULLprint(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.7475555cor(ok.cv$var1.pred,ok.cv$obs)## [1] 0.7798622cor(uk.cv$var1.pred,uk.cv$obs)## [1] 0.8227579#media de residuos, idealmente 0
mean(idw.cv$residual)## [1] -2.390472mean(ok.cv$residual)## [1] 0.0002767632mean(uk.cv$residual)## [1] -0.0009510359#Desviaciones estandar del error de interpolación (residuos), idealmente pequeño
sd(idw.cv$residual)## [1] 218.5966sd(ok.cv$residual)## [1] 0.1945751sd(uk.cv$residual)## [1] 0.1790096# MSPE (mean square predictor error), idealmente pequeño
mean(idw.cv$residual^2)## [1] 47302.6mean(ok.cv$residual^2)## [1] 0.03747323mean(uk.cv$residual^2)## [1] 0.03156708#Error medio cuadrático (RMSE) es una medida general. Idealmente pequeño
sqrt(sum(idw.cv$residual^2)/length(idw.cv$residual))## [1] 217.4916sqrt(sum(ok.cv$residual^2)/length(ok.cv$residual))## [1] 0.19358sqrt(sum(uk.cv$residual^2)/length(uk.cv$residual))## [1] 0.1776713#Varianza de residuos, idealmente cercanos a cero
var(idw.cv$residual, na.rm=T)## [1] 47784.48var(ok.cv$residual, na.rm=T)## [1] 0.03785947var(uk.cv$residual, na.rm=T)## [1] 0.03204445#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.99997Ha llegado al final de la lección!