CIENCIA DE DATOS PARA CIUDADES - CLASE 1
Angie Scetta
FEPP, 2022
ANĆLISIS Y GEOPROCESAMIENTO DE DATOS
En la Ćŗltima clase del MĆ³dulo 1 hicimos una introducciĆ³n a los datos espaciales y vimos los siguientes temas:
Abrir y analizar datos espaciales: El paquete sf.
Unir datos espaciales y datos tradicionales: La funciĆ³n left_join().
Pero todavĆa nos queda mucho 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 pero hay algo que si comparten: su ubicaciĆ³n en el espacio.
A unir datos a partir de sus ubicaciones espaciales 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. Cualquier duda pueden revisar la documentaciĆ³n de la librerĆa.
Pero vamos por partes ya que aĆŗn nos quedan aprender 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).
UNIR DATOS POR SUS UBICACIONES ESPACIALES
A continuaciĆ³n trabajaremos con los datos propiedades en venta de Properati, un portal web de compra, venta y alquiler de inmuebles en toda AmĆ©rica Latina. Estos datos son pĆŗblicos y pueden encontrarlos en el siguiente link.
Pero en este caso, para facilitar la manipulaciĆ³n de la informaciĆ³n durante la clase, usaremos un set de datos (en formato .csv) previamente procesado que contiene datos de propiedades en venta publicadas en AMBA en Junio de 2021. Pueden descargarlo desde este link
RecomendaciĆ³n: Al descargarlo del campus, moverlo de la carpeta āDescargasā a la carpeta llamada ādataā dentro de la carpeta del Proyecto donde estĆ©n trabajando.
Ahora si, Ā”Manos a la obra! Activemos nuestras 2 librerĆas:
#install.packages(tidyverse)
library(tidyverse)
#install.packages(sf)
library(sf)Y carguemos nuestros 2 datsets:
partidos_amba <- st_read("data/partidos_amba.geojson", stringsAsFactors = TRUE)## Reading layer `partidos_amba' from data source
## `C:\Users\27356214477\Documents\UTDT_FEPP_CDDPC\CLASE-01\data\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 84properati <- read.csv("data/amba_properati2021.csv", stringsAsFactors = TRUE)Veamos que informaciĆ³n contiene nuestro dataset de Properati:
names(properati)## [1] "created_on" "rooms" "surface_total" "surface_covered"
## [5] "price" "currency" "property_type" "operation_type"
## [9] "lat" "lon"Las columnas tienen fecha de publicaciĆ³n de la propiedad, cantidad de ambientes, superficie total y superficie cubierta, precio publicado, tipo de moneda (ARS o USD), tipo de propiedad, tipo de operaciĆ³n, y finalmente la ubicaciĆ³n del inmueble con sus coordenadas: latitud y longitud.
Veamos un resumen de nuestros datos:
summary(properati)## created_on rooms surface_total surface_covered
## Min. :202106 Min. : 1.000 Min. : 17.0 Min. : 16.0
## 1st Qu.:202106 1st Qu.: 2.000 1st Qu.: 55.0 1st Qu.: 49.0
## Median :202106 Median : 3.000 Median : 84.0 Median : 71.0
## Mean :202106 Mean : 3.332 Mean : 129.2 Mean :101.7
## 3rd Qu.:202106 3rd Qu.: 4.000 3rd Qu.: 154.0 3rd Qu.:122.0
## Max. :202106 Max. :10.000 Max. :4400.0 Max. :800.0
## price currency property_type operation_type
## Min. : 15000 USD:6111 Casa :1675 Venta:6111
## 1st Qu.: 94000 Departamento:3743
## Median : 145000 PH : 693
## Mean : 210965
## 3rd Qu.: 240000
## Max. :2900000
## lat lon
## Min. :-35.12 Min. :-59.02
## 1st Qu.:-34.64 1st Qu.:-58.54
## Median :-34.60 Median :-58.46
## Mean :-34.61 Mean :-58.47
## 3rd Qu.:-34.56 3rd Qu.:-58.41
## Max. :-34.15 Max. :-57.85Ya podemos ir generando algunos insights:
Los registros corresponden al mes de Junio 2021.
La propiedad de menor superficie cubierta tiene 16 m2 y la de mayor tiene 800 m2.
La mayorĆa de las propiedades publicadas son departamentos (3.743).
Como habrĆ”n notado, 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. AcĆ” vamos a tener que utilizar si o si la ubicaciĆ³n espacial de ambos datasets.
Comencemos utilizando la longitud y latitud en un grƔfico de puntos para ver como lucen en el espacio:
ggplot(properati) +
geom_point(aes(x=lon, y=lat))
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() y geom_point()), tenemos que asignarle a cada una el dataset correspondiente dentro de la capa y ggplot() 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 (properati).
Veamos esto en detalle:
ggplot()+
geom_sf(data=partidos_amba)+
geom_point(data=properati, aes(x=lon, y=lat))
Ahora si! Mucho mĆ”s clara la ubicaciĆ³n de las propiedades publicadas.
Parece que tienen algo en comĆŗn el dataset de partidos y el de Properati, Āæ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, y por suerte el paquete sf tiene una funciĆ³n muy sencilla de utilizar que se llama st_as_sf():
properati_geo <- properati %>%
st_as_sf(coords = c("lon", "lat"), crs = 4326)Revisemos que efectivamente haya ocurrido lo que querĆamos:
class(properati_geo)## [1] "sf" "data.frame"Perfecto, properati_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(properati_geo)## Simple feature collection with 6 features and 8 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -58.90955 ymin: -34.61551 xmax: -58.38517 ymax: -34.45655
## Geodetic CRS: WGS 84
## created_on rooms surface_total surface_covered price currency property_type
## 1 202106 1 26 20 62000 USD PH
## 2 202106 1 60 40 90000 USD Casa
## 3 202106 1 48 36 75000 USD Departamento
## 4 202106 1 40 36 99500 USD Departamento
## 5 202106 1 28 28 68000 USD Departamento
## 6 202106 1 91 53 135000 USD Departamento
## operation_type geometry
## 1 Venta POINT (-58.38517 -34.61551)
## 2 Venta POINT (-58.65276 -34.4987)
## 3 Venta POINT (-58.90955 -34.45655)
## 4 Venta POINT (-58.47684 -34.56799)
## 5 Venta POINT (-58.4227 -34.58355)
## 6 Venta POINT (-58.42886 -34.59237)Tal como esperĆ”bamos, nuestro dataset de Properati 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(properati_geo)+
geom_sf()
Y probemos agregando un poco de color y dejando los partidos por debajo:
ggplot()+
geom_sf(data=partidos_amba, fill="gray75", color="white")+
geom_sf(data=properati_geo, aes(color=property_type), alpha=0.5)+
labs(title="Propiedades publicadas en AMBA")
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 properati_geo sume una nueva columna que indique a que partido pertenece:
properati_geo <- st_join(properati_geo, partidos_amba)Cabe destacar que, en 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(properati_geo)## created_on rooms surface_total surface_covered
## Min. :202106 Min. : 1.000 Min. : 17.0 Min. : 16.0
## 1st Qu.:202106 1st Qu.: 2.000 1st Qu.: 55.0 1st Qu.: 49.0
## Median :202106 Median : 3.000 Median : 84.0 Median : 71.0
## Mean :202106 Mean : 3.332 Mean : 129.2 Mean :101.7
## 3rd Qu.:202106 3rd Qu.: 4.000 3rd Qu.: 154.0 3rd Qu.:122.0
## Max. :202106 Max. :10.000 Max. :4400.0 Max. :800.0
##
## price currency property_type operation_type
## Min. : 15000 USD:6111 Casa :1675 Venta:6111
## 1st Qu.: 94000 Departamento:3743
## Median : 145000 PH : 693
## Mean : 210965
## 3rd Qu.: 240000
## Max. :2900000
##
## nombre provincia area_km2 geometry
## Comuna 14: 503 CABA:3134 Min. : 6.30 POINT :6111
## La Plata : 447 GBA :2966 1st Qu.: 14.73 epsg:4326 : 0
## Comuna 13: 418 NA's: 11 Median : 22.14 +proj=long...: 0
## Comuna 12: 318 Mean :141.85
## Tigre : 290 3rd Qu.:120.68
## (Other) :4124 Max. :889.49
## NA's : 11 NA's :11A cada registro se le uniĆ³ una columna del nombre del partido con el que se solapa. Pero acĆ” aparece algo nuevo, algunos valores nulos o NAĀ“s que hacen referencia a aquellas propiedades ubicadas 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=properati_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:
properati_geo <- properati_geo %>%
filter(!is.na(nombre))Visualicemos nuevamente:
ggplot()+
geom_sf(data=partidos_amba)+
geom_sf(data=properati_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 propiedades publicadas que tienen o segĆŗn el valor del m2 promedio? Nuevamente manipulamos el set de datos con group_by() y luego utilizamos left_join():
properati_geo <- properati_geo %>%
group_by(nombre) %>%
summarise(cantidad=n(),
valor_m2=mean(price/surface_covered))head(properati_geo)## Simple feature collection with 6 features and 3 fields
## Geometry type: GEOMETRY
## Dimension: XY
## Bounding box: xmin: -58.52936 ymin: -34.93481 xmax: -57.8822 ymax: -34.58482
## Geodetic CRS: WGS 84
## # A tibble: 6 x 4
## nombre cantidad valor_m2 geometry
## <fct> <int> <dbl> <GEOMETRY [Ā°]>
## 1 Almirante Brown 22 1341. MULTIPOINT ((-58.35463 -34.81037), (-58.335~
## 2 Avellaneda 73 1310. MULTIPOINT ((-58.30226 -34.69496), (-58.301~
## 3 Berazategui 30 1546. MULTIPOINT ((-58.10421 -34.82882), (-58.170~
## 4 Berisso 1 1200 POINT (-57.8822 -34.93481)
## 5 Comuna 1 267 3122. MULTIPOINT ((-58.37084 -34.62929), (-58.370~
## 6 Comuna 10 95 1990. MULTIPOINT ((-58.52936 -34.61466), (-58.526~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 properati_geo a un dataset tradicional, quitĆ”ndole la geometrĆa con st_set_geometry():
properati_geo <- properati_geo %>%
st_set_geometry(NULL)head(properati_geo)## # A tibble: 6 x 3
## nombre cantidad valor_m2
## <fct> <int> <dbl>
## 1 Almirante Brown 22 1341.
## 2 Avellaneda 73 1310.
## 3 Berazategui 30 1546.
## 4 Berisso 1 1200
## 5 Comuna 1 267 3122.
## 6 Comuna 10 95 1990.Ahora si, solo tiene 3 columnas con el nombre del partido, la cantidad de propiedades y el valor promedio del m2, pero ya no tiene geometrĆa. Procedamos a realizar la uniĆ³n:
partidos_amba <- left_join(partidos_amba, properati_geo, by="nombre")Probemos mapear las cantidades:
ggplot()+
geom_sf(data=partidos_amba, aes(fill=cantidad))
Ya tenemos un mapa de cantidad de propiedades en venta por partido, y podemos ver como la mayor cantidad de concentra en la Comuna 14 (CABA), seguida por La Plata (GBA).
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), color="black")+
labs(title = "Propiedades en venta por km2",
subtitle = "Partidos de AMBA",
fill = "Cantidad/km2",
caption= "Fuente: Properati - Junio 2021") +
scale_fill_distiller(palette = "YlOrRd", direction = 1) +
theme_light()
Y ahora veamos un mapa coroplƩtico de valor de venta del m2:
ggplot()+
geom_sf(data=partidos_amba, aes(fill=valor_m2), color=NA)+
labs(title = "Propiedades en venta - Valor promedio del m2",
subtitle = "Partidos de AMBA",
fill = "Valor m2",
caption= "Fuente: Properati - Junio 2021") +
scale_fill_distiller(palette = "YlOrRd", direction = 1) +
theme_light()
AGREGAR MAPA BASE
Si bien hasta aquĆ hemos hecho mapas bastante completos, nunca estĆ” de mĆ”s agregar un fondo que nos ayude a interpretar mejor el contexto de los datos que estamos analizando.
Para lograr esto utilizaremos la librerĆa ggmap que nos permitirĆ” aƱadir una capa que funcione como base de nuestros mapas. Cabe destacar que, este paquete utiliza recursos de diferentes APIs (Google, Stamenmap, etc), y por lo tanto es necesario tener conexiĆ³n a internet al momento de utilizarlo. Cualquier duda pueden consultar la documentaciĆ³n de ggmap.
Pero antes de seguir debemos instalarla y activarla:
#install.packages(ggmap)
library(ggmap)Una vez activada nuestra librerĆa, el siguiente paso es determinar cuĆ”l es la ābounding boxā o caja de coordenadas de mi dataset espacial. Es decir, cuĆ”les son los lĆmites geogrĆ”ficos del cuadro que contiene todos los datos.
Para esto, utilizaremos la funciĆ³n st_bbox() sobre el dataset espacial que tiene mayor cobertura en el especio: los partidos del AMBA.
bbox_amba <- st_bbox(partidos_amba)bbox_amba## xmin ymin xmax ymax
## -59.33920 -35.23893 -57.70946 -34.23007Bien, ya tengo los 4 valores que delimitan mi bounding box, pero el formato es de tipo ābboxā y para utilizar ggmap necesitamos que sea de tipo numĆ©rico, asĆ que modifiquemoslo:
bbox_amba <- as.numeric(bbox_amba)Ahora vamos a usar get_stamenmap() para descargar de internet el āmapa baseā. En este caso utilizaremos un mapa de tipo ātoner-liteā pero hay otras opciones que pueden verse en este link
mapa_amba <- get_stamenmap(bbox = bbox_amba,
maptype = "toner-lite",
zoom=10)Con 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 iniciamos un objeto ggplot por una lĆnea que llame a nuestro mapa base con ggmap() y luego agregando el parĆ”metro āinherit.aes=FALSEā dentro del geom_sf() para anular la estĆ©tica predeterminada del objeto:
ggmap(mapa_amba)+
geom_sf(data=partidos_amba, aes(fill=valor_m2), color=NA, inherit.aes=FALSE, alpha=0.8)+
labs(title = "Propiedades en venta - Valor promedio del m2",
subtitle = "Partidos de AMBA",
fill = "Valor m2",
caption= "Fuente: Properati - Junio 2021") +
scale_fill_distiller(palette = "YlOrRd", direction = 1) +
theme_light()
Listo! Ya tenemos un mapa con fondo y con informaciĆ³n que generamos nosotros a partir de un cruce de 2 bases de datos.
Ahora les toca practicar a uds!