5.1 Prediccion

library(gstat)
library(stars)
idw = krige(formula=prcptcn ~ 1, locations= prec, newdata=r.sf)

## [inverse distance weighted interpolation]

idw.sp<- as(idw, "Spatial")
class(idw.sp)

## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"

gridded(idw.sp)<- TRUE
idw.r<- brick(idw.sp)

5.2 Visualizacion

library(rasterVis)

levelplot(idw.r[[1]], margin=FALSE, main="idw")+ latticeExtra::layer(sp.polygons(as(cuenca, "Spatial"), col='black'))

idw.r.masked<- mask(idw.r,as(cuenca, "Spatial"))

levelplot(idw.r.masked[[1]], margin=FALSE, main="idw")+ latticeExtra::layer(sp.polygons(as(cuenca, "Spatial"), col='black'))

library(RColorBrewer)
levelPal <- colorRampPalette(brewer.pal(n = 9, 'Blues'))
levelplot(idw.r.masked[[1]], margin=FALSE, main="idw",col.regions = levelPal)+ latticeExtra::layer(sp.polygons(as(cuenca, "Spatial"), col='black'))+ latticeExtra::layer(sp.points(as(prec, "Spatial"), pch=20, col='black'))

5.3 IDW validacion cruzada

idw.cv<- krige.cv(formula=prcptcn ~ 1, locations= prec)

head(idw.cv)

## Simple feature collection with 6 features and 6 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 705090 ymin: 9662396 xmax: 797962 ymax: 9727551
## epsg (SRID):    NA
## proj4string:    NA
##   var1.pred var1.var  observed  residual zscore fold               geometry
## 1  89.89335       NA 102.80000  12.90664     NA    1 POINT (717092 9662396)
## 2 127.70690       NA 251.20000 123.49310     NA    2 POINT (745529 9662895)
## 3 241.61738       NA 461.20001 219.58262     NA    3 POINT (797962 9699430)
## 4  86.22249       NA  73.50001 -12.72247     NA    4 POINT (705090 9665080)
## 5 211.78507       NA 157.49997 -54.28510     NA    5 POINT (748262 9713142)
## 6 227.24650       NA 340.70025 113.45376     NA    6 POINT (765780 9727551)

library(Metrics)
(idw.mse<- mse(idw.cv$var1.pred, idw.cv$observed))

## [1] 5852.184

(idw.rmse<- rmse(idw.cv$var1.pred, idw.cv$observed))

## [1] 76.49957

(idw.cor<- cor(idw.cv$var1.pred, idw.cv$observed))

## [1] 0.8251744

(idw.bias<- bias(idw.cv$var1.pred, idw.cv$observed))

## [1] -16.51191

plot(idw.cv$var1.pred, idw.cv$observed)

library(ggplot2)
ggplot(idw.cv, aes(observed, y = var1.pred, color = "red")) +
  geom_abline(colour = "#CCCCCC", linetype = 2) +
  geom_point(aes(y = var1.pred, col = "IDW")) +
  geom_smooth(aes(y = var1.pred, col = "IDW"), method = lm, se = FALSE) +
  xlab("observaciones")+
  ylab("predicciones")

6.0.1 Exploracion de atributos - Univariado- Analisis EDA (Exploratory Data Analysis)

prec$prcptcn

##  [1] 102.80000 251.20000 461.20001  73.50001 157.49997 340.70025  36.09999
##  [8] 264.20020 198.30017  75.19997 120.50000  91.20000  26.50000  32.70000
## [15] 271.80000 359.00000 246.80000 314.90000  50.40000  33.10000

sort(prec$prcptcn)

##  [1]  26.50000  32.70000  33.10000  36.09999  50.40000  73.50001  75.19997
##  [8]  91.20000 102.80000 120.50000 157.49997 198.30017 246.80000 251.20000
## [15] 264.20020 271.80000 314.90000 340.70025 359.00000 461.20001

hist(prec$prcptcn, breaks = 16) #Distribuci昼㸳n no sim攼㸹trica, sesgada hacia la derecha

summary(prec$prcptcn) # media mayor a la mediana (sesgada hacia la derecha, distribuci昼㸳n no-normal)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   26.50   67.72  139.00  175.38  266.10  461.20

#Al ser una distribuci昼㸳n no sim攼㸹trica, se aplica logaritmo para transformar
#los valores y obtener una distribuci昼㸳n simetrica (normal). 
#Esto, adem攼㸱s, reduce los posibles outliers.
prec$log.prcptcn <- log10(prec$prcptcn)
hist(prec$log.prcptcn, breaks = 16)

summary(prec$log.prcptcn)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.423   1.825   2.139   2.095   2.425   2.664

head(prec$log.prcptcn)

## [1] 2.011993 2.400020 2.663889 1.866287 2.197280 2.532372

head(prec$log.prcptcn)

## [1] 2.011993 2.400020 2.663889 1.866287 2.197280 2.532372

#Luego para transformar a los valores originales se utiliza antilog. antilog (x)=10^x
head(10^(prec$log.prcptcn))

## [1] 102.80000 251.20000 461.20001  73.50001 157.49997 340.70025

6.1 Tendencias espaciales

PREC.plt1 #visualizar

#Extraer coordenadas
prec$x<- st_coordinates(prec)[,1]
prec$y<- st_coordinates(prec)[,2]

#Graficos de Dispercion
plot(prec$x, prec$log.prcptcn)

plot(prec$y, prec$log.prcptcn)

#Regresion lineal
tendencia<- lm(log.prcptcn~x+y, prec)
summary(tendencia)

## 
## Call:
## lm(formula = log.prcptcn ~ x + y, data = prec)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.43634 -0.14576  0.06537  0.16496  0.35197 
## 
## Coefficients:
##                 Estimate   Std. Error t value Pr(>|t|)    
## (Intercept)  9.500226097 33.397991421   0.284 0.779498    
## x            0.000012693  0.000002673   4.749 0.000186 ***
## y           -0.000001736  0.000003573  -0.486 0.633228    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2447 on 17 degrees of freedom
## Multiple R-squared:  0.6643, Adjusted R-squared:  0.6248 
## F-statistic: 16.82 on 2 and 17 DF,  p-value: 0.00009354

tendencia.x<- lm(log.prcptcn~x, prec) #Trabajar con tendencia en x
summary(tendencia.x)

## 
## Call:
## lm(formula = log.prcptcn ~ x, data = prec)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.42893 -0.15221  0.05885  0.16619  0.32142 
## 
## Coefficients:
##                 Estimate   Std. Error t value  Pr(>|t|)    
## (Intercept) -6.711589197  1.492181295  -4.498  0.000278 ***
## x            0.000011861  0.000002008   5.906 0.0000137 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2394 on 18 degrees of freedom
## Multiple R-squared:  0.6596, Adjusted R-squared:  0.6407 
## F-statistic: 34.88 on 1 and 18 DF,  p-value: 0.00001368

tendencia.y<- lm(log.prcptcn~y, prec)
summary(tendencia.y)

## 
## Call:
## lm(formula = log.prcptcn ~ y, data = prec)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.57278 -0.20402 -0.00725  0.27740  0.59027 
## 
## Coefficients:
##                  Estimate    Std. Error t value Pr(>|t|)  
## (Intercept) -86.441275347  39.421379944  -2.193   0.0417 *
## y             0.000009133   0.000004066   2.246   0.0375 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3627 on 18 degrees of freedom
## Multiple R-squared:  0.2189, Adjusted R-squared:  0.1755 
## F-statistic: 5.044 on 1 and 18 DF,  p-value: 0.0375

r.sf$x<-st_coordinates(r.sf)[,1]
r.sf$y<-st_coordinates(r.sf)[,2]

r.sf$tendencia<- predict(tendencia.x,r.sf)

r.sf.tendencia<- as(r.sf, "Spatial")
gridded(r.sf.tendencia)<- TRUE

levelplot(brick(r.sf.tendencia)[[3]], margin=FALSE, main="Tendencia",col.regions = levelPal)+ latticeExtra::layer(sp.polygons(as(cuenca, "Spatial"), col='black'))+ latticeExtra::layer(sp.points(as(prec, "Spatial"), pch=20, col='black'))

6.2 Estructura espacial local: autocorrelacion

3 Pasos
- variograma empirico
- variograma teorico o inicial
- ajuste de variograma

La dependencia espacial local significa que cuanto mas cerca esten los puntos en el espacio geografico, mas cerca tambien lo estan en el espacio de atributos. Esto se llama autocorrelacion. Los valores de un atributo pueden estar correlacionados con si mismos, y la fuerza de esta correlacion depende de la distancia de separacion entre puntos. Cada par de puntos, estara separado por una distancia en el espacio geografico y una semivarianza en el espacio de atributos.

alt text

#Calcule cuantos pares de puntos hay en el dataset 
n <- length(prec$log.prcptcn)
n * (n - 1)/2

## [1] 190

#calcule la distancia y semivarianza entre los dos puntos primeros puntos del dataset
st_distance(prec[1:2, ])[2,1]#distancia

## 28441.38 [m]

0.5 * (prec$log.prcptcn[1] - prec$log.prcptcn[2])^2 # gamma semivarianza, unidades log(mm/mes)^2

## [1] 0.07528229

library(automap)
autovar=autofitVariogram(log.prcptcn ~ 1, as(prec, "Spatial"))
summary(autovar)

## Experimental variogram:
##   np      dist      gamma dir.hor dir.ver   id
## 1  7  6432.642 0.03291035       0       0 var1
## 2 14 13716.274 0.03956771       0       0 var1
## 3 17 19150.251 0.06966989       0       0 var1
## 4 21 25460.118 0.10157339       0       0 var1
## 5 22 32506.573 0.09946560       0       0 var1
## 6 30 38802.994 0.15749881       0       0 var1
## 
## Fitted variogram model:
##   model      psill    range kappa
## 1   Nug 0.02684691     0.00     0
## 2   Ste 0.51458242 93707.44     2
## Sums of squares betw. var. model and sample var.[1] 0.00000000002105057

plot(autovar, pch=19, col="black")

autovar=autofitVariogram(log.prcptcn ~ x, as(prec, "Spatial"), miscFitOptions = list(min.np.bin = 5, merge.small.bins = FALSE, cutoff=40000,width=4500))
summary(autovar)

## Experimental variogram:
##   np      dist        gamma dir.hor dir.ver   id
## 1  1  2134.063 0.0001485358       0       0 var1
## 2  1  2806.454 0.0005027130       0       0 var1
## 3  1  5145.038 0.1025893775       0       0 var1
## 4  4  8735.734 0.0603181729       0       0 var1
## 5 14 13716.274 0.0369368034       0       0 var1
## 6 17 19150.251 0.0792982207       0       0 var1
## 7 21 25460.118 0.0661758891       0       0 var1
## 8 22 32506.573 0.0696067143       0       0 var1
## 9 30 38802.994 0.0657149505       0       0 var1
## 
## Fitted variogram model:
##   model      psill    range kappa
## 1   Nug 0.00000000    0.000     0
## 2   Ste 0.06199022 5178.708    10
## Sums of squares betw. var. model and sample var.[1] 0.0000000002688331

plot(autovar, pch=19, col="black")

ve <- variogram(log.prcptcn ~ x, prec, cutoff=40000, width=4500) 
# ve 
plot(ve, plot.numbers = T, asp=1)

show.vgms()

vt <- vgm(psill = 0.065, model = "Sph", range = 15000,nugget = 0.00) 
vt

##   model psill range
## 1   Nug 0.000     0
## 2   Sph 0.065 15000

plot(ve, pl = T, model = vt)

va <- fit.variogram(ve, vt) 
va

##   model      psill    range
## 1   Nug 0.00000000     0.00
## 2   Sph 0.07445553 27251.84

plot(ve, pl = T, model = va)

vt.i <- vgm(nugget=0, model="Sph", range=sqrt(diff(as(prec,"Spatial")@bbox["coords.x1",])^2 +
                 diff(as(prec,"Spatial")@bbox["coords.x2",])^2)/4, 
                 psill=var(prec$log.prcptcn))
#ajuste
va <- fit.variogram(ve, vt.i) 
va

##   model      psill    range
## 1   Nug 0.00000000     0.00
## 2   Sph 0.07445686 27252.85

plot(ve, pl = T, model = va)

# Regla general para determinar par攼㸱metros automaticamente
# Nugget = 0
# partial-sill = varianza total de los datos
# Range = 1/4 de la distancia maxima, diagonal del bounding box

6.3 Interpolacion

kriging ordinario “ordinary” significa que
- (1) la variable es modelada a partir de si misma;
- (2) la media espacial no es conocida a priori, sino estimada de los datos y no existente tendencia espacial.

Kriging Universal significa que existe una tendencia espacial relacionada con las coordendas o la localizacion

Kriging con External Drif significa que se usan otros predictores ademas de las coordenadas.

El resultado de la prediccion proporciona:
- prediction (var1.pred)
- prediction variance (var1.var)

uk <- krige(formula=log.prcptcn ~ x, locations= prec, newdata=r.sf, model = va) 

## [using universal kriging]

uk$pred <- 10^(uk$var1.pred)#volver a valores originales

str(uk)

## Classes 'sf' and 'data.frame':   13054 obs. of  4 variables:
##  $ var1.pred: num  1.5 1.51 1.52 1.53 1.54 ...
##  $ var1.var : num  0.113 0.112 0.111 0.11 0.109 ...
##  $ geometry :sfc_POINT of length 13054; first list element:  'XY' num  682258 9744863
##  $ pred     : num  31.6 32.3 33.1 33.9 34.8 ...
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA
##   ..- attr(*, "names")= chr  "var1.pred" "var1.var" "pred"

uk.sp<- as(uk, "Spatial")
class(uk.sp)

## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"

gridded(uk.sp)<- TRUE
uk.r<- brick(uk.sp)

uk.r.masked<- mask(uk.r,as(cuenca, "Spatial"))

levelplot(uk.r.masked[[3]], margin=FALSE, main="uk",col.regions = levelPal)+ latticeExtra::layer(sp.polygons(as(cuenca, "Spatial"), col='black'))+ latticeExtra::layer(sp.points(as(prec, "Spatial"), pch=20, col='black'))

levelplot(uk.r.masked[[2]], margin=FALSE, main="uk - varianza")+ latticeExtra::layer(sp.polygons(as(cuenca, "Spatial"), col='black'))+ latticeExtra::layer(sp.points(as(prec, "Spatial"), pch=20, col='white'))

uk.cv <- krige.cv(formula=log.prcptcn ~ x, locations= prec, model = va)
uk.cv[1:5,]

## Simple feature collection with 5 features and 6 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 705090 ymin: 9662396 xmax: 797962 ymax: 9713142
## epsg (SRID):    NA
## proj4string:    NA
##   var1.pred   var1.var observed    residual      zscore fold
## 1  1.849873 0.06799553 2.011993  0.16211994  0.62172204    1
## 2  2.133035 0.08164801 2.400020  0.26698508  0.93435987    2
## 3  2.729451 0.11842431 2.663889 -0.06556134 -0.19051421    3
## 4  1.789338 0.05699885 1.866287  0.07694940  0.32230881    4
## 5  2.186670 0.07514738 2.197280  0.01061078  0.03870707    5
##                 geometry
## 1 POINT (717092 9662396)
## 2 POINT (745529 9662895)
## 3 POINT (797962 9699430)
## 4 POINT (705090 9665080)
## 5 POINT (748262 9713142)

uk.cv$obs <- 10^(uk.cv$observed)

(uk.mse<-mse(10^(uk.cv$var1.pred), uk.cv$obs))

## [1] 3402.075

(uk.rmse<-rmse(10^(uk.cv$var1.pred), uk.cv$obs))

## [1] 58.32731

(uk.cor<-cor(10^(uk.cv$var1.pred), uk.cv$obs))

## [1] 0.8993445

(uk.bias<-bias(10^(uk.cv$var1.pred), uk.cv$obs))

## [1] -7.214766

plot(10^(uk.cv$var1.pred), uk.cv$obs)

ggplot(uk.cv, aes(obs, y = var1.pred, color = "red")) +
  geom_abline(colour = "#CCCCCC", linetype = 2) +
  geom_point(aes(y = 10^var1.pred, col = "UK")) +
  geom_smooth(aes(y = 10^var1.pred, col = "UK"), method = lm, se = FALSE) +
  xlab("observaciones")+
  ylab("predicciones")

7.0.1 Media constante

• simple kriging = media conocida
• ordinary kriging = media desconocida
• indicator krigning = datos categoricos / above - under threshold

7.0.2 Tendencia en la media

• co-kriging = covariables, solo disponible en ciertos puntos
• universal kriging = tendencia en funcion de las coordenadas, covariables en todos los puntos de prediccion. deterministic part is modeled as a function of coordinates -kriging with external drift= tendencia en funcion de covariables, covariables en todos los puntos de prediccion. both components are predicted simultaneously.
• Regresion kriging = se maneja por separado la componente deterministica y estocastica. the deterministic (regression) and stochastic (kriging) predictions are done separately