Introducción

Este trabajo forma parte del módulo Estadística Espacial en el posgrado en Big Data e Inteligencia Territorial de FLACSO Argentina a cargo de Priscilla Minotti.

Se abordará un data set sobre puntos de elevación (originalmente medidos en pies pero convertidos a metros) de edificios en la ciudad de Nueva York con distintos métodos de interpolación espacial.

Los datos originales están disponibles acá mientras que la muestra con la que trabajamos junto con el proyecto del trabajo acá.

Paquetes

Cargamos los paquetes que vamos a utilizar:

library(gstat)
library(tidyverse)
library(sf)
library(sp)
library(rgdal)
library(raster)
library(osmdata)
library(leaflet)
library(tmap)

Ahora vamos a explicar paso por paso cómo es el procedimiento desde la carga de datos al procesamiento, las visualizaciones, los métodos de interpolación en sí y los resultados.

Carga de datos

Cargamos el data set con la muestra de los puntos de elevación que vamos a usar para interpolar.

elevaciones <- st_read("elevaciones_NY_muestra.shp") %>%
  st_as_sf() 
## Reading layer `elevaciones_NY_muestra' from data source `/Users/usernamemateo/Library/Mobile Documents/com~apple~CloudDocs/R/geostat/elevaciones_NY_muestra.shp' using driver `ESRI Shapefile'
## Simple feature collection with 10849 features and 2 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 913410.9 ymin: 121129.2 xmax: 1067130 ymax: 271166.4
## Projected CRS: NAD83 / New York Long Island (ftUS)

Para entender mejor los puntos de elevación vamos a usar summary()

summary(elevaciones$elevacion)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   3.682  14.914  22.315  25.085  31.188 430.803

Como vemos, el punto más de medición de altura de edificios mas bajo está elevado a 3.68 metros, en promedio la altura de los edificios es de 25.085 metros y la mediana es 22.315 metros. La altura máxima es de algo más de 430 metros.

Visualizaciones básicas

Para entender mejor los datos vamos a graficar los puntos de elevación ubicados en el mapa. Para eso vamos a usar tmap. Lo configuramos.

tmap_mode("view")

Y ahora, graficamos los puntos con el mapa base:

qtm(elevaciones, scale = 0.5)

Ahora vamos a ver específicamente cómo de distribuyen los puntos de elevación de nuestro data set y cual es su altura:

ggplot() + 
  geom_sf(data = elevaciones, aes(color = elevacion), size = 0.5) + 
  ggtitle("Puntos de elevación: edificios en la ciudad de Nueva York") +
  scale_color_viridis_c(name = "Altura (mts)", direction = -1, option = "magma") + 
  theme_dark()

Como podemos observar, en nuestro dataset hay unos pocos edificios con altura mayor a 300 metros. La mayoría los edificios tienen alturas menores de 150 metros.

Límites y Áreas de Interés

En este apartado vamos a ver tres métodos para definir el área de trabajo o de interés (AOI) que nos resulta útil para definir el área de análisis espacial a partir de nuestros datos.

Bounding box (bbox)

Bounding box (o simplemente bbox) es un área definida por dos latitudes y dos longitudes.En este caso, usamos bbox que nos sirve para setear los límites (en forma de caja).

elevaciones_bbox <- st_as_sfc(st_bbox(elevaciones))

Buffer

Cuando usamos la función st_buffer(), ésta nos devuelve una geometría que representa todos los puntos cuya distancia es menor o igual a la distancia indicada. En este caso vamos a usar como distancia 2000 metros para crear un polígono usable para interpolar.

elevaciones_buffer <- elevaciones %>% 
  st_buffer(dist= 2000) %>% 
  st_geometry() %>% 
  st_simplify() %>% 
  st_union()

elevaciones_buffer
## Geometry set for 1 feature 
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 911410.9 ymin: 119129.2 xmax: 1069130 ymax: 273166.4
## Projected CRS: NAD83 / New York Long Island (ftUS)
## MULTIPOLYGON (((919078.2 120248, 919065.4 12022...

Polígono máximo convexo o convex hull

Convex hull es el polígono convexo más pequeño que incluye a todos los puntos incluidos en un data set. Vamos a calcular el polígono máximo convexo de nuestros puntos de elevación:

elevaciones_convex_hull <- elevaciones %>% st_union()%>% st_convex_hull()
elevaciones_convex_hull
## Geometry set for 1 feature 
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 913410.9 ymin: 121129.2 xmax: 1067130 ymax: 271166.4
## Projected CRS: NAD83 / New York Long Island (ftUS)
## POLYGON ((917283.4 121129.2, 916440.4 121295.2,...

Límites políticos y administrativos de la ciudad

Cargamos un shapefile con los límiter políticos y administrativos de la ciudad (sin agua). Los límites están disponibles acá.

elevaciones_limites <- st_read("nybb.shp") %>%
  st_geometry() %>% 
  st_simplify() %>% 
  st_union()
## Reading layer `nybb' from data source `/Users/usernamemateo/Library/Mobile Documents/com~apple~CloudDocs/R/geostat/nybb.shp' using driver `ESRI Shapefile'
## Simple feature collection with 5 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 913175.1 ymin: 120121.9 xmax: 1067383 ymax: 272844.3
## Projected CRS: NAD83 / New York Long Island (ftUS)

Límites en el mapa

Ahora vamos a ver cada tipo de límite en el mapa con los puntos de elevación:

qtm(elevaciones, scale = 0.5, title = "Límites de la ciudad en el mapa") + 
  #en rojo los bordes de la bounding box
  qtm(elevaciones_bbox, fill = NULL, borders="red") +
  #en azul los bordes de buffer
  qtm(elevaciones_buffer, fill = NULL, borders = "blue") +
  #en verde los bordes del poligono convexo
  qtm(elevaciones_convex_hull, fill = NULL, borders = "green") +
  #en negro los bordes de los límites políticos y administrativos
  qtm(elevaciones_limites, fill = NULL, borders = "black")

De ahora en adelante vamos a trabajar con los límites políticos y administrativos.

Interpolación espacial

Vamos a desarrollar tres métodos de interpolación: Polígonos de Voronoi o Thiessen, Distancia Inversa Ponderada (de ahora en más IDW) y Kriging.

Polígonos de Voronoi o Thiessen

Los puntos de elevación conocidos son utilizados como puntos generadores para dividir nuestra área de interés en tantos polígonos o teselas como puntos conocidos. El espacio interno de estos polígonos contiene todos los puntos cuya distancia al punto generador es menor que la distancia a cualquier otro punto externo. La división es a través de bisectores perpendiculares (líneas) que se forman al calcular las distancias entre puntos.

voronoi <- elevaciones %>% 
          st_geometry() %>%
          st_union() %>%
          st_voronoi() %>%
          st_collection_extract()%>%
          st_sfc(crs = st_crs(elevaciones)) %>% 
          st_sf() %>% 
          st_join(elevaciones) %>%
          st_intersection(elevaciones_limites)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

Graficamos:

qtm(voronoi, fill = "elevacion", fill.palette = "-viridis", borders = "black", scale = 0.1,  title = "Polígonos de Voronoi: altura de edificios en Nueva York")

Como podemos ver, se generaron tantos polígonos como puntos tiene nuestro dataset. El polígono con el edificio más alto está ubicado en la isla de Manhattan.

Grilla de puntos

Para utilizar los dos métodos de interpolación restantes necesitamos una grilla de puntos. Vamos a usar el dataset de puntos de elevación como base y la grilla estará contenida en una bounding box. Para eso:

grilla_bbox <-  st_make_grid(elevaciones, cellsize = c(700, 700), what = "centers")

Nuestra grilla tiene 47300 puntos. Ahora vamos usar la función st_intersection() para quedarnos solo con los valores que se encuentran dentro de los límites de la ciudad.

elevaciones_grilla <- st_intersection(elevaciones_limites, grilla_bbox)

Vamos a renombrar la columna x como geometry.

elevaciones_grilla <- elevaciones_grilla %>%
  st_as_sf() %>%
  rename(geometry=x)

La grilla resultante tiene 17160 puntos ubicados dentro de los límites de la ciudad. Vamos a ver la grilla de puntos en el mapa:

qtm(elevaciones_grilla, scale = 0.5)

A partir de ahora, usamos esta grilla de puntos para interpolar con IDW y Kriging.

IDW

En la Distancia Inversa Ponderada cada punto tiene influencia local que disminuye con la distancia. El valor a predecir se calcula usando una combinación lineal ponderada de los vectores muestreados. A mayor proximidad al lugar a predecir, mayor peso de los puntos vecinos.

La potencia de la distancia inversa por defecto es uno (1). Vamos a calcular con esta configuración básica:

elevaciones_idw <- gstat::idw(elevacion ~ 1, elevaciones, elevaciones_grilla)
## [inverse distance weighted interpolation]

Algunas estadísticas sobre la altura de las predicciones:

summary(elevaciones_idw$var1.pred)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   5.478  18.985  24.646  25.994  30.340 178.847

Como se puede observar, el valor máximo de las predicciones de altura de los edificios es de 178.847 metros. La altura promedio de los edificios es de casi 26 metros.

Graficamos con ggplot:

ggplot() + 
  geom_sf(data = elevaciones_idw, aes(color = var1.pred), size = 0.35) + 
  ggtitle("IDW: puntos de elevación", subtitle = "Altura de edificios en la ciudad de Nueva York") +
  scale_color_viridis_c(name = "altura (mts)", direction = -1, option = "magma") + 
  theme_dark()

En este caso, el rango de valores es de 0 a 180 metros.

Kriging

Este método calcula la estructura de variación espacial a partir de muestras dispersas en el espacio. Genera una superficie estimada, es decir una variable regionalizada y una superficie de error.

Kriging calcula la variación de una variable asumiendo que la media no tiene grandes variaciones, entonces necesita conocer la varianza de los datos. Para conocer esto se calcula la semivarianza que describe el patron de variabilidad espacial entre pares de puntos a distintos intervalos de distancia.

Semivariogramas

Vamos a calcular el semivariograma muestreal y el teórico o ajustado. El objetivo del semivariograma es mostrar que a partir de cierto valor de distancia, el estimador de variación se nivela (no varía más o varía pero no en función de la distancia entre pares de puntos).

La distancia a la que el modelo se estabiliza se denomina rango o range. Los puntos muestrales separados por distancias más cortas que el rango están autocorrelacionadas espacialmente, mientras que los más alejados que el rango no lo están.

El valor de variación máxima se llama meseta o sill (C1).

Variograma muestreal

Los puntos resultantes del cálculo de la semivarianza se grafican en relación con la distancia en los que conocemos como semivariograma muestral o experimental.

#calculamos el variograma muestreal
lnz_vgm <- variogram(log(elevacion)~1, elevaciones)

#graficamos
plot(lnz_vgm, main ="Variograma muestreal")

Como se puede ver, las distancias (en metros) son grandes y la semivarianza es baja, pero sigue la distribución esperada. El rango de aproximadamente 20.000 metros.

Variograma teórico o ajustado

Los puntos del semivariograma muestral se ajustan según un modelo matemático (en este caso vamos a usar Sph) y el resultado es el semivariograma teórico o ajustado.

lnz_fit_sph <- fit.variogram(lnz_vgm, model=vgm("Sph"))
lnz_fit_sph
##   model      psill    range
## 1   Nug 0.03026392     0.00
## 2   Sph 0.23049692 26640.97

Como podemos ver: al ajustar, el rango se encuentra en 26640.97

plot(lnz_vgm, lnz_fit_sph, 
     main ="Variograma teórico o ajustado",
     sub = "Método: Sph")

Kriging ordinario

El kriging usa los valores muestreales como el log(elevacion) de nuestro dataset original para interpolar cada punto de la grilla elevaciones_grilla. Usa el modelo ajustado, que en este caso es el esférico ya que recién usamos sph.

lnz_pred <- gstat::krige(log(elevacion)~1, elevaciones, elevaciones_grilla, lnz_fit_sph)
## [using ordinary kriging]

Ahora vamos a graficar:

ggplot() + 
  geom_sf(data = lnz_pred, aes(color = exp(var1.pred)), size = 0.35) +
  ggtitle("Kriging: puntos de elevación", subtitle = "Altura de edificios en la ciudad de Nueva York") +
  scale_color_viridis_c(name = "altura (mts)", direction = -1, option = "magma") + 
  theme_dark()

A diferencia de IDW, en el mapa (y en la escala) podemos ver que el resultado de Kriging es algo mas “suavizado”. El rango de valores se encuentra entre 0 y 120 metros.

Además, podemos calcular el error:

ggplot() + 
  geom_sf(data = lnz_pred, aes(color = var1.var), size = 0.35) + 
  ggtitle("Varianza de Kriging en interpolación de puntos de elevación", subtitle = "Altura de edificios en la ciudad de Nueva York") +
  scale_color_viridis_c(name = "varianza", direction = -1, option = "magma") + 
  theme_dark()

Como podemos ver, la varianza es mayor en las zonas en las que contamos con menos puntos en el dataset original y es menor cuando contamos con más.

Conclusiones

Al trabajar con dos interpoladores exactos (Voronoi e IDW) es interesante poder medir el error de nuestra interpolación y en eso Kriging lleva ventaja.

La desventaja de los Polígonos de Voronoi es que no hay estimación de error y que lo que queda contenido por un polígono toma de manera arbitraria el valor del punto central. Tampoco se ponderan los vecinos.

Hay que destacar que el método de Voronoi es el que mayor fidelidad tiene respecto a los datos originales, ya que los polígonos toman valores de nuestro dataset. La desventaja en este sentido es que debemos tener cuidado con valores los extremos y los outliers.

Kriging nos permitió además de interpolar, conocer la varianza. La desventaja en este caso particular fue que la superficie generada es excesivamente suave. Además, al contar con diferencias entre alturas que producen quiebres, se perdió fidelidad de la distribución espacial del proceso.

La Distancia Inversa Ponderada tiene como ventaja que funciona bien cuando el conjunto de puntos muestrales es denso ya que captura la extensión de la variación de la superficie. Si bien no hay estimación de error es la que mejor funcionó en este caso.