INFORMACIÓN GEOGRÁFICA Y MAPAS

Llegó el momento de visualizar información en un mapa. Pero antes de empezar, ¿A que nos referimos cuando hablamos de Sistemas de Información Geográfica o SIG? Bueno, nos referimos a las herramientas informáticas que nos permiten ubicar y analizar un conjunto de datos en lugares específicos del territorio (georreferenciar).

Para poder hacer uso de los SIG necesitamos contar con información geográfica en nuestro dataset, es decir, que además de la información que ya vimos que puede haber en un dataset tradicional, se sume un componente espacial en cada registro (partido, barrio, manzana, calle o directamente las coordenadas X e Y).

Las herramientas que nos permitirán visualizar toda esta información serán los mapas, que son nada más y nada menos que representaciones planas, reducidas y simplificadas de la tierra que nos dan la posibilidad cruzar y relacionar datos en el espacio. Es decir que, mantienen una relación ordenada en el traspaso de puntos ubicados en la superficie curva de la tierra a puntos ubicados en la superficie plana de los mapas. Esto es posible a partir del uso de sistemas de coordenadas proyectadas.

En R hay varios paquetes de funciones que nos permiten manipular este tipo de información, entre los que se encuentra sf, que lo aprenderemos hoy. Para comenzar a utilizarlo vamos a tener que instalarlo y luego activarlo con library() al igual que lo veníamos haciendo con tidyverse:

#install.packages(tidyverse)
library(tidyverse)
#install.packages(sf)
library(sf)

Como verán, activamos los 2 paquetes porque ambos presentan funciones que son necesarias a la hora de mapear información.

Analizar datos espaciales

Arranquemos cargando nuestro dataset espacial o shape con las geometrías correspondientes a todos los partidos de AMBA.

Recomendación: Al descargar el shape deberán moverlo de la carpeta “Descargas” a la carpeta del Proyecto donde estan trabajando.

Para poder cargar nuestros datos espaciales en formato geoJSON utilizaremos la función st_read() de la siguiente forma:

partidos_amba <- st_read("partidos_amba.geojson", stringsAsFactors = TRUE)
## Reading layer `partidos_amba' from data source `D:\00-TRABAJO\04-FLACSO\partidos_amba.geojson' using driver `GeoJSON'
## Simple feature collection with 48 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -59.3392 ymin: -35.23893 xmax: -57.70946 ymax: -34.23007
## Geodetic CRS:  WGS 84

Veamos que información contiene:

summary(partidos_amba)
##              nombre   provincia    area_km2               geometry 
##  Almirante Brown: 1   CABA:15   Min.   :  6.30   MULTIPOLYGON :48  
##  Avellaneda     : 1   GBA :33   1st Qu.: 17.72   epsg:4326    : 0  
##  Berazategui    : 1             Median : 55.91   +proj=long...: 0  
##  Berisso        : 1             Mean   :140.59                     
##  Comuna 1       : 1             3rd Qu.:177.07                     
##  Comuna 10      : 1             Max.   :889.49                     
##  (Other)        :42

Las primeras 3 columnas presentan el nombre del partido, la provincia y el área en km2. Hasta acá son datos muy similares a las que ya veníamos encontrando en los dataset tradicionales; sin embargo, aparece una 4ta columna llamada “geometry” que hasta ahora no la habíamos visto y es donde se aloja la geometría de cada uno de los registros. La información de este campo es la que hace que el dataset sea espacial.

Si queda alguna duda, podemos utilizar la función class() para ver con que tipo de datos estamos trabajando:

class(partidos_amba)
## [1] "sf"         "data.frame"

Efectivamente partidos_amba, es un dataset espacial u objeto del tipo “sf”, que hace referencia a “simple features” por estar compuesto de geometrías bidimensionales (polígono, punto, línea, multipunto, multilínea, etc.).

Para poder visualizar toda esta información plasmada en un mapa vamos a utilizar nuevamente ggplot() pero como en esta oportunidad queremos sumar capas geográficas (sf) trabajaremos con geom_sf(). Veamos un ejemplo:

ggplot(partidos_amba)+
  geom_sf()

Como habrán notado, la lógica en la estructura del chunk que utilizamos para hacer este mapa es la misma del capítulo anterior, donde dentro del ggplot() asignamos el dataset y luego sumamos la capa a graficar/mapear (en este caso un objeto sf).

En el mapa podemos ver como las geometrías de cada registro del dataset espacial son polígonos que representan los partidos de AMBA. Probemos agregar algún atributo estético aes():

ggplot(partidos_amba)+
  geom_sf(aes(fill=area_km2))

Ya tenemos nuestro primer mapa coroplético. Probemos ajustando algunas cuestiones estéticas: cambiemos la paleta (y su escala) y el color del borde de los polígonos.

ggplot(partidos_amba)+
  geom_sf(aes(fill=area_km2), color="white")+
  scale_fill_viridis_c(breaks=c(0,250,500,750,1000))

Pero nuestro shape tiene poca información no? No podemos hacer muchos análisis si solo contamos con la variable de superficie, así que seamos proactivos y agreguemos más datos!

Cruzar datos tradicionales y espaciales

Es muy común que a la hora de trabajar con datos no encontremos toda la información que necesitamos en un solo dataset. Frente a esto, lo que se hace es unir/cruzar datos provenientes de diferentes fuentes de información. En este apartado veremos como realizar estos cruces entre datos tradicionales y espaciales.

Para poder unir 2 set de datos ambos tienen que tener algo en común (alguna columna como por ejemplo un ID), sino sería imposible que R entienda que tiene que unir con qué. Por lo tanto, antes de unirlos podríamos necesitar manipularlos y realizarles diferentes tipos de transformaciones que nos permitan llegar a generar esta variable en común.

Por ejemplo en nuestro caso, si queremos unir algún dato a nuestro dataset espacial de partidos, deberíamos tener un valor único por cada nombre de partido (que funciona como el ID de mi dataset espacial). Probemos esto con algunos datos del Censo 2010:

partidos_censo2010 <- read.csv("partidos_censo2010.csv", stringsAsFactors = TRUE)

Veamos que tenemos:

summary(partidos_censo2010)
##      codigo                 nombre      pob_2010          viv_2010     
##  Min.   :2001   Almirante Brown: 1   Min.   :  56729   Min.   : 19287  
##  1st Qu.:2013   Avellaneda     : 1   1st Qu.: 174013   1st Qu.: 64864  
##  Median :6319   Berazategui    : 1   Median : 223281   Median : 93389  
##  Mean   :5075   Berisso        : 1   Mean   : 301108   Mean   :104104  
##  3rd Qu.:6544   Comuna 1       : 1   3rd Qu.: 340723   3rd Qu.:123359  
##  Max.   :6861   Comuna 10      : 1   Max.   :1775816   Max.   :447306  
##                 (Other)        :42                                     
##     hog_2010     
##  Min.   : 17116  
##  1st Qu.: 59537  
##  Median : 81055  
##  Mean   : 95612  
##  3rd Qu.:109566  
##  Max.   :484909  
## 

Bien, es un dataset tradicional con 5 variables o columnas: codigo, nombre, población 2010, cantidad de viviendas 2010 y cantidad de hogares 2010. Podemos analizar un poco los valores y ver por ejemplo que, en promedio los partidos de AMBA tienen 301108 habitantes, y que el más poblado tiene 1775816 y el menos poblado 56729.

Pero el detalle más importante es que hay una columna llamada “nombre” que tiene el nombre de nuestros partidos, y eso es una buena noticia ya que ella será quien nos permitirá unir ambos dataset (el espacial y el tradicional). Para esto utilizaremos la función llamada left_join():

partidos_amba <- left_join(partidos_amba, partidos_censo2010, by="nombre")

Revisemos el resultado!

summary(partidos_amba)
##              nombre   provincia    area_km2          codigo    
##  Almirante Brown: 1   CABA:15   Min.   :  6.30   Min.   :2001  
##  Avellaneda     : 1   GBA :33   1st Qu.: 17.72   1st Qu.:2013  
##  Berazategui    : 1             Median : 55.91   Median :6319  
##  Berisso        : 1             Mean   :140.59   Mean   :5075  
##  Comuna 1       : 1             3rd Qu.:177.07   3rd Qu.:6544  
##  Comuna 10      : 1             Max.   :889.49   Max.   :6861  
##  (Other)        :42                                            
##     pob_2010          viv_2010         hog_2010               geometry 
##  Min.   :  56729   Min.   : 19287   Min.   : 17116   MULTIPOLYGON :48  
##  1st Qu.: 174013   1st Qu.: 64864   1st Qu.: 59537   epsg:4326    : 0  
##  Median : 223281   Median : 93389   Median : 81055   +proj=long...: 0  
##  Mean   : 301108   Mean   :104104   Mean   : 95612                     
##  3rd Qu.: 340723   3rd Qu.:123359   3rd Qu.:109566                     
##  Max.   :1775816   Max.   :447306   Max.   :484909                     
## 

Efectivamente, se sumaron 4 nuevas columnas a mi dataset espacial: codigo, pob_2010, viv_2010 y hog_2010. ¿Qué esperamos? ¡Veamos alguna de ellas en un mapa!

ggplot(partidos_amba)+
  geom_sf(aes(fill=pob_2010), color="white")+
  scale_fill_viridis_c()

Cómo veran, hay un partido que resalta ya que es el más poblado: La Matanza. Pero no nos quedemos con esa idea, y para que estos datos sean comparables entre todos los partidos, es necesario que calculemos una densidad relacionando la cantidad de habitantes y la superficie (km2) de cada uno. También aprovechemos y mejoremos la estética agregando etiquetas (título, subtítulo, etc) y cambiando el aspecto (theme):

ggplot(partidos_amba)+
  geom_sf(aes(fill=pob_2010/area_km2), color=NA)+
    labs(title = "Población por Partido",
         subtitle = "Censo 2010",
         fill = "Población/km2",
         caption= "Fuente: INDEC") +
  scale_fill_distiller(palette = "Spectral") +
  theme_light()

Ahora el panorama cambió un poco y vemos como algunas comunas de CABA son las más densas. ¿Cómo hacemos si queremos ver este mapa con un “zoom” solo de CABA?

Tenemos que filtrar el dataset de la siguiente forma:

ggplot(partidos_amba %>% filter(provincia=="CABA"))+
  geom_sf(aes(fill=pob_2010/area_km2), color=NA)+
    labs(title = "Población por Partido",
         subtitle = "Censo 2010",
         fill = "Población/km2",
         caption= "Fuente: INDEC") +
  scale_fill_viridis_c(option = "plasma", direction=-1) +
  theme_light()

Ahora probemos con otra variable: cantidad de viviendas por km2.

ggplot(partidos_amba %>% filter(provincia=="CABA"))+
  geom_sf(aes(fill=viv_2010/area_km2), color=NA)+
    labs(title = "Viviendas por Partido",
         subtitle = "Censo 2010",
         fill = "Vivienda/km2",
         caption= "Fuente: INDEC") +
  scale_fill_viridis_c(option = "plasma", direction=-1) +
  theme_light()

Pongamos algunas etiquetas para poder diferenciar las comunas rápidamente!

ggplot(partidos_amba %>% filter(provincia=="CABA"))+
  geom_sf(aes(fill=viv_2010/area_km2), color=NA)+
  geom_sf_text(aes(label = nombre), size=2) +
    labs(title = "Densidad de Viviendas por Partido",
         subtitle = "Censo 2010",
         fill = "Vivienda/km2",
         caption= "Fuente: INDEC") +
  scale_fill_viridis_c(option = "plasma", direction=-1) +
  theme_light()
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data

Podemos ver que aparece la Comuna 6 como la más densa seguida por la 5. Sin embargo hay las densidades van de 4000 a 12000 habitantes por km2. Veamos cuáles son las más densas, por ejemplo las que tienen más de 8000 habitantes por km2:

ggplot(partidos_amba %>% filter(provincia=="CABA" & viv_2010/area_km2>=8000))+
  geom_sf(aes(fill=viv_2010/area_km2), color=NA)+
  geom_sf_text(aes(label = nombre), size=2) +
    labs(title = "Densidad de viviendas por Partido",
         subtitle = "Partidos con más de 8000 viviendas por km2.",
         fill = "Vivienda/km2",
         caption= "Fuente: Censo 2010. INDEC") +
  scale_fill_viridis_c(option = "plasma", direction=-1) +
  theme_light()
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data

Y por último agreguémosle el mapa de fondo así las comunas no quedan flotando:

ggplot()+
  geom_sf(data=partidos_amba %>% filter(provincia=="CABA"))+
  geom_sf(data=partidos_amba %>% filter(provincia=="CABA" & viv_2010/area_km2>=8000), aes(fill=viv_2010/area_km2))+
  geom_sf_text(data=partidos_amba %>% filter(provincia=="CABA" & viv_2010/area_km2>=8000), aes(label = nombre), size=2) +
    labs(title = "Densidad de viviendas por Partido",
         subtitle = "Partidos con más de 8000 viviendas por km2.",
         fill = "Vivienda/km2",
         caption= "Fuente: Censo 2010. INDEC") +
  scale_fill_viridis_c(option = "plasma", direction=-1) +
  theme_light()
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data

Variables Categóricas en el Mapa

Hasta acá todos los mapas que desarrollamos tomaron color a partir de una variable numérica (cantidad o densidad por km2), pero ¿Qué pasa si lo que queremos mapear es una variable categórica, cómo por ejemplo diferenciar aquellas que tienen más de 300000 habitantes por km2 de las que tienen menos?

Para resolver esta incógnita debemos generar una nueva variable que contenga la información:

partidos_amba <- partidos_amba %>%
  mutate(CATEGORIA=as.factor(if_else(pob_2010>=300000, "MAYOR DE 300K", "MENOR DE 300K")))
summary(partidos_amba)
##              nombre   provincia    area_km2          codigo    
##  Almirante Brown: 1   CABA:15   Min.   :  6.30   Min.   :2001  
##  Avellaneda     : 1   GBA :33   1st Qu.: 17.72   1st Qu.:2013  
##  Berazategui    : 1             Median : 55.91   Median :6319  
##  Berisso        : 1             Mean   :140.59   Mean   :5075  
##  Comuna 1       : 1             3rd Qu.:177.07   3rd Qu.:6544  
##  Comuna 10      : 1             Max.   :889.49   Max.   :6861  
##  (Other)        :42                                            
##     pob_2010          viv_2010         hog_2010               geometry 
##  Min.   :  56729   Min.   : 19287   Min.   : 17116   MULTIPOLYGON :48  
##  1st Qu.: 174013   1st Qu.: 64864   1st Qu.: 59537   epsg:4326    : 0  
##  Median : 223281   Median : 93389   Median : 81055   +proj=long...: 0  
##  Mean   : 301108   Mean   :104104   Mean   : 95612                     
##  3rd Qu.: 340723   3rd Qu.:123359   3rd Qu.:109566                     
##  Max.   :1775816   Max.   :447306   Max.   :484909                     
##                                                                        
##          CATEGORIA 
##  MAYOR DE 300K:17  
##  MENOR DE 300K:31  
##                    
##                    
##                    
##                    
## 

Excelente ya tenemos nuestra variable categórica que nos indica que 17 partidos tienen más de 300000 habitantes y 31 tienen menos. A mapear!

ggplot(partidos_amba)+
  geom_sf(aes(fill=CATEGORIA))

Nada mal no? Modifiquemos su estética como ya aprendimos:

ggplot(partidos_amba)+
  geom_sf(aes(fill=CATEGORIA), color="black") +
    labs(title = "Población por Partido",
         subtitle = "Censo 2010",
         fill = "Categoría",
         caption= "Fuente: INDEC") +
  scale_fill_manual(values=c("darkseagreen1", "cyan4"))+
  theme_light()

Unir datos a partir de sus ubicaciones espaciales

Hasta acá hicimos una introducción a qué son los datos espaciales y vimos los siguientes temas:

  • Abrir y analizar datos espaciales: El paquete sf.

  • Unir datos espaciales y datos tradicionales: El left_join().

Pero todavía nos queda camino por recorrer en el mundo de los datos espaciales así que a continuación abordaremos algunos temas más. Por ejemplo, reiteradas veces nos vamos a encontrar con que queremos unir 2 dataset que no tienen columnas en común y es imposible generarlas; pero hay algo que si tienen en común y que nos permitirá unirlos: su ubicación en el espacio. A esto nos referimos cuando hablamos de “Spatial Join” o “Unión Espacial” y en R lo haremos a partir de la función st_join() que forma parte del ya conocido paquete sf.

Pero vamos por partes ya que aún nos quedan algunos pasos previos a la unión espacial. Continuemos trabajando con el dataset de partidos_amba y carguemos una nueva base de datos: un dataset tradicional que tenga como conexión al territorio las coordenadas geográficas x e y (longitud y latitud).

Para eso carguemos el dataset de equipamiento cultural (https://data.world/angie-scetta/equipamiento-cultural) que fue generado a partir de los datos descargados del Portal de Datos Abiertos de Argentina (https://datos.gob.ar/dataset/cultura-espacios-culturales):

cultura <- read.csv("equipamiento_cultural.csv", stringsAsFactors = TRUE)

Si queda alguna duda, podemos utilizar la función class() para ver con que tipo de datos estamos trabajando:

class(cultura)
## [1] "data.frame"

Tal como lo imaginábamos, equipamiento_cultural es un simple dataframe. Analicemos su contenido:

summary(cultura)
##  espacio_cultural_id                           categoria  
##  Min.   :     9      Bibliotecas Especializadas     :398  
##  1st Qu.:110756      Bibliotecas Populares          :142  
##  Median :120588      Casas del Bicentenario         :  2  
##  Mean   :117078      Monumentos y Lugares Historicos:256  
##  3rd Qu.:160381      Museos                         : 90  
##  Max.   :301038      Salas de Cine                  : 38  
##                      Salas de Teatro                :321  
##                               nombre_establecimiento    latitud      
##  Biblioteca Central                      :  20       Min.   :-35.77  
##  Centro de Documentación                 :  10       1st Qu.:-34.63  
##  Biblioteca Popular Juan Bautista Alberdi:   6       Median :-34.60  
##  Cinemark                                :   4       Mean   :-34.63  
##  Atlas Cine                              :   3       3rd Qu.:-34.59  
##  Biblioteca del Docente                  :   3       Max.   :-34.24  
##  (Other)                                 :1201                       
##     longitud     
##  Min.   :-59.90  
##  1st Qu.:-58.47  
##  Median :-58.40  
##  Mean   :-58.46  
##  3rd Qu.:-58.38  
##  Max.   :-58.00  
## 

Hay 5 columnas que refieren a ID del equipamiento, categoría a la que pertenece, nombre del establecimiento y ubicación en el espacio (longitud y latitud). Por ejemplo, podemos ver que la categoría que más abunda en nuestro dataset es Bibliotecas Especializadas y la que menos es Casas del Bicentenario. Sin embargo, no hay ninguna variable que tenga la misma información que el dataset espacial y por lo tanto no vamos a poder unirlos tan fácilmente con left_join() como hicimos la clase pasada. Pero aprovechemos que tenemos la longitud y latitud y hagamos un gráfico de puntos para ver como lucen en el espacio:

ggplot(cultura) +
  geom_point(aes(x=longitud, y=latitud))

Con mucha imaginación uno puede darse cuenta que los datos tienen la forma de AMBA, pero está faltando información que nos facilite la lectura del mismo. Para esto agreguemos nuestro dataset espacial de fondo así contextualizamos un poco la información, pero siempre teniendo en cuenta 2 cosas:

  1. Cuando utilizamos 2 capas (en este caso geom_sf() y geom_point()), tenemos que asignarle a cada una el dataset correspondiente dentro de la capa y ggplot() queda “vacío”.

  2. Las capas se grafican en el mismo orden que se escribe el código, es decir que la que agregamos primero (en este caso partidos_amba) será el fondo de la que agreguemos luego (cultura).

Veamos esto en detalle:

ggplot()+
  geom_sf(data=partidos_amba, fill="gray75", color="white")+
  geom_point(data=cultura, aes(x=longitud, y=latitud, color=categoria), alpha=0.5)+
  labs(title="Área Metropolitana de Buenos Aires")

Ahora si! Mucho más clara la ubicación de los equipamientos culturales. Parece que tienen algo en común el dataset de partidos y el de culutra, ¿No? Se solapan en el espacio y eso es justo lo que necesitábamos para poder realizar nuestra unión espacial!

Pero antes de seguir, vale aclarar que el segundo dataset que cargamos sigue siendo “tradicional” y para poder hacer una unión espacial, ambos deben ser espaciales. Por lo tanto, tenemos que transformarlo a este formato, pero no desesperen que por suerte el paquete sf tiene una función muy sencilla de utilizar que se llama st_as_sf():

cultura_geo <- cultura %>% 
    st_as_sf(coords = c("longitud", "latitud"), crs = 4326)

Revisemos que efectivamente haya ocurrido lo que queríamos:

class(cultura_geo)
## [1] "sf"         "data.frame"

Perfecto, cultura_geo es un dataset espacial u objeto del tipo “sf”, que hace referencia a “simple features” por estar compuesto de geometrías bidimensionales (en este caso puntos).

Veamos si se agregó el campo geometry:

head(cultura_geo)
## Simple feature collection with 6 features and 3 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -58.56202 ymin: -34.67005 xmax: -58.50309 ymax: -34.44614
## Geodetic CRS:  WGS 84
##   espacio_cultural_id                  categoria
## 1              160850 Bibliotecas Especializadas
## 2              160860 Bibliotecas Especializadas
## 3              160019 Bibliotecas Especializadas
## 4              160021 Bibliotecas Especializadas
## 5              160023 Bibliotecas Especializadas
## 6              160056 Bibliotecas Especializadas
##                                                       nombre_establecimiento
## 1                              Biblioteca del Instituto de Tecnología Minera
## 2                                                          Biblioteca - Unlm
## 3                                                   Departamento Informática
## 4 Biblioteca - Centro de Investigaciones En Láseres y Aplicaciones - Conicet
## 5                Biblioteca - Centro de Investigaciones En Sólidos - Conicet
## 6                                                    Biblioteca Max Von Buch
##                      geometry
## 1  POINT (-58.5097 -34.57088)
## 2 POINT (-58.56202 -34.67005)
## 3 POINT (-58.50309 -34.55512)
## 4 POINT (-58.50309 -34.55512)
## 5  POINT (-58.5091 -34.55351)
## 6 POINT (-58.52993 -34.44613)

Tal como esperábamos, nuestro dataset de cultura se transformó en espacial y eso lo vemos con el reemplazo de sus columnas X e Y por una única columna llamada “geometry”. Hagamos un mapa y veamos que pasa:

ggplot(cultura_geo)+
  geom_sf()

Y probemos dejando los partidos por debajo:

ggplot()+
  geom_sf(data=partidos_amba, fill="gray75", color="white")+
  geom_sf(data=cultura_geo, aes(color=categoria), alpha=0.5)+
  labs(title="Área Metropolitana de Buenos Aires")

A simple vista los datos son iguales que cuando hicimos el geom_point() pero la diferencia es que ahora son espaciales y por lo tanto… ¡Ya estamos en condiciones de unirlos con otro dataset espacial!

Llegó momento, hagamos nuestra unión espacial con st_join() para que, cada registro de cultura_geo sume una nueva columna que indique a que partido pertenece:

cultura_geo <- st_join(cultura_geo, partidos_amba)

Cabe destacar que, el st_join(), al igual que en el left_join(), del lado izquierdo se pone el dataset al que queremos unirle información nueva y del lado derecho el dataset del que queremos extraer la información.

Revisemos el resultado!

summary(cultura_geo)
##  espacio_cultural_id                           categoria  
##  Min.   :     9      Bibliotecas Especializadas     :398  
##  1st Qu.:110756      Bibliotecas Populares          :142  
##  Median :120588      Casas del Bicentenario         :  2  
##  Mean   :117078      Monumentos y Lugares Historicos:256  
##  3rd Qu.:160381      Museos                         : 90  
##  Max.   :301038      Salas de Cine                  : 38  
##                      Salas de Teatro                :321  
##                               nombre_establecimiento       nombre    provincia 
##  Biblioteca Central                      :  20       Comuna 1 :306   CABA:828  
##  Centro de Documentación                 :  10       Comuna 2 :135   GBA :368  
##  Biblioteca Popular Juan Bautista Alberdi:   6       Comuna 14: 75   NA's: 51  
##  Cinemark                                :   4       Comuna 3 : 67             
##  Atlas Cine                              :   3       Comuna 13: 55             
##  Biblioteca del Docente                  :   3       (Other)  :558             
##  (Other)                                 :1201       NA's     : 51             
##     area_km2          codigo        pob_2010          viv_2010     
##  Min.   :  6.30   Min.   :2001   Min.   :  87185   Min.   : 29309  
##  1st Qu.: 14.32   1st Qu.:2001   1st Qu.: 187537   1st Qu.: 82926  
##  Median : 18.20   Median :2007   Median : 205886   Median :107967  
##  Mean   : 59.05   Mean   :3386   Mean   : 278731   Mean   :112325  
##  3rd Qu.: 50.30   3rd Qu.:6270   3rd Qu.: 276190   3rd Qu.:130771  
##  Max.   :889.49   Max.   :6861   Max.   :1775816   Max.   :447306  
##  NA's   :51       NA's   :51     NA's   :51        NA's   :51      
##     hog_2010              CATEGORIA            geometry   
##  Min.   : 24926   MAYOR DE 300K:265   POINT        :1247  
##  1st Qu.: 76455   MENOR DE 300K:931   epsg:4326    :   0  
##  Median : 84468   NA's         : 51   +proj=long...:   0  
##  Mean   :100337                                           
##  3rd Qu.:102918                                           
##  Max.   :484909                                           
##  NA's   :51

Como verán, a cada registro se le unió una columna del nombre del partido con el que solapan. Pero acá aparece algo nuevo, muchos valores nulos o NA´s que hacen referencia a aquellos equipamientos culturales ubicados fuera de los partidos de AMBA (es decir que no solaparon). Veamos esto en un mapa así queda más claro:

ggplot()+
  geom_sf(data=partidos_amba)+
  geom_sf(data=cultura_geo, aes(color=nombre), alpha=0.75, show.legend = FALSE)

Se ve como todos los puntos que se ubican en AMBA tienen color según sus partidos, pero los puntos fuera de AMBA quedaron grises porque su valor es NA.

Para “limpiar” esto filtremos todos los registros que tienen partido asignado:

cultura_geo <- cultura_geo %>%
  filter(!is.na(nombre))

Visualicemos nuevamente:

ggplot()+
  geom_sf(data=partidos_amba)+
  geom_sf(data=cultura_geo, aes(color=nombre), alpha=0.75, show.legend = FALSE)

Ya tenemos solo los registros de AMBA. ¿Pero ahora cómo hacemos si queremos realizar un mapa coroplético como los de la clase anterior pero coloreando los partidos según la cantidad de equipamiento cultural que tienen? Nuevamente manipulamos el set de datos con group_by() y luego utilizamos left_join():

cultura_geo <- cultura_geo %>%
  group_by(nombre) %>%
  summarise(cantidad=n())
head(cultura_geo)
## Simple feature collection with 6 features and 2 fields
## Geometry type: MULTIPOINT
## Dimension:     XY
## Bounding box:  xmin: -58.5093 ymin: -34.89443 xmax: -58.17487 ymax: -34.58117
## Geodetic CRS:  WGS 84
## # A tibble: 6 x 3
##   nombre       cantidad                                                 geometry
##   <fct>           <int>                                         <MULTIPOINT [°]>
## 1 Almirante B~       13 ((-58.4024 -34.80811), (-58.39952 -34.89443), (-58.3967~
## 2 Avellaneda         42 ((-58.39811 -34.66852), (-58.38907 -34.66918), (-58.382~
## 3 Berazategui         4 ((-58.25725 -34.78558), (-58.18802 -34.82783), (-58.174~
## 4 Comuna 1          306 ((-58.39269 -34.60922), (-58.39266 -34.60933), (-58.392~
## 5 Comuna 10           4 ((-58.5093 -34.61965), (-58.50573 -34.61458), (-58.4881~
## 6 Comuna 11          15 ((-58.50863 -34.60288), (-58.5065 -34.60821), (-58.5055~

Pero como vimos al principio de la clase, si usamos left_join(), uno de los dataset no tiene que ser espacial. Para esto, transformemos nuestro dataset espacial cultura_geo a un dataset tradicional, quitándole la geometría con st_set_geometry():

cultura_geo <- cultura_geo %>%
  st_set_geometry(NULL)
head(cultura_geo)
## # A tibble: 6 x 2
##   nombre          cantidad
##   <fct>              <int>
## 1 Almirante Brown       13
## 2 Avellaneda            42
## 3 Berazategui            4
## 4 Comuna 1             306
## 5 Comuna 10              4
## 6 Comuna 11             15

Ahora si, solo tiene 2 columnas con el nombre del partido y la cantidad de equipamientos, pero ya no tiene geometría. Procedamos a realizar la unión:

partidos_amba <- left_join(partidos_amba, cultura_geo, by="nombre")

Y a mapear las cantidades:

ggplot()+
  geom_sf(data=partidos_amba, aes(fill=cantidad))

Ya tenemos un mapa de cantidad de equipamientos culturales por partido, y podemos ver como la mayor cantidad de concentra en la Comuna 1 de CABA. Mejoremos un poco su estética y calculemos la densidad de equipamientos para ver como da:

ggplot()+
  geom_sf(data=partidos_amba, aes(fill=cantidad/area_km2))+
    labs(title = "Equipamientos Culturales por km2",
         subtitle = "Partidos de AMBA",
         fill = "Cantidad/km2",
         caption= "Fuente: Portal de Datos Abiertos Argentina") +
  scale_fill_distiller(palette = "YlOrRd", direction = 1) +
  theme_light()

Los partidos de AMBA se ven todos del mismo color, saquemos AMBA para que se ajuste la escala y podamos ver más detalle en los partidos:

ggplot()+
  geom_sf(data=filter(partidos_amba, provincia=="GBA"), aes(fill=cantidad/area_km2))+
    labs(title = "Equipamientos Culturales por km2",
         subtitle = "Partidos de GBA",
         fill = "Cantidad/km2",
         caption= "Fuente: Portal de Datos Abiertos Argentina") +
  scale_fill_distiller(palette = "YlOrRd", direction = 1) +
  theme_light()

Ahora si, vemos que los partidos del primer cordon son los que tienen la mayor densidad y va disminuyendo a medida que nos alejamos de CABA.

Y miremos solo también solo CABA:

ggplot()+
  geom_sf(data=filter(partidos_amba, provincia=="CABA"), aes(fill=cantidad/area_km2))+
    labs(title = "Equipamientos Culturales por km2",
         subtitle = "Comunas de CABA",
         fill = "Cantidad/km2",
         caption= "Fuente: Portal de Datos Abiertos Argentina") +
  scale_fill_distiller(palette = "YlOrRd", direction = 1) +
  theme_light()

En el mapa anterior se puede observar claramente que las comunas con mayor densidad de equipamiento cultural son la 5 y la 1.

Agregar geometrías extra

Muchas veces nos vamos a encontrar con que queremos sumar más shapes o capas para darle contexto a nuestro mapa. ¿Cómo lo hacemos? Fácil! agregando más geom_sf() al ggplot. Veamos un ejemplo y carguemos el geoJSON de las líneas de FFCC:

lineas_ffcc <- st_read("lineas_ffcc.geojson")
## Reading layer `lineas_ffcc' from data source `D:\00-TRABAJO\04-FLACSO\lineas_ffcc.geojson' using driver `GeoJSON'
## Simple feature collection with 23 features and 5 fields
## Geometry type: MULTILINESTRING
## Dimension:     XY
## Bounding box:  xmin: -59.2763 ymin: -34.98144 xmax: -57.94964 ymax: -34.2941
## Geodetic CRS:  WGS 84

Y agreguemos las líneas sobre el último mapa de AMBA que hicimos:

ggplot()+
  geom_sf(data=partidos_amba, aes(fill=cantidad/area_km2))+
    labs(title = "Equipamientos Culturales por km2",
         subtitle = "Partidos de AMBA",
         fill = "Cantidad/km2",
         caption= "Fuente: Portal de Datos Abiertos Argentina") +
  scale_fill_distiller(palette = "YlOrRd", direction = 1) +
  geom_sf(data=lineas_ffcc)+
  theme_light()

Ahora mejoremos un poco la estética y resaltemos las líneas de FFCC:

ggplot()+
  geom_sf(data=partidos_amba) +
  geom_sf(data=lineas_ffcc, aes(color=as.factor(LINEA)), size=1) +
    labs(title = "Líneas de FFCC",
         subtitle = "AMBA",
         color = "Línea",
         caption= "Fuente: Min de Transporte") +
  theme_light()

Y listo! Ahora les toca practicar a uds!

Ejercicios de práctica

  1. Elegir una ciudad de cualquier parte del mundo que les interese y que disponga de un portal de datos abiertos que ofrezca un shapefile (formato geojson o shp) de polígonos con alguna división geográfica: barrios, comunas, etc.

  2. Realizar todas las transformaciones que sean necesarias para poder desarrollar al menos un ejemplo de cada uno de los siguientes mapas: 2.a. Mapa coroplético que muestre la distribución geográfica de una variable numérica. 2.b. Mapa coroplético que muestre la distribución geográfica de una variable categórica.

  3. Del mismo portal de datos donde descargaron el shapefile, o de otra fuente si la tienen, elegir un dataset con registros geo-referenciados. Por ejemplo, las escuelas de la ciudad (o las comisarías, o las propiedades en alquiler, etc) con sus coordenadas. Puede ser en formato csv o shapefile. 3.a. Realizar un “join espacial”, asignando a cada registro geo-referenciado la unidad geográfica (barrio, comuna, etc) que le corresponda y utilizar ggplot() para realizar un mapa coroplético con los límites de los barrios, cuyo color de relleno indique la cantidad encontrada en cada uno.

ACLARACIÓN

Pueden elegir una ciudad de esta planilla, donde ya se han identificado portales con shapefiles. Si eligen otra, por favor agreguen los datos a la planilla así crece de forma colaborativa.