Métodos determinísticos

Análisis espacial con R

Dr. Francisco Zambrano

Preparar datos para interpolación

library(tidyverse)
library(sf)
library(tmap)
region <- read_rds('data/regs_utm_sf.rds')
data_temp <- read_rds('data/data_temp_enero.rds')[c('ema','temp_enero')]
set.seed(432)
data_temp_mod <- data_temp |> sample_frac(0.8)
data_temp_val <- data_temp |> filter(!(ema %in% data_temp_mod$ema))

calc_met <- function(ras,object_sf,var = 'temp_enero') {
  df_val <- terra::extract(ras,terra::vect(object_sf)) |> 
    cbind(object_sf[,var]) 
  rmse <- function(obs,pred) sqrt(sum((obs-pred)^2)/length(obs))
  rmse_v <- rmse(df_val[,3],df_val[,2])
    r2 <- cor(df_val[,3],df_val[,2])^2
  c('rmse' = rmse_v,'rsq' = r2)
}


tm_shape(region) + 
  tm_borders() +
  tm_shape(data_temp) +
  tm_dots(col='temp_enero')

Preparar datos para interpolación

Grilla

En R, espacializar temperatura promedio del mes de enero 2023.

library(stars)
library(terra)
bb <- st_bbox(region)
grilla <- rast(xmin = bb$xmin,xmax = bb$xmax,
               ymin = bb$ymin,ymax = bb$ymax, 
               res=10000,crs="EPSG:32719") 
grilla_st <- grilla |>   st_as_stars() #como clase stars

grilla |> as.polygons() |> plot()

Voronoi (Polígonos de Thiessen)

Zonas de influencia más cercana a cada punto. Caso particular de vecino más próximo con k=1.

temp_vor <- voronoi(vect(data_temp_mod)) |> 
  terra::rasterize(grilla,'temp_enero') |> 
  mask(vect(region)) 
(met_vor <- calc_met(temp_vor,data_temp_val))

     rmse       rsq 
2.1243421 0.3790889 

Vecino más próximo

Vecino más póximo con k=1.

kNN Interactivo

library(gstat)
temp_nn <- idw(temp_enero~1, data_temp_mod, grilla_st,nmax=1)

[inverse distance weighted interpolation]

temp_nn <- temp_nn |> rast() |> mask(vect(region)) |> trim() 
(met_nn1 <- calc_met(temp_nn[[1]],data_temp_val))

     rmse       rsq 
2.1243421 0.3790889 

Vecino más próximo

Vecino más póximo con k=5.

library(gstat)
temp_nn5 <- idw(temp_enero~1, data_temp_mod, grilla_st,nmax=5)

[inverse distance weighted interpolation]

temp_nn5 <- temp_nn5 |> rast() |> mask(vect(region)) |> trim() 
(met_nn5 <- calc_met(temp_nn5[[1]],data_temp_val))

    rmse      rsq 
1.620426 0.516043 

Vecino más próximo

Vecino más póximo con k=10.

library(gstat)
temp_nn10 <- idw(temp_enero~1, data_temp_mod, grilla_st,nmax=10)

[inverse distance weighted interpolation]

temp_nn10 <- temp_nn10 |> rast() |> mask(vect(region)) |> trim() 
(met_nn10 <- calc_met(temp_nn10[[1]],data_temp_val))

     rmse       rsq 
1.5513405 0.5455748 

Inverso de la Distancia

• Quizas uno de los métodos más antiguos de predicción espacial.

\[\hat z(s_0) = \sum_{i=1}^n \lambda_i(s_0)\cdot z(s_i)\]

Notación de matriz

\[\hat z(s_0) = \lambda_0^T \cdot z\]

Los pesos se determinan por el inverso de la distancia

\[\lambda(s_0) = \frac{\frac{1}{d^\beta (s_0,S_i)}}{\sum_{i=1}^n \frac{1}{d^\beta (s_0,s_i)}}\] \(\beta >0\)

Inverso de la Distancia

Inverso de la Distancia

En R, espacializar temperatura promedio del mes de enero 2023.

temp_idw <- idw(temp_enero~1, data_temp_mod, grilla_st)

[inverse distance weighted interpolation]

temp_idw <- temp_idw |> rast() |> mask(vect(region)) |> trim() 
met_idw <- calc_met(temp_idw[[1]],data_temp_val)

tm_shape(temp_idw[[1]]) +
  tm_raster(style = 'cont',palette = viridis::viridis(20)) +
  tm_shape(region) + 
  tm_borders() +
  tm_shape(data_temp) +
  tm_dots(col='temp_enero',title = "Temperatura (°C)",style = 'jenks',palette = viridis::viridis(10))

Trend surface / regresión en coordenadas

• Asumimos que los valores de la variable objetivo (geográficamente ubicada) es función de sus coordenadas.

• Podemos determinar sus valores encontrando una función a través de los puntos observados.

• Se basa en el siguiente modelo:

\[Z(s) = f(x,y) + \epsilon\]

Trend surface / regresión en coordenadas

En R

\[Temp=\beta_0+\beta_1\cdot x+\beta_2\cdot y\]

data_temp_mod <- cbind(data_temp_mod,st_coordinates(data_temp_mod))
mod.lm1 <- lm(temp_enero ~ X + Y, data=data_temp_mod)

\[Temp = \beta_0+\beta_1\cdot x + \beta_1 \cdot y + \beta_2\cdot x^2+\beta_3\cdot y^2 + \beta_4\cdot x\cdot y\]

mod.lm2 <- lm(temp_enero ~ X + Y + I(X^2) + I(Y^2) + I(X*Y), data=data_temp_mod)

Trend surface / regresión en coordenadas

\[Temp=\beta_0+\beta_1\cdot x+\beta_2\cdot y\]

summary(mod.lm1)$r.squared

[1] 0.1376663

\[Temp = \beta_0+\beta_1\cdot x + \beta_1 \cdot y + \beta_2\cdot x^2+\beta_3\cdot y^2 + \beta_4\cdot x\cdot y\]

summary(mod.lm2)$adj.r.squared

[1] 0.3496565

summary(mod.lm2)$r.squared

[1] 0.3646414

Trend surface / regresión en coordenadas

temp_ts <- c(grilla,grilla)
xy <- xyFromCell(grilla,1:ncell(grilla)) |> 
  as_tibble() |> set_names(c('X','Y'))
values(temp_ts[[1]])<-  predict(mod.lm1, newdata=xy)

values(temp_ts[[2]]) <- predict(mod.lm1, newdata=xy,se.fit=TRUE,inerval='prediction')$se.fit
names(temp_ts) <- c('temperatura','error_estandar')
temp <- mask(temp_ts, vect(region))
(met_ts1 <- calc_met(temp_ts[[1]],data_temp_val))

     rmse       rsq 
1.9647131 0.2579367 

Trend surface / regresión en coordenadas

temp_ts2 <- c(grilla,grilla)
values(temp_ts2[[1]])<-  predict(mod.lm2, newdata=xy)

values(temp_ts2[[2]]) <- predict(mod.lm2, newdata=xy,se.fit=TRUE,inerval='prediction')$se.fit
names(temp_ts2) <- c('temperatura','error_estandar')
temp_ts2 <- mask(temp_ts2, vect(region))
(met_ts2 <- calc_met(temp_ts2[[1]],data_temp_val))

     rmse       rsq 
1.8479598 0.3549047 

Evaluación de los modelos