Taller: Distancia Inversa Ponderada

EstimaciĆ³n de un valor usando IDW

El mĆ©todo de la distancia inversa ponderada o IDW por sus siglas en inglĆ©s, es un interpolador exacto, ya que respeta Ć­ntegramente los valores de entrada, asumiendoles como observaciones veraces del fenĆ³meno a estudiar.

IDW puede expresarse de diferentes maneras, para este caso se desarrolla como

\[\hat{z}(s_0) = \frac{\sum_{i=1}^{n} z(s_{i})d_{i0}^{-r}}{\sum_{i=1}^{n}d_{i0}^{-r}}\] En donde:

  • \(\hat{z}(s_0)\) es el valor interpolado en la ubicaciĆ³n \(s_{0}\),
  • \(z(s_{i})\) es el valor conocido en la ubicaciĆ³n \(s_{i}\),
  • \(d_{i0}\) es la distancia entre la ubicaciĆ³n \(s_{i}\) y la ubicaciĆ³n \(s_{0}\),
  • \(r\) es el parĆ”metro de potencia que controla el peso de la distancia,
  • \(n\) es el nĆŗmero de puntos conocidos.

La aplicaciĆ³n del mĆ©todo IDW resulta sencilla, dado que no requiere de la inspecciĆ³n de supuestos, su implementaciĆ³n requiere solo la ubicaciĆ³n de los valores en un espacio coordenado, sobre los cuales se realiza el cĆ”lculo de la distancia a partir del cual se estima el peso optimo, elevĆ”ndole al inverso del parĆ”metro de potencia \(r\).

Datos espaciales

Para el ejemplo se considerarĆ” la informaciĆ³n de temperatura para 6 estaciones climatolĆ³gicas

# Crear el data.frame con los datos proporcionados
df <- data.frame(
  Nombre = c("DoƱa Juana", "El Dorado", "Las Mercedes", "La Florida", "El Japon", "Pasca"),
  Elevacion = c(2700, 2546, 810, 1915, 280, 2256), # NA para el valor faltante
  Temp = c(9.3,14,25.1,16.9,26.1,15.8),
  x = c(-74.16666667, -74.15, -74.52661111, -74.43763889, -73.30130556, -74.31175),
  y = c(4.5, 4.7, 4.581888889, 4.770888889, 4.376972222, 4.310111111),
  Este = c(4870620,4872504,4830714,4840625,4966583,4854491),
  Norte = c(2055349,2077449,2064474,2085344,2041659,2034390))
# Imprimir las primeras 5 filas
knitr::kable(df,
             caption = "InformaciĆ³n de estaciones climatolĆ³gicas",
             align=c('c','c','c','c','c','c'),col.names = c('EstaciĆ³n','ElevaciĆ³n','Temperatura','X','Y','Este','Norte'))
InformaciĆ³n de estaciones climatolĆ³gicas
EstaciĆ³n ElevaciĆ³n Temperatura X Y Este Norte
DoƱa Juana 2700 9.3 -74.16667 4.500000 4870620 2055349
El Dorado 2546 14.0 -74.15000 4.700000 4872504 2077449
Las Mercedes 810 25.1 -74.52661 4.581889 4830714 2064474
La Florida 1915 16.9 -74.43764 4.770889 4840625 2085344
El Japon 280 26.1 -73.30131 4.376972 4966583 2041659
Pasca 2256 15.8 -74.31175 4.310111 4854491 2034390

Estos valores se distribuyen en el espacio de la siguiente manera

# Crear paleta de colores
fun_color_range <- colorRampPalette(c("blue2", "red"))  # Create color generating 
# Crear el grƔfico con ggplot2 y ggspatial
temp_plot=ggplot(df, aes(x = x, y = y)) +
  geom_point(aes(color = Temp), size = 4) +   # Usar color para la variable Temp
  scale_colour_gradientn(colors = fun_color_range(16))+
  labs(title = "Temperaturas por EstaciĆ³n",
       x = "Longitud", 
       y = "Latitud", 
       color = "Temperatura (Ā°C)") +
  theme_gray() +
  theme(legend.position = "none")+
    # Agregar la flecha norte
  annotation_north_arrow(location = "ur", which_north = "true", 
                         pad_x = unit(.05, "in"), pad_y = unit(3, "in"),
                         style = north_arrow_fancy_orienteering)
# Convertir el grƔfico de ggplot en interactivo con plotly
ggplotly(temp_plot)

A partir de estos datos se busca estimar cuĆ”l serĆ­a el valor de la temperatura para la estaciĆ³n Sutatausa ubicada en las coordenadas \(X = -73.85\) y \(Y = 5.25\)

# Crear objeto basado en coordenadas planas
df_planas = df
coordinates(df) <- ~x+y
coordinates(df_planas) <- ~Este+Norte

# PosiciĆ³n de valor a estimar
new_data <- data.frame(x = -73.85, y = 5.25) # GeogrƔficas
new_data_planas = data.frame(Este = 4905845, Norte = 2138180) # Proyectadas

#Convertir en objeto sp
coordinates(new_data) <- ~x + y  # Definir las coordenadas
coordinates(new_data_planas) <- ~Este + Norte  # Definir las coordenadas

Estimar valores usando IDW

Para este caso se estimarĆ” el valor de temperatura en la ubicaciĆ³n seleccionada usando diferentes valores de \(r\), inicialmente serĆ”n considerados los siguientes:

  • 0, 2, 5, 10, 80

Los resultados obtenidos se muestran separados a continuaciĆ³n, diferenciando por el objeto espacializado usando las coordenadas geogrĆ”ficas y proyectadas, en este caso el MAGNA-SIRGAS Origen-Nacional.

Estimaciones de valor desconocido

Usando los valores de \(r\) mencionados se obtienen los siguientes resultados

# Estimar en posiciĆ³n dada
# Realizar la interpolaciĆ³n IDW
idp = c(0,2,5,10,80)

for (p in idp) {
  #Estimar con geogrƔficas
  estimado = idw(formula = Temp ~ 1, 
                 locations = df, 
                 newdata = new_data, 
                 idp = p,
                 maxdist=Inf)["var1.pred"]$var1.pred
  # Estimar con planas
  estimado_planas =idw(formula = Temp ~ 1, 
                       locations = df_planas, 
                       newdata = new_data_planas, 
                       idp = p,
                       maxdist=Inf)["var1.pred"]$var1.pred
  # Ver resultados en geogrƔficas
  print(paste('Con un r =',p,'Se obtiene: ',estimado,' Estimando en geogrƔficas'))
  # Ver resultados con planas
  print(paste('Con un r =',p,'Se obtiene: ',estimado_planas,' Estimando en proyectadas'))
  }
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [1] "Con un r = 0 Se obtiene:  17.8666666666667  Estimando en geogrƔficas"
## [1] "Con un r = 0 Se obtiene:  17.8666666666667  Estimando en proyectadas"
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [1] "Con un r = 2 Se obtiene:  16.6559304323074  Estimando en geogrƔficas"
## [1] "Con un r = 2 Se obtiene:  16.6531270456695  Estimando en proyectadas"
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [1] "Con un r = 5 Se obtiene:  15.2207011697664  Estimando en geogrƔficas"
## [1] "Con un r = 5 Se obtiene:  15.2149125490283  Estimando en proyectadas"
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [1] "Con un r = 10 Se obtiene:  14.2822067689526  Estimando en geogrƔficas"
## [1] "Con un r = 10 Se obtiene:  14.27697316753  Estimando en proyectadas"
## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
## [1] "Con un r = 80 Se obtiene:  14.0000006775029  Estimando en geogrƔficas"
## [1] "Con un r = 80 Se obtiene:  NaN  Estimando en proyectadas"

AdemĆ”s de las sutiles diferencias en la estimaciĆ³n resultante, llama la atenciĆ³n que, usando un objeto espacializado con coordenadas proyectadas, no es posible obtener un valor de temperatura usando un \(r\) de 80.

Interpretando resultados

A partir de una iteraciĆ³n que considere valores de \(r\) entre 0 y 100 es posible identificar cĆ³mo estos afectan la estimaciĆ³n generada considerando la distancia entre los puntos, con ello es posible visualizar que la estimaciĆ³n usando un valor \(r = 0\) es igual a la media de los valores observados, mientras que con valores \(r\) muy grandes, la influencia del punto mĆ”s cercano se vuelve mĆ”s fuerte y condiciona el valor resultante, generando asĆ­ que el mĆ©todo se comporte mĆ”s como un estimador local que como uno global

# Graficar
# Scatter plot by group
iteraciones = ggplot(estimados_idw, aes(x =  p_value , y = Estimado, color = Tipo)) +
  geom_point(aes(shape=Tipo, color=Tipo),size=2.5)+
  labs(title = "EstimaciĆ³n usando coordenadas proyectadas y geogrĆ”ficas", 
       subtitle = "IDP Vs Temperatura estiamada",
       x='Potencia de ponderaciĆ³n de distancia inversa',
       y='Temperatura CĀ°')+
  theme_grey()+
  theme(legend.position = "top")
# Convertir el grƔfico de ggplot en interactivo con plotly
ggplotly(iteraciones)

EstimaciĆ³n de superficies

Una vez observada la influencia que tienen tanto los diferentes valores de \(r\) en la estimaciĆ³n con IDW, como tambiĆ©n el sistema de coordenadas implemetado, es posible generar superficies continuas de interpolaciĆ³n, incluyendo ahora el valor observado en la posiciĆ³n que anteriormente fue estimada, la estaciĆ³n Sutatausa, con registro de temperatura media de 13.6 CĀ°.

# Crear el data.frame con los datos proporcionados
df <- data.frame(
  Nombre = c("DoƱa Juana", "El Dorado", "Las Mercedes", "La Florida", "El Japon", "Pasca","Sutatausa"),
  Elevacion = c(2700, 2546, 810, 1915, 280, 2256,3600), # NA para el valor faltante
  Temp = c(9.3,14,25.1,16.9,26.1,15.8,13.6),
  x = c(-74.16666667, -74.15, -74.52661111, -74.43763889, -73.30130556, -74.31175,-73.85),
  y = c(4.5, 4.7, 4.581888889, 4.770888889, 4.376972222, 4.310111111,5.25),
  Este = c(4870620,4872504,4830714,4840625,4966583,4854491,4905845),
  Norte = c(2055349,2077449,2064474,2085344,2041659,2034390,2138180))
# Espacializar
coordinates(df) <- ~Este + Norte

Una estimaciĆ³n en superficie continua tambiĆ©n permite aplicar una validaciĆ³n cruzada, de tal forma que sea posible calcular un error medio con respecto a los valores observados para la variable en cuestiĆ³n.

# Crear una reticula de pixeles vacĆ­a
grd              <- as.data.frame(spsample(df, "regular", n=50000))
names(grd)       <- c("X", "Y")
coordinates(grd) <- c("X", "Y")
gridded(grd)     <- TRUE  # Crear Objeti Spatialpixel
fullgrid(grd)    <- TRUE  # Crear objeto SpatialGrid

# crear objeto nuevo almacenando la informaciĆ³n interpolada

# Estimar valores usando un r = 2
temp.idw <- idw(Temp ~ 1, df, newdata=grd, idp = 2.0)
## [inverse distance weighted interpolation]
# Convertir a raster obejto grid interpolado
r_temp <- rast(temp.idw)

# Graficar
# Convertir raster a dataframe para ggplot
r_df <- as.data.frame(r_temp, xy = TRUE)
# Convertir a objeto sf
df_sf <- st_as_sf(df, crs = 9377)
# Graficar el raster y los puntos usando ggplot
temp_idw=ggplot() +
  geom_raster(data = r_df, aes(x = x, y = y, fill = var1.pred)) + # Graficar raster
  scale_fill_distiller(palette = "RdBu", direction = -1, name = "Temperatura CĀ°") + # Color scale
  geom_sf(data = df_sf, color = "black", size = 0.5) +  # Graficar puntos sobre raster
  labs(title = "Temperatura estimada",
       x = "Este", y = "Norte") +
  theme_minimal() +
  # Eliminar leyenda fuera del grƔfico
  theme(legend.position = "none")
# Convertir a dinƔmico
ggplotly(temp_idw)

Siguiendo el mismo procedimiento y utilizando un valor \(r = 2\) sugerido por la literatura, es posible generar tambiĆ©n una superficie de interpolaciĆ³n a partir de los valores de altura.

# Crear una reticula de pixeles vacĆ­a
grd              <- as.data.frame(spsample(df, "regular", n=50000))
names(grd)       <- c("X", "Y")
coordinates(grd) <- c("X", "Y")
gridded(grd)     <- TRUE  # Crear Objeti Spatialpixel
fullgrid(grd)    <- TRUE  # Crear objeto SpatialGrid

# crear objeto nuevo almacenando la informaciĆ³n interpolada

# Estimar valores usando un r = 2
elev.idw <- idw(Elevacion ~ 1, df, newdata=grd, idp = 2.0)
## [inverse distance weighted interpolation]
# Convertir a raster obejto grid interpolado
r_elev <- rast(elev.idw)
# Graficar
# Convertir raster a dataframe para ggplot
r_df <- as.data.frame(r_elev, xy = TRUE)
# Convertir a objeto sf
df_sf <- st_as_sf(df, crs = 9377)
# Graficar el raster y los puntos usando ggplot
temp_idw=ggplot() +
  geom_raster(data = r_df, aes(x = x, y = y, fill = var1.pred)) + # Graficar raster
  scale_fill_distiller(palette = "BrBG", direction = -1, name = "Altura sobe el nivel del mar") + # Color scale
  geom_sf(data = df_sf, color = "black", size = 0.5) +  # Graficar puntos sobre raster
  labs(title = "Altura estimada",
       x = "Longitude", y = "Latitude") +
  theme_minimal() +
  # Eliminar leyenda fuera del grƔfico
  theme(legend.position = "none")
# Convertir a dinƔmico
ggplotly(temp_idw)

Las superficies generadas presentan grandes similitudes, fenĆ³meno que se debe en gran parte a la fuerte relaciĆ³n existente entre Temperatura y Altitud. De igual forma es posible observa como IDW genera contornos alrededor de los valores observados, evidenciando su caracterĆ­stica de estimador exacto al preservar los valores de entrada. Esto puede ser un problema en casos en donde la concentraciĆ³n de puntos en un espacio reducido, es amplia, generando transiciones que no necesariamente corresponden a la realidad.

# Add the regression line
ggplot(df_sf, aes(x=Elevacion, y=Temp)) + 
  geom_point()+
  labs(title = "Temperatura CĀ° Vs ElevaciĆ³n", 
       x='Altura m.s.n.m',
       y='Temperatura CĀ°')+
  geom_smooth(method=lm)
## `geom_smooth()` using formula = 'y ~ x'

Para IDW es necesario considerar que los valores estimados se mueven en el rango de los valores muestreados, lo cual hace necesario que el muestreo que alimenta el modelo tenga el mejor cubrimiento espacial y temƔtico posible.

ValidaciĆ³n cruzada

Validar los resultados es un paso importante dentro de las estimaciones espaciales, en este caso se hace uso de la validaciĆ³n cruzada para obtener una idea sobre la eficiencia del mĆ©todo a la hora de estimar valores que previamente se conocen. Confrontando el valor observados con el valor estimado. Esto permite resumir en un indicador numĆ©rico cuĆ”l es la eficacia del mĆ©todo para modelar la variable.

# Graficar las diferencias
# Crear un dataframe con las variables
data <- data.frame(Observed = idw.temp.model$observed, Predicted = round(idw.temp.model$var1.pred,1), Error =idw.temp.model$var1.pred -  idw.temp.model$observed)

# Crear el grĆ”fico de dispersiĆ³n usando ggplot2
ggplot(data, aes(x =  Observed , y =Predicted)) +
  ylim(0,27)+
  geom_point(color = rgb(0, 0, 0, 0.5), size = 3) +  # Puntos
  geom_smooth(method = "lm", color = "red", linetype = "dashed", se = FALSE) +  # LĆ­nea de regresiĆ³n
  geom_abline(intercept = 4.4966, slope = 0.7763 , linetype = "solid", color = "blue") +  # lm omitiendo valores extremos
  labs(x = "Observado", y = "Estimado", title = "Valor Observado Vs Estimado") +
  theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'

En este caso se identifica una correlaciĆ³n baja entre los valores estimados y los valores observados, sugiriendo la necesidad de densificar la muestra con el objetivo de mejorar la eficiencia del modelo espacial.