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.888889

Capitales 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       0

Para 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 9293710

El 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 9293710

A 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)] <- NA

Cuando el resultado es de clase units 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()

Buy Me A Coffee

Dr. Dominic Royé
Dr. Dominic Royé
Investigador postdoctoral senior

Relacionado