library(tidyverse)
library(leaflet)
library(sf)
library(units)



Siempre que queremos hacer operaciones matemáticas en base a objetos espaciales como por ejemplo calcular áreas o medir distancias entre puntos debemos constatar que esos objetos se correspondan con una proyección cartográfica. A la hora de buscar datos georeferenciados es usual encontrarlos con datos de coordenadas geográficas (latitud/longitud) y bajo el sistema de referencia WGS84 por lo que debemos convertirlos para poder hacer esos cálculos.

Sistema de referencia

Un sistema de referencia es un recurso matemático que permite asignar coordenadas a puntos localizados en la superficie terrestre. El sistema de referencia tiene dos componentes básicos:

  • el elipsoide: es la geometría que se asemeja al geoide. Se define por un semi-eje ecuatorial y un semi-eje polar.


  • el dátum: es el punto de coincidencia entre el geoide y el elipsoide y habilita el cálculo de coordenadas.

Existen varios sistemas de referencia y proyecciones cartográficas. Estos están definidos por parámetros y codificados bajo las siglas EPSG - SRID o SRC - CRS. Podemos consultarlos aquí

EPSG: son las iniciales de un ex-organismo (European Petroleum Survey Group) que desarrolló un repositorio codificado de parámetros geodésicos que contiene información sobre sistemas de referencia y proyecciones cartográficas de todo el mundo.

SRID: es un identificador de referencia espacial asociado con un sistema de coordenadas específico. Hay varios SRID estándar reconocidos como los define la EPSG.

SRC/CRS: son las siglas de “Coordinate Reference System”.

El código y los parámetros los leemos con la función st_crs().

WGS84

Usamos como ejemplo el dataset de conos volcánicos que lo podemos leer o descargar desde la IDECAT (Infraestructura de Datos Espaciales de la provincia de Catamarca).

conos <- st_read("http://200.43.169.149/geoserver/wfs?srsName=EPSG%3A4326&typename=geonode%3Acono&outputFormat=json&version=1.0.0&service=WFS&request=GetFeature")
## Reading layer `OGRGeoJSON' from data source `http://200.43.169.149/geoserver/wfs?srsName=EPSG%3A4326&typename=geonode%3Acono&outputFormat=json&version=1.0.0&service=WFS&request=GetFeature' using driver `GeoJSON'
## Simple feature collection with 76 features and 6 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -68.88783 ymin: -27.77598 xmax: -67.2974 ymax: -25.18257
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs

Al leer el dataset de los conos volcánicos nos dice el código EPSG/SRID y en “proj4string” sus parámetros. Éstos mismos datos los podemos consultar con la función st_crs():

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

Según el ejemplo, se usa el Sistema WGS84 con dátum con la misma denominación, define el código 4326 y que los datos de coordenadas están en latitud/longitud.

Sistema de coordenadas

Es un sistema de representación matemático que permite definir la posición de un punto sobre la superficie terrestre. Puede ser:

  • de coordenadas geográficas, como el caso de los conos volcánicos. Se basa en coordenadas angulares de latitud (en el sentido de los meridianos) y longitud (en el sentido de los paralelos) para determinar la posición de un punto. Lo ideal es que se expresen en grados decimales (-58.435667; -34.606513) aunque también se suelen expresar en grados sexagesimales ( 58°26’8.40“O; 34°36’23.45”S).

  • de coordenadas planas: se basa en sistemas de proyección cartográfica, es decir, transformaciones matemáticas que permiten convertir los puntos del elipsoide en un plano.

WGS84 a POSGAR

POSGAR 07:
Es el Marco de Referencia Geodésico Nacional de la República Argentina. Un marco de referencia es la materialización de un sistema de referencia a partir de la construcción, la medición y el posterior cálculo de las coordenadas de una serie de puntos o pilares localizados sobre la superficie terrestre. Dichos puntos conforman una Red Geodésica. Este marco de referencia es el punto de partida para la confección de la cartografía, desarrollo del catastro, planos de obras civiles y otras actividades.
POSGAR es compatible con la proyección Gauss Krüger que divide al país en 7 fajas de 3° de longitud:



Parámetros según faja:
Para cada faja varían los valores de, obviamente el del código EPSG y algunos parámetros definidos en proj4string: el valor de la longitud de origen (“lon_0”) que es el valor del meridiano central de la faja y el valor del falso este (“x_0”):

Faja 1 - 5343: “+proj=tmerc +lat_0=-90 +lon_0=-72 +k=1 +x_0=1500000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs”.
Faja 2 - 5344: “+proj=tmerc +lat_0=-90 +lon_0=-69 +k=1 +x_0=2500000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs”.
Faja 3 - 5345: “+proj=tmerc +lat_0=-90 +lon_0=-66 +k=1 +x_0=3500000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs”.
Faja 4 - 5346: “+proj=tmerc +lat_0=-90 +lon_0=-63 +k=1 +x_0=4500000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs”.
Faja 5 - 5347: “+proj=tmerc +lat_0=-90 +lon_0=-60 +k=1 +x_0=5500000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs”.
Faja 6 - 5347: “+proj=tmerc +lat_0=-90 +lon_0=-57 +k=1 +x_0=6500000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs”.
Faja 7 - 5347: “+proj=tmerc +lat_0=-90 +lon_0=-54 +k=1 +x_0=7500000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs”.


Para saber a qué faja corresponde transformar nuestro dataset, observamos el promedio de las longitudes en el “bbox”, o si fueran puntos, las longitudes de esos puntos. En el caso de los conos volcánicos, el promedio de las longitudes se acerca al valor de -69, es decir Faja 2.

conos_metros = st_transform(conos, 5344)
st_crs(conos_metros)
## Coordinate Reference System:
##   EPSG: 5344 
##   proj4string: "+proj=tmerc +lat_0=-90 +lon_0=-69 +k=1 +x_0=2500000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"



Cálculo de área

area_conos <- conos_metros %>%
  mutate(area = round(set_units(st_area(conos_metros), km^2), 1))

area_conos$area
## Units: [km^2]
##  [1]   8.9  41.8  14.0  14.0   3.1  18.0   0.4 201.4  44.8  69.4 193.6
## [12]   4.5  11.4   2.4  21.1   1.0   3.1   7.7   3.5  10.7   3.2   1.4
## [23]  81.0  22.7   9.2   0.1   3.3  10.5   4.3  17.3   7.4  80.7   5.2
## [34] 142.7  75.5   0.7   0.7   1.3   0.4   0.2   2.0   0.4   0.6   4.1
## [45]   1.3   0.3   1.2   0.2   6.2  15.1  64.4  73.6  38.2  59.9  58.7
## [56]  49.9  15.8  11.5   5.3   5.5  51.5  90.6   0.1  12.3  23.4  70.0
## [67]  17.1 205.2   0.9 294.6  17.1  80.6   9.2  20.1   2.2  47.9



Consulta del área en el mapa

Consultamos el área del cono volcánico haciendo click en el mapa:



WGS84 a UTM

Cuando trabajamos con datasets que no son de Argentina, debemos conocer el sistema que aplica un determinado país o bien, lo que nos conviene es convertir los datos a UTM (Universal Transverse Mercator).

Con el sistema UTM, las coordenadas de un punto se definen por una zona UTM a partir de un grillado de 60 husos numerados, cada uno de los cuales abarca 6° de longitud y en sentido N-S, abarca 20 zonas identificadas con letras desde la C hasta la X. Las bandas C a M están en el hemisferio sur y las bandas N a X están en el hemisferio norte. Aquí una representación de las zonas UTM:



Cálculo de distancia

Tenemos los datos de coordenadas geográficas de dos terminales de ómnibus de Brasil y queremos saber la distancia lineal entre esos dos puntos:

p1 <- data.frame(
  nombre = c("Rodoviária de Belo Horizonte"),
  lat = c(-19.913838),
  long = c(-43.941907))

p2 <- data.frame(
  nombre = c("Tietê Bus Terminal de Sao Paulo"),
  lat = c(-23.51620),
  long = c(-46.62347))

Asignamos a los puntos el código 4326 del sistema WGS84

p1_sf <- st_as_sf(p1, coords = c("long", "lat"), crs = 4326)
p2_sf <- st_as_sf(p2, coords = c("long", "lat"), crs = 4326)

Proyectamos a la faja UTM 22S y calculamos la distancia. Por la proyección, la distancia la calcula en metros por lo que debemos cambiar la unidad según la escala.

distancia <- st_distance(
  st_transform(p1_sf, crs = 32723),
  st_transform(p2_sf, crs = 32723),
  by_element = TRUE)

distancia
## 485692.8 [m]
round(set_units(distancia, km), 1)
## 485.7 [km]



Gauss Krüger Buenos Aires a WGS84

Otro de los sistemas locales en los que podemos encontrar datasets, en este caso de la Ciudad Autónoma de Buenos Aires, es la proyección llamada “Gauss Krüger Buenos Aires”. Sus parámetros son:

  • DATUM: CAMPO INCHAUSPE
  • ESFEROIDE: INTERNATIONAL 1924
  • PROYECCION: TRANSVERSE_MERCATOR
  • FALSO ESTE: 100000
  • FALSO NORTE: 100000
  • MERIDIANO CENTRAL: -58.4627
  • FACTOR DE ESCALA: 0.999998
  • LATITUD DE ORIGEN: -34.6297166
  • UNIDAD: METROS

En este ejemplo usamos un shape (shp.) de veredas descargado de BA Data. Para ejemplificar, seleccionamos una parte del dataset.

Si consultamos los parámetros de la proyección nos dice “No EPSG code”:

st_crs(veredas)
## Coordinate Reference System:
##   No EPSG code
##   proj4string: "+proj=tmerc +lat_0=-34.6297166 +lon_0=-58.4627 +k=0.9999980000000001 +x_0=100000 +y_0=100000 +ellps=intl +units=m +no_defs"

Si cambiamos el dataset al sistema WGS84 (4326) directamente nos va a pasar lo siguiente:

veredas_to_wgs84 <- st_transform(veredas, 4326)

Como Gauss Krüger Buenos Aires usa el dátum Campo Inchauspe, lo que hacemos es convertir al viejo sistema Campo Inchauspe de Argentina, luego al sistema POSGAR en la correspondiente faja y finalmente al sistema WGS84:

El paso a paso aquí

Mundo

Algunas proyecciones en las que podemos presentar los mapas

WGS84

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



plot(mundo_wgs84)



World Mercator

mundo_mercator <- st_transform(mundo_wgs84, 54004)
st_crs(mundo_mercator)
## Coordinate Reference System:
##   EPSG: 54004 
##   proj4string: "+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"



Si ploteamos así:

Aquí limitamos las coordenadas:

ggplot() +  
    geom_sf(data = mundo_mercator)+ 
    labs(title = "Mundo Mercator") +
    coord_sf(ylim = c(18000000, -10000000)) + 
    theme_light()



World Mollweide

mundo_mollweide <- st_transform(mundo_wgs84, 54009)
st_crs(mundo_mollweide)
## Coordinate Reference System:
##   EPSG: 54009 
##   proj4string: "+proj=moll +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
plot(mundo_mollweide)