Introducción a datos geográficos en R

Práctica Guiada

Reiniciar R

Cargamos las librerías. En esta ocasión usaremos por primera vez el paquete sf.

library(tidyverse)
library(ggthemes)
library(sf)

Representación de la información en mapas

Trabajaremos con datos del portal de datos abiertos de la Ciudad de Buenos Aires. En este caso, usamos los radios censales de la ciudad y los hemos descargado de https://bitsandbricks.github.io/data/CABA_rc.geojson.

La función st_read nos permite levantar un archivo de tipo .geojson. (también permite cargar otros con capas, pero no lo veremos en esta ocasión).

Vemos que el mismo cuenta con 3.554 features y 8 campos.
epsg (SRID): 4326 y proj4string: +proj=longlat +datum=WGS84 +no_defs refieren a que nuestros datos usan el sistema de coordenadas WGS84, también conocido por su código EPSG 4326 . Es el mismo que usan los sistemas GPS, Google Maps, y las aplicaciones de internet en general.

radios <- st_read("../fuentes/CABA_rc.geojson")

Reading layer `CABA_rc' from data source `/home/diego/Documents/0_GIT/2_Ed/intro_ds/fuentes/CABA_rc.geojson' using driver `GeoJSON'
Simple feature collection with 3554 features and 8 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -58.53092 ymin: -34.70574 xmax: -58.33455 ymax: -34.528
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs

Como los 8 campos son equivalentes a columnas/variables, podemos pedir un summary de la información contenida en radios.

summary(radios)

    RADIO_ID          BARRIO         COMUNA       POBLACION        VIVIENDAS     
 1_1_1  :   1   PALERMO  : 295   1      : 329   Min.   :   0.0   Min.   :   0.0  
 1_10_1 :   1   CABALLITO: 215   13     : 305   1st Qu.: 646.2   1st Qu.: 311.2  
 1_10_10:   1   RECOLETA : 198   14     : 295   Median : 786.0   Median : 377.0  
 1_10_11:   1   BALVANERA: 191   3      : 254   Mean   : 813.2   Mean   : 401.4  
 1_10_12:   1   FLORES   : 183   4      : 252   3rd Qu.: 928.0   3rd Qu.: 462.0  
 1_10_13:   1   BELGRANO : 170   7      : 250   Max.   :3945.0   Max.   :1405.0  
 (Other):3548   (Other)  :2302   (Other):1869                                    
    HOGARES        HOGARES_NBI        AREA_KM2                 geometry   
 Min.   :   0.0   Min.   :  0.00   Min.   :0.004468   MULTIPOLYGON :3554  
 1st Qu.: 259.0   1st Qu.:  2.00   1st Qu.:0.018626   epsg:4326    :   0  
 Median : 310.0   Median :  6.00   Median :0.035548   +proj=long...:   0  
 Mean   : 323.6   Mean   : 19.35   Mean   :0.057350                       
 3rd Qu.: 371.0   3rd Qu.: 23.00   3rd Qu.:0.062847                       
 Max.   :1093.0   Max.   :403.00   Max.   :3.804422                       
                                                                          

A la hora de graficar, podemos representar la información utilizando la combinación del famoso ggplot() en conjunto con geom_sf, que identifica la variable con la geometría y la grafica.

ggplot() + 
  geom_sf(data = radios)

Podemos colorear los radios censales de acuerdo a la cantidad de viviendas que hay en cada una de ellas, convocando a la variable en el parámetro de relleno. A su vez, podemos “jugar” con el color y ancho de los bordes.

ggplot(data = radios, aes(fill = VIVIENDAS)) + 
  geom_sf() 

ggplot(data = radios, aes(fill = VIVIENDAS)) + 
  geom_sf(color = "white") 

ggplot(data = radios, aes(fill = VIVIENDAS)) +
  geom_sf(color = "white", size=.1) 

O bien podemos omitir los bordes.

ggplot(data = radios, aes(fill = POBLACION)) + 
  geom_sf(color = NA)

También podemos usar el campo del relleno para realizar una transformación de los datos. Por ejemplo, podemos calcular la densidad de la población utilizando los datos de POBLACION y AREA_KM2. Aprovechamos para mejorar otras cuestiones estéticas 🌈

ggplot(data = radios, aes(fill = POBLACION/AREA_KM2)) + 
    geom_sf(color = NA) +
    scale_fill_viridis_c() +
    labs(title = "Densidad de población",
         subtitle = "Ciudad Autónoma de Buenos Aires",
         fill = "hab/km2") +
  theme_void()

También podemos rellenar de acuerdo a la variable BARRIO, que es una forma de graficar agregados. De todas formas, luego veremos como agruparlos con group_by.

ggplot(data = radios, aes(fill = BARRIO)) + 
  geom_sf(color = NA) +
  theme(legend.position = 'none')

También podemos representar rápidamente la información de alguna de las variables realizando un select() de la misma y la geometría, y pidiendo un plot() (el parámetro lwd define el ancho de las líneas de cada geometría, en este caso, los radios censales). En términos visuales, es como haber agregado por BARRIO la información, no?

radios %>% 
  select(BARRIO, geometry) %>% 
  plot(lwd=0.06)

Agrupando polígonos

Sin embargo, si lo que queremos es reconstruir los polígonos de los barrios, sólo necesitamos hacer un group_by y un summarise, y automáticamente en la columna geometry se crea el polígono combinado. Además, podemos hacer una agregación de las demás variables.

radios %>% 
  sample_n(5)

Simple feature collection with 5 features and 8 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -58.49039 ymin: -34.64358 xmax: -58.37329 ymax: -34.56128
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
  RADIO_ID        BARRIO COMUNA POBLACION VIVIENDAS HOGARES HOGARES_NBI
1  15_11_6     CHACARITA     15       798       288     266          22
2   2_20_7      RECOLETA      2       569       380     247           5
3   5_12_4       ALMAGRO      5       735       407     328           8
4   12_9_6 VILLA URQUIZA     12       960       452     377           3
5  4_12_11      BARRACAS      4       632       296     260           7
    AREA_KM2                       geometry
1 0.07511176 MULTIPOLYGON (((-58.43938 -...
2 0.02277754 MULTIPOLYGON (((-58.40963 -...
3 0.01722292 MULTIPOLYGON (((-58.42623 -...
4 0.07053921 MULTIPOLYGON (((-58.48944 -...
5 0.03611037 MULTIPOLYGON (((-58.3737 -3...

barrios_geo <- radios %>% 
    group_by(BARRIO) %>% 
    summarise(POBLACION = sum(POBLACION),
              VIVIENDAS = sum(VIVIENDAS),
              HOGARES = sum(HOGARES),
              HOGARES_NBI = sum(HOGARES_NBI),
              AREA_KM2 = sum(AREA_KM2))

barrios_geo %>% 
  sample_n(5)

Simple feature collection with 5 features and 6 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -58.53092 ymin: -34.67492 xmax: -58.43793 ymax: -34.55955
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs

Si pedimos un plot() del objeto, nos devuelve la representación geográfica con los datos de cada una de las variables.

plot(barrios_geo)

Si deseamos agregar otra capa al gráfico que incluya texto o etiquetas de los datos, contamos con geom_sf_text(), en la cual podemos también solicitar que se etiqueten sólo aquellos datos que cumplen con alguna condición.

ggplot(data = barrios_geo, aes(fill = POBLACION/AREA_KM2)) +
    geom_sf(color = NA) +
    geom_sf_text(data = barrios_geo %>% 
                   filter(POBLACION/AREA_KM2>25000), 
                 aes(label = BARRIO), size=3)+
  theme(legend.position = 'bottom')+
  scale_fill_viridis_c(option = 'C') +
  labs(title = "Densidad de población",
        subtitle = "Ciudad Autónoma de Buenos Aires",
        fill = "hab/km2")+
  theme_void()

Volcando en el mapa información de múltiples fuentes: Subtes

Ahora utilizaremos la información de líneas y estaciones de subte. Hemos descargado los .geojson de http://bitsandbricks.github.io/data/subte_lineas.geojson y http://bitsandbricks.github.io/data/subte_estaciones.geojson.

subte_lineas <- st_read('../fuentes/subte_lineas.geojson')

Reading layer `subte_lineas' from data source `/home/diego/Documents/0_GIT/2_Ed/intro_ds/fuentes/subte_lineas.geojson' using driver `GeoJSON'
Simple feature collection with 80 features and 2 fields
geometry type:  MULTILINESTRING
dimension:      XY
bbox:           xmin: -58.48639 ymin: -34.64331 xmax: -58.36993 ymax: -34.55564
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs

subte_lineas %>% 
  sample_n(5)

Simple feature collection with 5 features and 2 fields
geometry type:  MULTILINESTRING
dimension:      XY
bbox:           xmin: -58.46165 ymin: -34.64331 xmax: -58.37804 ymax: -34.57842
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
  ID LINEASUB                       geometry
1 16  LINEA C MULTILINESTRING ((-58.37953...
2  6  LINEA D MULTILINESTRING ((-58.4212 ...
3 70  LINEA A MULTILINESTRING ((-58.44865...
4 59  LINEA E MULTILINESTRING ((-58.45789...
5  5  LINEA D MULTILINESTRING ((-58.42571...

subte_estaciones <- st_read('../fuentes/subte_estaciones.geojson')

Reading layer `subte_estaciones' from data source `/home/diego/Documents/0_GIT/2_Ed/intro_ds/fuentes/subte_estaciones.geojson' using driver `GeoJSON'
Simple feature collection with 86 features and 3 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: -58.48639 ymin: -34.64331 xmax: -58.36993 ymax: -34.55564
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs

subte_estaciones %>% 
  sample_n(5)

Simple feature collection with 5 features and 3 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: -58.48101 ymin: -34.62192 xmax: -58.37992 ymax: -34.5778
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
  ID                 ESTACION LINEA                    geometry
1 74 DE LOS INCAS -PQUE. CHAS     B POINT (-58.47424 -34.58125)
2 58                 SAN JUAN     C POINT (-58.37992 -34.62192)
3 18            INDEPENDENCIA     C POINT (-58.38017 -34.61813)
4 80               ECHEVERR�A     B  POINT (-58.48101 -34.5778)
5 46            INDEPENDENCIA     E POINT (-58.38153 -34.61794)

Con geom_sf agregamos diferentes capas de información geográfica:

• los barrios: coloreados según su proporción de hogares con necesidades básicas insatisfechas, sin color de borde.
• las líneas de subte: líneas en amarillo
• las estaciones de subte: puntos en naranja

Alguna reflexión sobre los resultados?

ggplot() +
    geom_sf(data = barrios_geo, aes(fill = HOGARES_NBI/HOGARES), color = NA) +
    geom_sf(data = subte_lineas, color = "yellow") +
    geom_sf(data = subte_estaciones, color = "orange") +
  theme(legend.position = 'bottom')+
  scale_fill_viridis_c()+
    labs(title = "Sistema de transporte subterráneo (SUBTE)",
         subtitle = "Ciudad de Buenos Aires")

Mapeo de Palos borrachos rosados en Buenos Aires

Ahora utilizaremos información del portal de datos abiertos de la Ciudad de Buenos Aires sobre arbolado. Nos hemos quedado con una selección de dicha base de datos, conservando sólo a los árboles de tipo “palo borracho rosado”.

Notemos que en este caso la información está en formato .rds, es un dataframe donde los ejemplares están georreferenciados con variables long-lat, no hay una geometría como en los casos anteriores.

palos_borrachos <- read_rds('../fuentes/arbolado_palo_borracho.rds')

palos_borrachos %>% 
  sample_n(5)

Transformación de datos tabulados a datos sf

Para transformar los datos a tipo sf, tenemos la función st_as_sf(), a la cual le indicamos en el parámetro coords las variables que contienen las coordenadas. Con la función st_set_crs() podemos indicar el sistema de coordenadas de referencia.

palos_borrachos <- st_as_sf(palos_borrachos, coords = c('long','lat')) %>% 
    st_set_crs(4326)

palos_borrachos %>% 
  sample_n(5)

Simple feature collection with 5 features and 18 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: -58.51574 ymin: -34.66754 xmax: -58.4335 ymax: -34.54007
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs

Ahora podemos mapear los palos borrachos en los barrios de CABA que ya estuvimos trabajando!

ggplot() +
  geom_sf(data = barrios_geo) +
  geom_sf(data = palos_borrachos, color = "darkgreen", size=.5) +
  labs(title = 'localización de los palos borrachos de Buenos Aires')+
  theme_void()+
  theme(legend.position = 'none')

Join espacial

Ahora probemos la forma de unir la información contenida en palos_borrachos con la de barrios_geo. Para eso usamos st_join(), pero… en qué orden deberíamos unirlos?

palos_borrachos %>% 
  st_join(barrios_geo) %>% 
  sample_n(5)

Simple feature collection with 5 features and 24 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: -58.4673 ymin: -34.63316 xmax: -58.36219 ymax: -34.59118
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs

barrios_geo %>% 
  st_join(palos_borrachos) %>% 
  sample_n(5)

Simple feature collection with 5 features and 24 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -58.50844 ymin: -34.62557 xmax: -58.33925 ymax: -34.53199
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs

Si queremos visualizar la densidad de borrachos por barrio:

palos_borrachos_barrio <- barrios_geo %>% 
  st_join(palos_borrachos)

palos_borrachos_barrio %>% 
  group_by(BARRIO) %>% 
  summarise(densidad=n()/unique(AREA_KM2)) %>% 
  ggplot(aes(fill=densidad, label=BARRIO))+
  geom_sf()+
  geom_sf_text( color='white')+
  theme_void()