library(tidyverse)
library(tidygeocoder)
library(sf)11 Geocoding
11.1 Introducción
La geocodificación (geocoding) es un proceso muy utilizado en el ámbito de la ciencia de datos geoespaciales, mediante el cual se convierten direcciones postales o nombres de lugares en coordenadas geográficas (latitud y longitud). Este proceso permite mapear datos físicos en un contexto geográfico, facilitando una amplia gama de aplicaciones, desde el análisis de patrones de distribución hasta la optimización de rutas y la planificación urbana. La geocodificación no solo es fundamental para visualizar datos sobre mapas, sino también para realizar análisis espaciales que pueden revelar correlaciones y tendencias ocultas basadas en la ubicación. En esta práctica, exploraremos en profundidad la librería tidygeocoder (Cambon et al., 2021), una herramienta poderosa y flexible diseñada para integrarse perfectamente con el ecosistema tidyverse (Wickham, 2022).
tidygeocoder simplifica notablemente el proceso de obtener coordenadas a partir de direcciones postales, permitiéndonos acceder a múltiples servicios de geocodificación a través de una única interfaz. En esta sesión aprenderemos a utilizar esta librería para transformar direcciones postales en datos espaciales precisos (coordenadas), explorando sus capacidades para realizar tanto geocodificación hacia adelante (forward geocoding) como geocodificación inversa (reverse geocoding), es decir, obtener direcciones postales a partir de coordenadas.
11.2 Callejero de Madrid
En primer lugar, obtendremos un conjunto de datos que podamos utilizar para geocodificar. El callejero de Madrid, disponible en el Portal de datos abiertos del Ayuntamiento de Madrid (también aquí), contiene una gran cantidad de direcciones postales, por lo que nos será de gran utilidad para ver el funcionamiento del paquete tidygeocoder. Pero además, este conjunto de datos también incluye las coordenadas de cada dirección postal, lo que nos facilitará la comprobación de los resultados obtenidos.
Cargamos librerías:
Cargamos y analizamos el dataset (descripción del fichero aquí):
callejero_madrid <- readr::read_delim(file = "data/DireccionesVigentes_20250218.csv", delim = ";",
locale = readr::locale(encoding = "ISO-8859-2", decimal_mark = ","))
dim(callejero_madrid)[1] 213318 20dplyr::glimpse(callejero_madrid)Rows: 213,318
Columns: 20
$ COD_VIA <dbl> 31001337, 31001337, 31001337, 31001337,…
$ VIA_CLASE <chr> "AUTOVÍA", "AUTOVÍA", "AUTOVÍA", "AUTOV…
$ VIA_PAR <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
$ VIA_NOMBRE <chr> "A-1", "A-1", "A-1", "A-1", "A-1", "A-1…
$ VIA_NOMBRE_ACENTOS <chr> "A-1", "A-1", "A-1", "A-1", "A-1", "A-1…
$ CLASE_APP <chr> "KILÓMETRO", "KILÓMETRO", "KILÓMETRO", …
$ NUMERO <dbl> 10000, 10000, 11000, 11000, 12000, 1200…
$ CALIFICADOR <chr> " EN", " SA", " EN", " SA", " EN", " SA…
$ TIPO_NDP <chr> "PARCELA", "PARCELA", "PARCELA", "PARCE…
$ COD_NDP <dbl> 31031089, 31031088, 31031091, 31031090,…
$ DISTRITO <dbl> 8, 16, 8, 16, 8, 16, 8, 16, 15, 15, 16,…
$ BARRIO <dbl> 6, 6, 6, 6, 6, 6, 6, 6, 6, 5, 2, 8, 2, …
$ COD_POSTAL <dbl> 28050, 28050, 28050, 28050, 28050, 2805…
$ UTMX_ED <dbl> 443056.3, 443122.5, 443675.2, 443737.5,…
$ UTMY_ED <dbl> 4482503, 4482491, 4483309, 4483252, 448…
$ UTMX_ETRS <dbl> 442947.0, 443013.1, 443565.8, 443628.2,…
$ UTMY_ETRS <dbl> 4482296, 4482283, 4483102, 4483044, 448…
$ LATITUD <chr> "40°29'21.84'' N", "40°29'21.45'' N", "…
$ LONGITUD <chr> "3°40'23.56'' W", "3°40'20.75'' W", "3°…
$ ANGULO_ROTULACION <chr> "78.87", "77.25", "41.35", "47.56", "67…
Ejercicio
- Crea un objeto sf para cada formato de coordenadas.
- ¿Qué código EPSG has utilizado para cada caso?
- Comprueba, por pares, que las coordenadas geográficas/proyectadas son equivalentes.
11.3 Forward Geocoding
El paquete tidygeocoder tiene 2 funciones que permiten realizar forward geocoding: geo() y geocode(). La principal diferencia entre ambas funciones es que la primera trabaja con un vector y la segunda con un data.frame (tibble).
Antes de ver con detalle cada función, vamos a crear en nuestro dataset una nueva variable la dirección completa:
callejero_madrid <- callejero_madrid %>%
mutate(DIRECCION = paste0(VIA_CLASE," ", VIA_NOMBRE, ", ", NUMERO, ", MADRID"))
set.seed(123)
(ids_muestra <- sample(1:nrow(callejero_madrid), 5))[1] 182735 188942 134058 124022 160997(muestra <- callejero_madrid[ids_muestra,])11.3.1 tidygeocoder::geo()
crds1 <- geo(address = muestra$DIRECCION)Passing 5 addresses to the Nominatim single address geocoderQuery completed in: 5.1 secondscrds1crds2 <- geo(street = paste(muestra$VIA_NOMBRE,muestra$NUMERO, sep = ", "),
city = rep("Madrid",nrow(muestra)),
country = rep("Spain",nrow(muestra)),
postalcode = muestra$COD_POSTAL)Passing 5 addresses to the Nominatim single address geocoder
Query completed in: 5.1 secondscrds2Comprobemos las diferencias en base a las distancias entre puntos:
crds1_sf <- sf::st_as_sf(crds1, coords=c("lat","long"), crs=4258, na.fail=F)
crds2_sf <- sf::st_as_sf(crds2, coords=c("lat","long"), crs=4258, na.fail=F)
sf::st_distance(crds1_sf,crds2_sf,by_element = T)Units: [m]
[1] 0.000 0.000 0.000 0.000 3427.641También podemos obtener resultados más completos, que variarán en función del proveedor de servicios utilizado (por defecto method = "osm"):
(crds3 <- geo(address = muestra$DIRECCION, full_results = T))Passing 5 addresses to the Nominatim single address geocoderQuery completed in: 5.1 seconds11.3.2 tidygeocoder::geocode()
Hagamos lo mismo que en el apartado anterior, pero esta vez con la función geocode():
crds1b <- tidygeocoder::geocode(muestra, address = DIRECCION)Passing 5 addresses to the Nominatim single address geocoderQuery completed in: 5.1 secondscrds1bcrds2b <- tidygeocoder::geocode(muestra,
street = VIA_NOMBRE,
postalcode = COD_POSTAL)Passing 5 addresses to the Nominatim single address geocoder
Query completed in: 5.1 secondscrds2bcrds1b_sf <- sf::st_as_sf(crds1b, coords=c("lat","long"), crs=4258, na.fail=F)
crds2b_sf <- sf::st_as_sf(crds2b, coords=c("lat","long"), crs=4258, na.fail=F)
sf::st_distance(crds1_sf,crds1b_sf,by_element = T)Units: [m]
[1] 0 0 0 0 011.4 Métodos
En los ejemplos anteriores hemos utilizado el proveedor de servicios por defecto method = "osm", aunque realmente el servicio de geocodificación se llama Nominatim, un servicio de geocodificación gratuito basado en datos de OpenStreetMap (OSM).
Como hemos podido comprobar, en ocasiones, la precisión de OSM no es todo lo buena que desearíamos. Para tratar de resolver esta restricción, tidygeocoder permite trabajar con otros métodos (servicios) como arcgis, google, opencage, mapbox, here, tomtom o mapquest, entre otros (Pérez y Aybar, 2024). Estos métodos están basados en proveedores de servicios privados, por lo que en la mayoría de los casos necesitamos una clave (API key) para poder usarlos. Antes de continuar, veamos algunos ejemplos para comprobar esto:
crds_arcgis <- tidygeocoder::geocode(muestra, address = DIRECCION, method = "arcgis")Passing 5 addresses to the ArcGIS single address geocoderQuery completed in: 2 secondscrds_google <- tidygeocoder::geocode(muestra, address = DIRECCION, method = "google")Passing 5 addresses to the Google single address geocoderQuery completed in: 0.4 secondscrds_bing <- tidygeocoder::geocode(muestra, address = DIRECCION, method = "bing")Passing 5 addresses to the Bing single address geocoderQuery completed in: 1 seconds
Atención
Para poder utilizar las API key obtenidas de los diferentes proveedores, hay que añadirlas en fichero .Renviron. Una forma sencilla de acceder y modificar a este fichero es mediante la sentencia usethis::edit_r_environ().
A continuación, se ofrecen los enlaces a las diferentes plataformas:
11.5 Reverse geocoding
Ya hemos visto cómo podemos obtener coordenadas a partir de direcciones postales (forward geocoding). Veamos ahora el proceso inverso. Las funciones que usaremos en este caso reverse_geo() y reverse_geocode().
11.5.1 tidygeocoder::reverse_geo()
coords <- mad3[ids_muestra,]
dirs1 <- tidygeocoder::reverse_geo(lat = coords$lat, long = coords$long)Passing 5 coordinates to the Nominatim single coordinate geocoderQuery completed in: 5.1 seconds11.5.2 tidygeocoder::reverse_geocode
dirs2 <- coords %>%
tidygeocoder::reverse_geocode(lat = lat, long = long)Passing 5 coordinates to the Nominatim single coordinate geocoderQuery completed in: 5.1 secondsdirs3 <- coords %>%
tidygeocoder::reverse_geocode(lat = lat, long = long, full_results=T)Passing 5 coordinates to the Nominatim single coordinate geocoder
Query completed in: 5.1 seconds11.6 Alternativas
Si bien es cierto que tidygeocoder es el paquete más versátil de R para hacer geocoding (Pérez y Aybar, 2024), existen otros paquetes alternativos que ofrecen ciertas ventajas. Algunos son: ggmap (Kahle y Wickham, 2013), nominatimlite (Hernangómez, 2024), tmaptools (Tennekes, 2021), hereR (Unterfinger, 2023), opencage (Possenriede et al., 2021), mapboxapi (Walker, 2023), mapquestr (Chiou, 2024).
library(ggmap)
library(nominatimlite)
library(tmaptools)11.6.1 ggmap
Este paquete de R permite la visualización espacial combinando información de mapas estáticos de Google Maps, OpenStreetMap, Stamen Maps o CloudMade Maps con la implementación de la gramática de gráficos de ggplot2. Además, introduce varias funciones útiles que permiten acceder a las API de Geocodificación, Matriz de Distancia y Direcciones de Google, ofreciendo un marco fácil, consistente y modular para gráficos espaciales y herramientas convenientes para el análisis de datos espaciales.
crds_ggmap_1 <- ggmap::geocode(location = muestra$DIRECCION)
crds_ggmap_2 <- ggmap::mutate_geocode(muestra,
location = DIRECCION,
output = "more",
source = "google")
mapa_base <- ggmap::get_map(location = c(lon = mean(crds_ggmap_1$lon),
lat = mean(crds_ggmap_1$lat)))
ggmap(mapa_base) +
geom_point(data = crds_ggmap_1,
aes(x = lon, y = lat),
color = "red", size = 3)
11.6.2 nominatimlite
El objetivo de nominatimlite es proporcionar una interfaz ligera para geocodificar direcciones, basada en la API de Nominatim. También permite cargar objetos espaciales usando el paquete sf.
crds_nominatim_1 <- nominatimlite::geo_lite(address = muestra$DIRECCION)
|
| | 0%
|
|========== | 20%
|
|==================== | 40%
|
|============================== | 60%
|
|======================================== | 80%
|
|==================================================| 100%crds_nominatim_2 <- nominatimlite::geo_lite_sf(address = muestra$DIRECCION)
|
| | 0%
|
|========== | 20%
|
|==================== | 40%
|
|============================== | 60%
|
|======================================== | 80%
|
|==================================================| 100%class(crds_nominatim_2)[1] "sf" "tbl_df" "tbl" "data.frame"dirs_nominatim_1 <- nominatimlite::reverse_geo_lite(lat = coords$lat, long = coords$long)
|
| | 0%
|
|========== | 20%
|
|==================== | 40%
|
|============================== | 60%
|
|======================================== | 80%
|
|==================================================| 100%dirs_nominatim_2 <- nominatimlite::reverse_geo_lite_sf(lat = coords$lat, long = coords$long)
|
| | 0%
|
|========== | 20%
|
|==================== | 40%
|
|============================== | 60%
|
|======================================== | 80%
|
|==================================================| 100%class(dirs_nominatim_2)[1] "sf" "tbl_df" "tbl" "data.frame"dirs_nominatim_2 %>%
leaflet::leaflet() %>%
leaflet::addTiles() %>%
leaflet::addMarkers(popup = ~address)También facilita la extracción de información de OpenStreetMaps:
GVA <- nominatimlite::geo_lite_sf("Comunidad Valenciana", points_only = FALSE)
class(GVA)[1] "sf" "tbl_df" "tbl" "data.frame"mcdonalds <- nominatimlite::geo_lite_sf("McDonald's, Comunidad Valenciana",
limit = 50,
custom_query = list(countrycodes = "es"))
GVA %>%
ggplot() +
geom_sf(data = GVA) +
geom_sf(data = mcdonalds, col = "red")
11.6.3 tmaptools
Este paquete proporciona un conjunto de herramientas que complementan el paquete tmap[R-tmap], facilitando la lectura, manipulación y procesamiento de datos espaciales. Incluye funciones para geocodificación, cálculo de áreas y distancias, y preparación de datos para la creación de mapas temáticos.
crds_tmaptools_1 <- tmaptools::geocode_OSM(muestra$DIRECCION, projection = 25830)
cbind(muestra$UTMX_ETRS, crds_tmaptools_1$x) [,1] [,2]
[1,] 443860.8 443872.7
[2,] 444925.1 444678.7
[3,] 434028.4 434009.5
[4,] 445717.2 445712.6
[5,] 441641.7 441605.8dirs_tmaptools_1 <- tmaptools::rev_geocode_OSM(x=coords$long, y=coords$lat)
dirs_tmaptools_2 <- tmaptools::rev_geocode_OSM(mad3_sf_4326[ids_muestra,])11.7 Caso de uso: Guardia Civil
La Guardia Civil es uno de los cuerpos de seguridad más antiguos de España, con una historia que se remonta a 1844. Se trata de una institución de naturaleza militar que forma parte de las Fuerzas y Cuerpos de Seguridad del Estado, cuya principal misión es proteger el libre ejercicio de los derechos y libertades de los ciudadanos y garantizar la seguridad ciudadana en todo el territorio español. Su ámbito de actuación se extiende por toda España, incluyendo áreas urbanas y rurales, aunque tradicionalmente ha tenido una presencia más marcada en las zonas rurales y en las pequeñas localidades.
La Guardia Civil desempeña un amplio abanico de funciones que abarcan desde la vigilancia de carreteras y del tráfico, la protección de la naturaleza, la lucha contra el narcotráfico y el terrorismo, hasta la seguridad ciudadana y la investigación criminal. Su estructura organizativa se divide en varias unidades especializadas para atender las distintas áreas de seguridad y protección.
En su compromiso con la comunidad, la Guardia Civil cuenta con dependencias que ofrecen un área de atención al ciudadano abierta las 24 horas del día. Estas instalaciones están diseñadas para brindar apoyo y asistencia a los ciudadanos en cualquier momento, reforzando así la cercanía y la accesibilidad de este cuerpo de seguridad a la población a la que sirve.
En la página web oficial de la Guardia Civil se ofrecen varios datasets con las direcciones postales de determinadas Dependencias. El dataset también está disponible aquí.
Ejercicio
- Descarga el dataset que contiene la direcciones postales de los puestos “24 horas” de la Guardia Civil.
- Selecciona una muestra aleatoria de 15 observaciones (fija la semilla en 1234).
- Obtén las coordenadas usando el servicio de Nominatim.
- Obtén las coordenadas usando el servicio de ArcGIS.
- ¿Qué cantidad de valores faltantes (NA) has obtenido para cada dataset?
- ¿Cuánto tiempo has tardado en obtener las coordenadas para cada dataset?
- ¿Qué servicio consideras que ha reportado mejores resultados? ¿Por qué?