Distancias geográficas
El primer post del año 2020, lo dedicaré a una consulta que me hicieron recientemente. Me plantearon la pregunta de cómo se podría calcular la distancia más corta entre diferentes puntos y cómo saber cúal es el punto más próximo a uno dado. Cuando trabajamos con datos espaciales en R, en la actualidad lo más fácil es usar el paquete sf en combinación con la colección de paquetes tidyverse. Además usamos el paquete units que es muy útil para trabajar con unidades de medida.
Paquetes
| Paquete | Descripción |
|---|---|
| tidyverse | Conjunto de paquetes (visualización y manipulación de datos): ggplot2, dplyr, purrr,etc. |
| sf | Simple Feature: importar, exportar y manipular datos vectoriales |
| units | Proporciona unidades de medida para vectores R: conversión, derivación, simplificación |
| maps | Mapas geográficos y conjuntos de datos |
| rnaturalearth | Mapas vectoriales del mundo ‘Natural Earth’ |
# instalamos los paquetes necesarios
if(!require("tidyverse")) install.packages("tidyverse")
if(!require("units")) install.packages("units")
if(!require("sf")) install.packages("sf")
if(!require("maps")) install.packages("maps")
if(!require("rnaturalearth")) install.packages("rnaturalearth")
# cargamos los paquetes
library(maps)
library(sf)
library(tidyverse)
library(units)
library(rnaturalearth)Unidades de medida
El uso de vectores y matrices de clase units nos permite calcular y transformar unidades de medida.
# longitud
l <- set_units(1:10, m)
l## Units: [m]
## [1] 1 2 3 4 5 6 7 8 9 10# convertir a otras unidades
set_units(l, cm)## Units: [cm]
## [1] 100 200 300 400 500 600 700 800 900 1000# sumar diferentes unidades
set_units(l, cm) + l## Units: [cm]
## [1] 200 400 600 800 1000 1200 1400 1600 1800 2000# area
a <- set_units(355, ha)
set_units(a, km2)## 3.55 [km2]# velocidad
vel <- set_units(seq(20, 50, 10), km/h)
set_units(vel, m/s)## Units: [m/s]
## [1] 5.555556 8.333333 11.111111 13.888889Capitales del mundo
Vamos a usar las capitales de todo el mundo con el objetivo de calcular la distancia a la capital más próxima y indicar el nombre de la ciudad.
# conjunto de ciudades del mundo con coordenadas
head(world.cities) # proviene del paquete maps## name country.etc pop lat long capital
## 1 'Abasan al-Jadidah Palestine 5629 31.31 34.34 0
## 2 'Abasan al-Kabirah Palestine 18999 31.32 34.35 0
## 3 'Abdul Hakim Pakistan 47788 30.55 72.11 0
## 4 'Abdullah-as-Salam Kuwait 21817 29.36 47.98 0
## 5 'Abud Palestine 2456 32.03 35.07 0
## 6 'Abwein Palestine 3434 32.03 35.20 0Para convertir puntos con longitud y latitud en un objeto espacial de clase sf, empleamos la función st_as_sf(), indicando las columnas de las coordenadas y el sistema de referencia de coordenadas (WSG84, epsg:4326).
# convertimos los puntos en un objeto sf con CRS WSG84
cities <- st_as_sf(world.cities, coords = c("long", "lat"), crs = 4326)
cities## Simple feature collection with 43645 features and 4 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -178.8 ymin: -54.79 xmax: 179.81 ymax: 78.93
## Geodetic CRS: WGS 84
## First 10 features:
## name country.etc pop capital geometry
## 1 'Abasan al-Jadidah Palestine 5629 0 POINT (34.34 31.31)
## 2 'Abasan al-Kabirah Palestine 18999 0 POINT (34.35 31.32)
## 3 'Abdul Hakim Pakistan 47788 0 POINT (72.11 30.55)
## 4 'Abdullah-as-Salam Kuwait 21817 0 POINT (47.98 29.36)
## 5 'Abud Palestine 2456 0 POINT (35.07 32.03)
## 6 'Abwein Palestine 3434 0 POINT (35.2 32.03)
## 7 'Adadlay Somalia 9198 0 POINT (44.65 9.77)
## 8 'Adale Somalia 5492 0 POINT (46.3 2.75)
## 9 'Afak Iraq 22706 0 POINT (45.26 32.07)
## 10 'Afif Saudi Arabia 41731 0 POINT (42.93 23.92)En el próximo paso, simplemente filtramos por las capitales codificadas en la columna capital con 1. La ventaja del paquete sf es la posibilidad de aplicar funciones de la colección tidyverse para manipular los atributos. Además, añadimos una columna con nuevas etiquetas usando la función str_c() del paquete stringr, la cúal es similar a la de R Base paste().
# filtramos por las capitales
capitals <- filter(cities, capital == 1)
# creamos una nueva etiqueta combinando nombre y país
capitals <- mutate(capitals, city_country = str_c(name, " (", country.etc, ")"))
capitals ## Simple feature collection with 230 features and 5 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -176.13 ymin: -51.7 xmax: 179.2 ymax: 78.21
## Geodetic CRS: WGS 84
## First 10 features:
## name country.etc pop capital geometry
## 1 'Amman Jordan 1303197 1 POINT (35.93 31.95)
## 2 Abu Dhabi United Arab Emirates 619316 1 POINT (54.37 24.48)
## 3 Abuja Nigeria 178462 1 POINT (7.17 9.18)
## 4 Accra Ghana 2029143 1 POINT (-0.2 5.56)
## 5 Adamstown Pitcairn 51 1 POINT (-130.1 -25.05)
## 6 Addis Abeba Ethiopia 2823167 1 POINT (38.74 9.03)
## 7 Agana Guam 1041 1 POINT (144.75 13.47)
## 8 Algiers Algeria 2029936 1 POINT (3.04 36.77)
## 9 Alofi Niue 627 1 POINT (-169.92 -19.05)
## 10 Amsterdam Netherlands 744159 1 POINT (4.89 52.37)
## city_country
## 1 'Amman (Jordan)
## 2 Abu Dhabi (United Arab Emirates)
## 3 Abuja (Nigeria)
## 4 Accra (Ghana)
## 5 Adamstown (Pitcairn)
## 6 Addis Abeba (Ethiopia)
## 7 Agana (Guam)
## 8 Algiers (Algeria)
## 9 Alofi (Niue)
## 10 Amsterdam (Netherlands)Calcular distancias
La distancia geográfica (euclidiana o de gran círculo) se calcula con la función st_distance(), o bien entre dos puntos, entre un punto y otros múltiples o entre todos. En el último caso obtenemos una matriz simétrica de distancias (NxN), tomados por pares de un conjunto. En la diagonal encontramos las combinaciones entre los mismos puntos dando todas nulas.
| A | B | C | |
| A | 0 | 310 | 323 |
| B | 310 | 0 | 133 |
| C | 323 | 133 | 0 |
Cuando queremos saber, por ejemplo, la distancia de Amsterdam a Abu Dhabi, Washington y Tokyo pasamos dos objetos espaciales.
# calcular la distancia
dist_amsterdam <- st_distance(slice(capitals, 10),
slice(capitals, c(2, 220, 205)))
dist_amsterdam # distancia en metros## Units: [m]
## [,1] [,2] [,3]
## [1,] 5163124 6187634 9293710El resultado es una matriz de una fila o de una columna (en función del orden de los objetos) con clase de units. Así es posible cambiar fácilmente a otra unidad de medida. Si queremos obtener un vector sin clase units, únicamente aplicamos la función as.vector().
# cambiamos de m a km
set_units(dist_amsterdam, "km")## Units: [km]
## [,1] [,2] [,3]
## [1,] 5163.124 6187.634 9293.71# units class a vector
as.vector(dist_amsterdam)## [1] 5163124 6187634 9293710A continuación estimamos la matriz de distancia entre todas las capitales. Es importante convertir los valores nulos a NA para obtener posteriormente el índice correcto de la matriz.
# calcular la distancia
m_distance <- st_distance(capitals)
# matriz
dim(m_distance)## [1] 230 230# cambiamos de m a km
m_distance_km <- set_units(m_distance, km)
# reemplazamos la distance de 0 con NA
m_distance_km[m_distance_km == set_units(0, km)] <- NAunits es necesario usar la misma clase para poder hacer consultas logicas. Por ejemplo, set_units(1, m) == set_units(1, m) vs. set_units(1, m) == 1.
Con el objetivo de obtener la distancia más corta, además de la posición de la misma, usamos la función apply() que a su vez nos permite aplicar la función which.min() y min() sobre cada fila. También sería posible emplear la función sobre columnas que daría el mismo resulado. Para finalizar, añadimos los resultados como nuevas columnas con la función mutate(). Las posiciones en pos nos permiten obtener los nombres de las ciudades más próximas.
# obtenemos la posición de la ciudad y la distancia
pos <- apply(m_distance_km, 1, which.min)
dist <- apply(m_distance_km, 1, min, na.rm = TRUE)
# añadimos la distancia y obtenemos el nombre de la ciudad
capitals <- mutate(capitals, nearest_city = city_country[pos],
geometry_nearest = geometry[pos],
distance_city = dist)Mapa de distancias a la próxima capital
Por último, construimos un mapa representando la distancia en circulos proporcionales. Para ello, usamos la gramática habitual de ggplot() añadiendo la geometría geom_sf(), primero para el mapamundi de fondo y después para los circulos de las ciudades. En aes() indicamos con el argumento size = distance_city la variable que debe ser mapeado proporcionalmente. La función theme_void() elimina todos los elementos de estilo. Además, definimos con la función coord_sf() una nueva proyección indicando el formato proj4.
# mapamundi
world <- ne_countries(scale = 10, returnclass = "sf")
# mapa
ggplot(world) +
geom_sf(fill = "black", colour = "white") +
geom_sf(data = capitals,
aes(size = distance_city),
alpha = 0.7,
fill = "#bd0026",
shape = 21,
show.legend = 'point') +
coord_sf(crs = "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs") +
labs(size = "Distance (km)", title = "Distance to the next capital") +
theme_void()
