Los diagramas de Voronoi, también conocidos como Polígonos de Thiessen son representaciones de la influencia espacial lineal de entidades. La influencia se calcula en función de la distancia y se trazan semiplanos a partir de las mediatrices de la Triangulación de Delaunay. Una técnica de interpolación espacial muy simple y muy antigua que sigue siendo útil para la exploración de datos espaciales.
Se construye a partir de puntos donde cada punto tendrá influencia hasta la mitad de la distancia euclideana medida hasta sus vecinos más próximos. Tomando como base la definición formal de distancia euclideana entre los puntos \(a\) y \(b\) como:
\[ D_{a - b}=\sqrt{(a_x - b_x)^2 + (a_y - b_y)^2} \]
…podemos dividir el espacio en semiplanos usando como referencia el punto ubicado a \(0.5D_{a-b}\) de cada par de puntos. La intersección de todos los \(n\) semiplanos alrededor del punto de interés será la región Voronoi del ese punto. El proceso iterativo sobre todos los \(N\) puntos del dataset permitirá generar todas las regiones que conformarán el Diagrama de Voronoi.
Por ejemplo, para la región de Voronoi del punto gris en la imagen (recuadro grande) resulta de la intersección de todos los semiplanos de color gris en los recuadros más pequeños
Intersección de semiplanos para construcción de una región de voronoi
Datos
Este mismo procedimiento lo aplicaremos a unos datos de ejemplo obtenidos a partir del WFS del Instituto Geográfico Militar de Ecuador. Primero vamos a visualizar los puntos disponibles de la siguiente forma:
#···· Para descargar desde el WFS del IGM ....
# library(httr)
# library(ows4R)
# wfs_igm_ec <- "http://www.geoportaligm.gob.ec/p_geoinformacion/wms"
# url <- parse_url(wfs_igm_ec)
# url$query <- list(service = "WFS",
# version = "1.0.0",
# request = "GetFeature",
# typename = "Infraestructura_Servicios:IS_EstablecimientoSalud",
# bbox = "695435.4, 9553697, 702943.2, 9560856",
# outputFormat = "application/json")
# request <- build_url(url)
# salud_igm <- st_read(request) %>%
# mutate(institucion = if_else(str_detect(descripcio, "IESS"),
# "SEGURO SOCIAL", "PÚBLICO"))
# Si ya los tienes descargados se pueden leer localmente
salud_igm <- st_read(dsn = "data/geodata.gpkg",
layer = "voronoi_data")
mapview(salud_igm, zcol = "institucion")Cálculo de las los polígonos de Voronoi con el paquete deldir
El paquete deldir es una de las mejores opciones para hacer la triangulación y teselación. Requiere las coordenadas de los puntos en una matriz de x e y. Para lo cual necesitamos extraer la matriz de coordenadas de los puntos de ejemplo. También vamos a crear un vector de colores, basado en una variable. Para el ejemplo usaremos el atributo gid de los datos originales.
# extraemos la matriz de coordenadas
coordenadas <- sf::st_coordinates(salud_igm)
# creamos la paleta de colores en función del atributo `gid`
pal <- colorRampPalette(RColorBrewer::brewer.pal(9,"YlGnBu")) # función que genera colores
colores <- pal(length(unique(salud_igm$gid))) # Lista de colores
colores_variable <- recode(salud_igm$gid, # Seleccionamos los colores en función de `gid`
!!!purrr::set_names(colores,
sort(unique(salud_igm$gid))))El cálculo de las regiones de Voronoi para cada punto se hace con la función deldir::deldir. Esta función creará un objeto de lista que contiene entre otras cosas dos data frame: Uno para la Triangulación de Delaunay y otro para los segmentos de línea que delimitan regiones de Voronoi de los puntos . Y con ayuda de la función deldir::tile.list() vamos a convertir los segmentos a polígonos.
# Generar los segmentos límite de las regiones de voronoi
voronoi_bordes <- deldir::deldir(x = coordenadas[,1],
y = coordenadas[,2],
eps = 1e-09,
sort = TRUE,
plotit = FALSE,
digits = 6,
suppressMsge = TRUE)
# Convertir los bordes a polígonos
voronoi_poligonos <- deldir::tile.list(voronoi_bordes)Veamos el resultado, las áreas en rellenas de color son la regiones de Voronoi obtenidas a partir de los puntos, y las líneas grises corresponden a la triangulación de Delaunay. Los números son los valores de gid de los datos originales.
# Hacer el plot
{
plot(voronoi_poligonos, fillcol = colores_variable, close = TRUE,
showpoints = TRUE, pch = 20, cex = 4, col.pts = "grey70")
plot(voronoi_bordes, wlines = "triang", add = TRUE,
cmpnt_col = c("grey50", "black", NA),
cmpnt_lty = 2)
text(coordenadas, labels = salud_igm$gid, col = "grey20")
}
Figure 1: Diagrama de voronoi con base y ggplot2
Calcular regiones de voronoi con dismo
Una forma corta de obtener los polígonos que corresponden a las regiones de Voronoi es mediante la función dismo::voronoi(). Esta función es un atajo para quienes trabajamos con datos espaciales. Permite generar las teselas a partir de un objeto de clase “Spatial”. El proceso a partir de nuestros datos que son de clase “sf” empieza con la coerción a “Spatial” y luego a la función dismo::voronoi(). Aprovecharemos la cadena de procesos para hacer el plot con ggplot2. Así:
# ggplot2 + dismo + Colores Brewer: Paleta "YlGnBu"
as(salud_igm, "Spatial") %>%
dismo::voronoi() %>%
st_as_sf() %>%
ggplot() +
geom_sf(aes(fill = gid)) +
geom_sf_label(aes(label = gid)) +
scale_fill_gradientn(colors = RColorBrewer::brewer.pal(9,"YlGnBu"))
Figure 2: Diagrama de voronoi con dismo y ggplot2
Una desventaja de usar dismo es que perdemos la posibilidad de conservar la triangulación de Delaunay como objeto y solo se conserva la teselación. Aunque esto no supone un gran inconveniente.
Calcular regiones de voronoi con el paquete sf
La última forma que exploraremos es calcular los polígonos usando la función sf::st_voronoi. Dependiendo del tipo de geometrías de tus datos (en nuestro caso ‘sfc_POINT’ ‘sfc’), se requerirá convertir en ‘XY’ ‘MULTIPOINT’ ‘sfg’, para lo cual haremos uso de la función base::do.call(). El cálculo requiere definir un marco límite para lo cual se usará el bbox de los datos. Hasta aquí ya habrás obtenido una colección de geometrías de polígonos que conforman las regiones de Voronoi. Solo falta procesar esta colección para transformar en un objeto de clase “sf data.frame” que conserve los datos originales. Esto se hará con las funciones de coerción del propio paquete sf y usando una unión espacial con los datos originales. El código siguiente contiene todos los pasos mencionados y el gráfico hecho con ggplot2
# ggplot2 + sf + Colores Brewer: Paleta "YlGnBu"
g <- sf::st_geometry(salud_igm)
bb <- sf::st_as_sfc(st_bbox(salud_igm))
voronoi_sf <- st_voronoi(do.call(c, g), st_geometry(bb)) %>%
st_collection_extract() %>%
st_sfc(crs = st_crs(salud_igm)) %>%
st_sf() %>%
st_join(salud_igm)
ggplot(voronoi_sf) +
geom_sf(aes(fill = gid)) +
geom_sf_label(aes(label = gid)) +
scale_fill_gradientn(colors = RColorBrewer::brewer.pal(9, "YlGnBu")) +
scale_color_gradientn(colors = RColorBrewer::brewer.pal(9, "YlGnBu")) +
coord_sf(xlim = range(st_coordinates(bb)[,1]),
ylim = range(st_coordinates(bb)[,2]))
Figure 3: Diagrama de voronoi con sfy ggplot2
Los diagramas de Voronoi son una forma muy rápida de saber la influencia espacial de los puntos monitoreados. Sus aplicaciones en el tiempo han variado desde explorarción visual hasta imputación de datos en base al vecino más próximo de una coordenada. En fin, espero que te sea útil y que puedas tener una mejor comprensión de cómo se hacen estos diagramas en R.
Si vez algún error o quieres sugerir algo, por favor déjanos un comentario!
