Distancia Punto a Linea

Paso 1. Descargar la información.

Primero había que descargar la información de la frontera norte de México a mi computadora. La primera opción que pensé hacer fué filtrar las líneas superiores de un shape de los estados del país, pero si lo piensas un poco… ¿como podría decirle a R cual es la línea superior de un estado? ¿Seleccionando los puntos más extremos en el eje Y?

Pues tal vez esa podría ser una opción, pero el algoritmo de Googlear la base de la línea de la frontera y descargar la base resultó más efectivo.

Si eres de los que les dá flojera ir a la página, descargar archivos, ir a la carpeta de Descargas, y mover el archivo a la carpeta de trabajo, puedes copipegar el código que viene a continuación (obviamente, descargando primero el paquete curl del CRAN en tu sesión de R).

# Descarga datos automaticamente, si te da flojera darle cl
curl::curl_download("https://opendata.arcgis.com/datasets/e735940321bd4383bab528a91bf526f8_0.zip?outSR=%7B%22latestWkid%22%3A3857%2C%22wkid%22%3A102100%7D",
destfile = "frontera.zip")

Esto me descargó un archivo *.zip, que yo nombré como frontera.zip, en mi directorio de trabajo. Como ya caducó mi licencia de WinRar 😔 y no puedo descomprimirlo en la computadora, utilizo el siguiente código para extraer la información:

# Deszipeado
unzip(zipfile = "frontera.zip")

Una vez descomprimido, tengo ahora un montón de archivos en mi directorio de trabajo llamados Mexico_and_US_Border con distintas extensiones, tales como *.shp, *.shx o *.dbf. Estos archivos son una de las maneras en que se guardan las bases de datos geográficas (formato shapefile), y estos archivos se pueden abrir desde R con la librería sf.

Fig. 2 Archivos que conforman en shapefile de la frontera.

Paso 2. Abriendo los datos.

Para abrir la base de datos de la frontera norte, usamos la función sf::st_read() como se ve a continuación:

# Librerias
library(sf)
library(tidyverse)
# Apertura del archivo (solo abrimos el archivo *.shp)
f <- st_read("Mexico_and_US_Border.shp")

Y ya. así lo abrimos. Para corroborar que sea efectivamente un archivo geográfico, lo ploteamos (dibujamos):

Nota. Utiliza siempre el argumento max.plot=1 para explorar tus datos (si no, tu compu lo va a pasar muy mal).

# Ploteado (te recomiendo que no le muevas al max.plot = 1)
plot(f, max.plot = 1)

Fig. 3 Polilínea de la frontera norte

Efectivamente, tiene los datos que necesitamos.

3. Definiendo el punto.

Ahora hay que determinar el punto al cual le vamos a sacar la distancia respecto de la frontera norte: como mi amigo es estudiante del Colmex, vamos a calcular la distancia de El Colegio de México a un punto (el más cercano) a la frontera norte.

Para sacar las coordenadas del Colmex, lo buscamos en Google Maps, y le picamos encima de la biblioteca para sacar sus coordenadas:

Fig. 4 Coordenadas del Colmex en Google Maps.

Ahora, creamos una tabla que almacene estas coordenadas latitud/longitud:

# Las coordenadas del Colmex en formato tabla-dataframe
pto <- data.frame(x = -99.20775, y = 19.303741) %>% 
  st_as_sf(coords = c("x", "y")) 

# Homologamos el Sistema de Coordenadas de Referencia con la base de la línea de la Frontera Norte
st_crs(pto) <- st_crs(f)

Nota. Cuando se va a trabajar con la interacción de dos bases de datos geográficas, o cuando se van a realizar mapas en leaflet, es siempre necesario homologar el Sistema de Coordenadas de Referencia, escogiendo, en este caso, el 4326 (WGS84).

4. Calculamos la distancia.

Ahora calculamos la distancia. La librería/biblioteca/paquetería sf tiene ya pre-programada una función llamada st_distance() para calcular la distancia geodésica que hay entre dos elementos geográficos del mapa. Utilizaremos esta función para calcular la distancia de un punto a la frontera.

st_distance(pto, f) 

## Units: [m]
##          [,1]
## [1,] 746695.6

Y ya, la distancia de El Colmex a el punto más cercano de la frontera norte (quién sabe cual) es de 746.6889 km ( 746,688.9 m).

Y ahí está la respuesta a la distancia, y al reto de mi amigo. Fin.

5. ¿A qué punto se calculó la distancia?

Pero ahora pueden surgir preguntas muy lógicas sobre esta respuesta… Una de ellas es: ¿A qué punto de la frontera se le calculó la distancia?

Entonces, hagamos el cálculo más detenidamente. Para esto, propongo descomponer la línea de la frontera norte en sus vértices o componentes:

# Descomponemos la linea en sus coordenadas
ptos_linea <- st_coordinates(f) %>%
  as.data.frame() %>%
  st_as_sf(coords = c("X", "Y")) 

st_crs(ptos_linea) <- st_crs(f)

class(ptos_linea)

## [1] "sf"         "data.frame"

st_crs(ptos_linea)

## Coordinate Reference System:
##   EPSG: 4326 
##   proj4string: "+proj=longlat +datum=WGS84 +no_defs"

Y ahora calculamos las distancias del Colmex a todos los vértices de la frontera norte:

# Calculo de las distancias
distancias <- st_distance(ptos_linea, pto)
# Observamos los primeros registros
head(distancias)

## Units: [m]
## [1] 775122.1 766724.1 766723.6 766720.5 766717.8 766712.2

Y ahora, sacamos la línea con la menor longitud:

# Obtencion de la distancia minima
distancia_minima <- min(distancias) 

# En kilometros
(distancia_minima / 1000) %>% as.numeric()

## [1] 746.6956

Como podemos ver, la distancia mínima es la misma distancia que calcula la función st_distance(). Ahora, ¿cuál punto de todos esos es el punto más cercano al Colmex? Para conocer esto, filtramos las distancias hasta quedarnos aquella observación que tenga la distancia igual a la distancia mínima.

# Punto minimo
punto_frontera <- 
  ptos_linea[distancias == distancia_minima,]

Y para visualizarlo, creamos una línea que va del Colmex a dicho punto:

# Linea de distancia mínima
linea <- st_linestring(matrix(c(pto[,"geometry"] %>% 
                              st_coordinates(), 
                              punto_frontera[,"geometry"] %>% st_coordinates()), ncol = 2, byrow = TRUE))

Y lo plasmamos en un mapa de leaflet:

library(leaflet)

# Hacemos el mapa 
leaflet() %>% 
  addTiles() %>% 
  addCircleMarkers(data = pto) %>% 
  addCircleMarkers(data = punto_frontera) %>% 
  addPolylines(data = linea) %>% 
  addPolylines(data = f, color = "red")