GEOPROCESAMIENTO DE DATOS
Angie Scetta
Mayo 2021
UNIR DATOS A PARTIR DE SUS UBICACIONES ESPACIALES
En la clase anterior (https://rpubs.com/angiescetta/info-espacial) 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().Agregar geometrías extra a nuestro mapa: La magia de
ggplot().
Pero todavía nos queda mucho por recorrer en el mundo de los datos espaciales y hoy 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. Comencemos cargando nuestras librerías:
#install.packages("tidyverse")
library(tidyverse)## Warning: package 'tidyverse' was built under R version 4.0.5## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --## v ggplot2 3.3.3 v purrr 0.3.4
## v tibble 3.1.0 v dplyr 1.0.5
## v tidyr 1.1.3 v stringr 1.4.0
## v readr 1.4.0 v forcats 0.5.1## Warning: package 'ggplot2' was built under R version 4.0.5## Warning: package 'tibble' was built under R version 4.0.5## Warning: package 'tidyr' was built under R version 4.0.5## Warning: package 'readr' was built under R version 4.0.5## Warning: package 'purrr' was built under R version 4.0.5## Warning: package 'dplyr' was built under R version 4.0.5## Warning: package 'stringr' was built under R version 4.0.5## Warning: package 'forcats' was built under R version 4.0.5## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()#install.packages("sf")
library(sf)## Warning: package 'sf' was built under R version 4.0.5## Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1Y el mismo dataset espacial que utilizamos la clase anterior. Recuerden que para poder cargar nuestros datos espaciales en formato shp utilizamos la función st_read() de la siguiente forma:
partidos_amba <- st_read("partidos_amba.shp", stringsAsFactors=TRUE)## Reading layer `partidos_amba' from data source `C:\Users\27356214477\Documents\FEPP-CDDPC-CLASE-1\partidos_amba.shp' using driver `ESRI Shapefile'
## 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 84Repasemos su contenido:
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) :42Tiene 4 columnas que contienen el nombre del partido, la provincia, el área en km2 y la geometría. Para poder visualizar toda esta información plasmada en un mapa vamos a utilizar nuevamente ggplot() y geom_sf() que nos permite sumar capas geográficas (sf).
ggplot(partidos_amba)+
geom_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 hacer algunos ajustes estéticos:
ggplot(partidos_amba)+
geom_sf(fill="gray75", color="white")+
labs(title="Área Metropolitana de Buenos Aires")+
theme_void()
Ahora carguemos nuestra segunda 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:
Cuando utilizamos 2 capas (en este caso
geom_sf()ygeom_point()), tenemos que asignarle a cada una el dataset correspondiente dentro de la capa yggplot()queda “vacío”.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)## although coordinates are longitude/latitude, st_intersects assumes that they are planar
## although coordinates are longitude/latitude, st_intersects assumes that they are planarCabe 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 geometry
## Min. : 6.30 POINT :1247
## 1st Qu.: 14.32 epsg:4326 : 0
## Median : 18.20 +proj=long...: 0
## Mean : 59.05
## 3rd Qu.: 50.30
## Max. :889.49
## NA's :51Como 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())## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_union assumes that they are planarhead(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 15Ahora 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 MAPA BASE
Si bien hasta aquí hemos hecho mapas bastante completos, nunca está de más agregarles un fondo que nos ayude a interpretar mejor el contexto. Para lograr esto utilizaremos la librería ggmap() que cómo ya sabemos, primero debemos instalarla y luego activarla:
#install.packages(ggmap)
library(ggmap)## Warning: package 'ggmap' was built under R version 4.0.5## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.## Please cite ggmap if you use it! See citation("ggmap") for details.Para poder obtener un mapa de fondo con ggmap, primero tengo que delimintar cuál es mi “bounding box”, que hace referencia a los límites de un cuadro en el cual se encuentran todos mis datos. Esto lo haremos de la siguiente forma:
bbox_amba <- as.numeric(st_bbox(partidos_amba))Una vez que ya tenemos el bbox, vamos a usar get_stamenmap() para descargar de internet el mapa. En este caso utilizaremos un mapa de tipo “toner-lite” pero hay otras opciones que pueden verse aquí: http://maps.stamen.com/
mapa_amba <- get_stamenmap(bbox = bbox_amba,
maptype = "toner-lite",
zoom=10)## Source : http://tile.stamen.com/toner-lite/10/343/615.png## Source : http://tile.stamen.com/toner-lite/10/344/615.png## Source : http://tile.stamen.com/toner-lite/10/345/615.png## Source : http://tile.stamen.com/toner-lite/10/346/615.png## Source : http://tile.stamen.com/toner-lite/10/347/615.png## Source : http://tile.stamen.com/toner-lite/10/343/616.png## Source : http://tile.stamen.com/toner-lite/10/344/616.png## Source : http://tile.stamen.com/toner-lite/10/345/616.png## Source : http://tile.stamen.com/toner-lite/10/346/616.png## Source : http://tile.stamen.com/toner-lite/10/347/616.png## Source : http://tile.stamen.com/toner-lite/10/343/617.png## Source : http://tile.stamen.com/toner-lite/10/344/617.png## Source : http://tile.stamen.com/toner-lite/10/345/617.png## Source : http://tile.stamen.com/toner-lite/10/346/617.png## Source : http://tile.stamen.com/toner-lite/10/347/617.png## Source : http://tile.stamen.com/toner-lite/10/343/618.png## Source : http://tile.stamen.com/toner-lite/10/344/618.png## Source : http://tile.stamen.com/toner-lite/10/345/618.png## Source : http://tile.stamen.com/toner-lite/10/346/618.png## Source : http://tile.stamen.com/toner-lite/10/347/618.png## Source : http://tile.stamen.com/toner-lite/10/343/619.png## Source : http://tile.stamen.com/toner-lite/10/344/619.png## Source : http://tile.stamen.com/toner-lite/10/345/619.png## Source : http://tile.stamen.com/toner-lite/10/346/619.png## Source : http://tile.stamen.com/toner-lite/10/347/619.pngCon ggmap() veamos que nos hemos descargado:
ggmap(mapa_amba)
Ahora podemos reutilizar el código de alguno de los mapas anteriores, simplemente reemplazando la primera línea, donde inicializamos un objeto ggplot por una línea que llame a nuestro mapa base de la siguiente forma con ggmap():
ggmap(mapa_amba)+
geom_sf(data=partidos_amba, aes(fill=cantidad/area_km2), alpha=0.75, inherit.aes = FALSE)+
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()## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
Además de cambiar la primer línea de código, nótese que agregamos un inherit.aes=FALSE dentro del geom_sf() para anular la estética predeterminada del objeto ggmap.
Y listo! Ya tenemos un mapa con fondo y con información que generamos nosotros a partir de un cruce espacial.
Ahora les toca practicar a uds!