Combinando shapefiles + datasets con campos de coordenadas (lat y lon)

En éste ejercicio vamos a practicar la creación de mapas que combinan shapefiles y tablas con coordenadas de latitud y longitud.

Datos en tabla

Como ejemplo, cargamos un dataset con información sobre obra pública en la CABA.

obra_publica <- read.csv("https://data.buenosaires.gob.ar/api/files/observatorio-de-obras-urbanas.csv/download/csv", sep = ";")

obra_publica

Datos en shapefile

Ahora un shapefile con los límites de las comunas de la CABA

library(tidyverse)
library(sf)

comunas <- st_read('https://bitsandbricks.github.io/data/CABA_comunas.geojson')
## Reading layer `CABA_comunas' from data source `https://bitsandbricks.github.io/data/CABA_comunas.geojson' using driver `GeoJSON'
## Simple feature collection with 15 features and 4 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -58.53152 ymin: -34.70529 xmax: -58.33514 ymax: -34.52754
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs

Todo junto al mapa

La idea es proyectar todo junto: una capa con las comunas, y una capa con la posición de los arboles.

En este caso, una de nuestras fuentes de datos (las fronteras de las comunas) es un shapefile. Cuando se realizan gráficos que incluyen un shapefile, conviene pasar al mismo formato las otras fuentes de datos a proyectar. Eso evita problemas relacionados con mala alineacion de proyecciones, sistemas de coordenadas, y otras vicisitudes de la información geográfica.

Es fácil convertir un dataset “tradicional” en uno geográfico, cuando incluye columnas con la posición de latitud y de longitud de cada fila.

Usamos la función st_as_sf, que requiere:

  • un dataframe
  • los nombres de las columnas en el datafame con longitud y latitud (en ese orden)
  • el sistema de coordenadas bajo el cual se definieron las latitudes y longitudes

Los dos primeros requisitos son fáciles de cumplir; sólo es cuestión de cargar el dataframe y revisar su contenido para encontrar las columnas geográficas, que suelen llamarse lon y lat o X e Y. En nustro caso, se llaman “LNG” y “LAT”:

names(obra_publica)
##  [1] "ID"                        "ENTORNO"                  
##  [3] "NOMBRE"                    "ETAPA"                    
##  [5] "TIPO"                      "AREA_RESPONSABLE"         
##  [7] "DESCRIPCION"               "MONTO_CONTRATO"           
##  [9] "COMUNA"                    "BARRIO"                   
## [11] "CALLE_1"                   "SECCION"                  
## [13] "MANZANA"                   "PARCELA"                  
## [15] "DIRECCION"                 "LAT"                      
## [17] "LNG"                       "FECHA_INICIO"             
## [19] "FECHA_FIN_INICIAL"         "PLAZO_MESES"              
## [21] "PORCENTAJE_AVANCE"         "IMAGEN_1"                 
## [23] "IMAGEN_2"                  "IMAGEN_3"                 
## [25] "IMAGEN_4"                  "LICITACION_OFERTA_EMPRESA"
## [27] "LICITACION_ANIO"           "BENFICIARIOS"             
## [29] "MANO_OBRA"                 "COMPROMISO"               
## [31] "LINK_INTERNO"              "PLIEGO_DESCARGA"

Saber el sistema de referencia de coordenadas, o CRS, es mucho más difícil si no fue explicitado en el sitio donde conseguimos los datos. Pero hay un truco: en la gran mayoría de los casos, las coordenadas están expresadas con CRS WG84, alias “Mercator”, cuyo código de identificación es 4326.

Y por último, debemos quitar filas sin coordenadas (si las hubiera), ya que un shapefile no acepta valores de posición faltantes.

Por lo tanto necesitamos algo como

obra_publica <- obra_publica %>% 
    filter(!is.na(LNG), !is.na(LAT)) %>% 
    st_as_sf(coords = c("LNG", "LAT"), crs = 4326)

Resultado final

Con todo listo, es solo cuestión de sumar capas geom_sf().

Por ejemplo, mostrando los límites de las comunas y la posición (y tipo) de las obras públicas.

ggplot() + 
    geom_sf(data = comunas) +
    geom_sf(data = obra_publica, aes(color = TIPO))