Introducción

Me llamo Andy y mi dataset público preferido es el Relevamiento de Usos de Suelo.

Sin más demoras, arranquemos…

Desarrollo

Cargo los paquetes necesarios.

library(sf)
library(dplyr)
library(here)
library(data.table)
library(ggplot2)
library(ggpointdensity)
library(viridis)
library(ggmap)

Importo el RUS 2017. Pueden descomentar el comando fread si no lo tienen descargado.

rus <- read.csv(here("data", "relevamiento-usos-del-suelo-2017.csv"),
                         encoding = "UTF-8")

# rus <- fread("http://cdn.buenosaires.gob.ar/datosabiertos/datasets/relevamiento-usos-del-suelo/relevamiento-usos-del-suelo-2017.csv",
#                   encoding = 'UTF-8',
#                  sep = ",",
#                  header = TRUE,
#                  showProgress = TRUE,
#                  data.table = FALSE)

nrow(rus)
## [1] 561722

Como vemos, tiene más de 550.000 filas.

Miro las primeras filas.

rus %>% head()
##           X         Y           SMP          CALLE  NUM TIPO1_16
## 1 -58.45563 -34.59689 047-040A-001C     MONTENEGRO  124        E
## 2 -58.45570 -34.59694  047-040A-002 NUEVA ZELANDIA 4820       EP
## 3 -58.45575 -34.59701  047-040A-003 NUEVA ZELANDIA 4830        E
## 4 -58.45581 -34.59707  047-040A-004 NUEVA ZELANDIA 4840      GAP
## 5 -58.45588 -34.59713  047-040A-005 NUEVA ZELANDIA 4850      GAP
## 6 -58.45588 -34.59713  047-040A-005 NUEVA ZELANDIA 4850        E
##         TIPO2_16 PISOS_16      NOMBRE OBSERVACIO    BARRIO COMUNA X5_DIG
## 1       VIVIENDA        1                        CHACARITA     15      3
## 2     METALURGIA        1 RODOCA S.A.            CHACARITA     15  28900
## 3           CASA        1                E X GAP CHACARITA     15      4
## 4 GARAGE PRIVADO        1                        CHACARITA     15  63312
## 5 GARAGE PRIVADO        1                        CHACARITA     15  63312
## 6           CASA        2                        CHACARITA     15      4
##   X4_DIG X3_DIG X2_DIG
## 1      0      0      0
## 2   2890    289     28
## 3      0      0      0
## 4   6331    633     63
## 5   6331    633     63
## 6      0      0      0
##                                                                         RAMA
## 1                                                                  EDIFICIOS
## 2 FABRICACIÓN DE PRODUCTOS ELABORADOS DE METAL, EXCEPTO MAQUINARIA  Y EQUIPO
## 3                                                                  EDIFICIOS
## 4                                    SERVICIOS COMPLEMENTARIOS AL TRANSPORTE
## 5                                    SERVICIOS COMPLEMENTARIOS AL TRANSPORTE
## 6                                                                  EDIFICIOS
##                                             SUBRAMA
## 1                                          VIVIENDA
## 2  FABRICACIÓN DE PRODUCTOS ELABORADOS DE METAL NCP
## 3                                              CASA
## 4 SERVICIOS COMPLEMENTARIOS AL TRANSPORTE TERRESTRE
## 5 SERVICIOS COMPLEMENTARIOS AL TRANSPORTE TERRESTRE
## 6                                              CASA
##                                 SSRAMA
## 1                                     
## 2                                     
## 3                                     
## 4 ESTACIONAMIENTOS, COCHERAS Y GARAJES
## 5 ESTACIONAMIENTOS, COCHERAS Y GARAJES
## 6

Veo los principales rubros.

rus %>%
    select(TIPO1_16) %>%
    table() %>%
    as.data.frame %>%
    arrange(desc(Freq))
##       .   Freq
## 1     E 270677
## 2     L 116862
## 3   GAP 115083
## 4   EDU  15389
## 5    LG  13186
## 6    EP   7973
## 7  LOTE   5209
## 8     G   4788
## 9    SU   3477
## 10    U   3274
## 11   GA   2818
## 12    P   2646
## 13   ES    340
rus %>%
    select(TIPO2_16) %>%
    table() %>%
    as.data.frame() %>%
    arrange(desc(Freq)) %>%
    top_n(20)
## Selecting by Freq
##                                 .   Freq
## 1                        VIVIENDA 149060
## 2                  GARAGE PRIVADO 115083
## 3                            CASA  97541
## 4                   LOCAL CERRADO  27937
## 5             VU - VIVIENDA UNICA  15738
## 6             EDIFICIO PRODUCTIVO   6622
## 7                        DEPOSITO   6544
## 8                        OFICINAS   6403
## 9                      I FEMENINA   6098
## 10 TALLER MECANICO DE AUTOMOTORES   3523
## 11                     PELUQUERIA   3187
## 12                              I   2951
## 13                         GARAGE   2783
## 14                   SUPERMERCADO   2628
## 15                   INMOBILIARIA   2387
## 16                           LOTE   2311
## 17                         KIOSCO   2307
## 18                            EDU   2176
## 19                     ABANDONADO   2146
## 20                        EN OBRA   2019
table(rus$TIPO1_16, rus$TIPO2_16) %>%
    as.data.frame() %>%
    filter(Freq > 0) %>%
    arrange(desc(Freq)) %>%
    top_n(20)
## Selecting by Freq
##    Var1                           Var2   Freq
## 1     E                       VIVIENDA 149060
## 2   GAP                 GARAGE PRIVADO 115083
## 3     E                           CASA  97541
## 4     L                  LOCAL CERRADO  25903
## 5     E            VU - VIVIENDA UNICA  15738
## 6    EP            EDIFICIO PRODUCTIVO   6622
## 7     L                     I FEMENINA   4139
## 8     L TALLER MECANICO DE AUTOMOTORES   3050
## 9     L                     PELUQUERIA   3029
## 10    G                       DEPOSITO   2895
## 11   GA                         GARAGE   2783
## 12    L                   SUPERMERCADO   2621
## 13    L                       OFICINAS   2582
## 14    E                       OFICINAS   2466
## 15    L                   INMOBILIARIA   2339
## 16 LOTE                           LOTE   2310
## 17    L                         KIOSCO   2254
## 18  EDU                            EDU   2176
## 19    E                     ABANDONADO   2120
## 20  EDU                       DEPOSITO   2096

Acá vemos diversidad de categorías, algunas más intuitivas que otras. ¿L es local? ¿I FEMENINA es indumentaria?

Veamos las principales categorías del tipo L.

rus %>%
    filter(TIPO1_16 == "L") %>%
    select(TIPO2_16) %>%
    table() %>%
    as.data.frame() %>%
    filter(Freq > 0) %>%
    arrange(desc(Freq)) %>%
    top_n(20)
## Selecting by Freq
##                                 .  Freq
## 1                   LOCAL CERRADO 25903
## 2                      I FEMENINA  4139
## 3  TALLER MECANICO DE AUTOMOTORES  3050
## 4                      PELUQUERIA  3029
## 5                    SUPERMERCADO  2621
## 6                        OFICINAS  2582
## 7                    INMOBILIARIA  2339
## 8                          KIOSCO  2254
## 9                      MAXIKIOSCO  1965
## 10                    RESTAURANTE  1670
## 11                              I  1650
## 12                     VERDULERIA  1411
## 13                            BAR  1364
## 14                        CALZADO  1332
## 15               LAVADERO DE ROPA  1227
## 16                       DEPOSITO  1224
## 17                        LOTERIA  1179
## 18                     FERRETERIA  1084
## 19                        ALMACEN  1065
## 20                MUEBLES (VENTA)  1065

Me quedo entonces con las principales 20 categorías de locales, quitando los locales cerrados, los depósitos y las oficinas.

categorias <- rus %>%
    filter(TIPO1_16 == "L" & !(TIPO2_16 %in% c("LOCAL CERRADO", "DEPOSITO", "OFICINAS"))) %>%
    select(TIPO2_16) %>%
    table() %>%
    as.data.frame() %>%
    filter(Freq > 0) %>%
    arrange(desc(Freq)) %>%
    top_n(20)
## Selecting by Freq
names(categorias)[1] <- "categoria"

locales <- rus %>%
    filter(TIPO1_16 == "L" & TIPO2_16 %in% categorias$categoria)

locales %>% head()
##           X         Y          SMP              CALLE  NUM TIPO1_16
## 1 -58.45545 -34.59703 047-040A-019         MONTENEGRO  110        L
## 2 -58.45545 -34.59703 047-040A-019         MONTENEGRO  112        L
## 3 -58.45620 -34.59678 047-040B-004 ARENAL, CONCEPCION 1826        L
## 4 -58.45645 -34.59701 047-040B-008 ARENAL, CONCEPCION 4860        L
## 5 -58.45629 -34.59713 047-040B-017     NUEVA ZELANDIA 4877        L
## 6 -58.45610 -34.59696 047-040B-020     NUEVA ZELANDIA 4851        L
##                         TIPO2_16 PISOS_16 NOMBRE OBSERVACIO    BARRIO
## 1 TALLER MECANICO DE AUTOMOTORES        1                   CHACARITA
## 2 TALLER MECANICO DE AUTOMOTORES        1                   CHACARITA
## 3 TALLER MECANICO DE AUTOMOTORES        1                   CHACARITA
## 4                     MAXIKIOSCO        1                   CHACARITA
## 5                     FERRETERIA        2                   CHACARITA
## 6                         KIOSCO        1                   CHACARITA
##   COMUNA X5_DIG X4_DIG X3_DIG X2_DIG                                  RAMA
## 1     15  50209   5020    502     50 VENTA Y REP. DE VEHICULOS AUTOMOTORES
## 2     15  50209   5020    502     50 VENTA Y REP. DE VEHICULOS AUTOMOTORES
## 3     15  50209   5020    502     50 VENTA Y REP. DE VEHICULOS AUTOMOTORES
## 4     15  52119   5211    521     52                              COMERCIO
## 5     15  52363   5236    523     52                              COMERCIO
## 6     15  52119   5211    521     52                              COMERCIO
##                                      SUBRAMA
## 1                 MANTENIMIENTO Y REPARACION
## 2                 MANTENIMIENTO Y REPARACION
## 3                 MANTENIMIENTO Y REPARACION
## 4 VENTA AL POR MENOR EXEPTO LA ESPECIALIZADA
## 5       VENTA AL POR MENOR ESPECIALIZADO NCP
## 6 VENTA AL POR MENOR EXEPTO LA ESPECIALIZADA
##                                                        SSRAMA
## 1                                                            
## 2                                                            
## 3                                                            
## 4                                                  MAXIKIOSCO
## 5 MATERIALES DE CONSTRUCCION, FERRETERIA, PINTURERIA, VIDRIOS
## 6                                                      KIOSCO

Dato interesante: aún existen las loterías.

Veo con cuántos locales me estoy quedando.

nrow(locales)
## [1] 35436

¡Hora de mapear!

Veamos la distribución espacial de estas categorías.

Primero debo pasar el objeto locales a tipo geográfico. Y también me traigo los barrios para tener una referencia en el espacio.

(“¿te gustan los mapas?” es una gran pregunta de primera cita que podría tomar de lxs geografxs)

locales <- locales %>%
    st_as_sf(coords = c("X", "Y"),
             crs = 4326,
             remove = FALSE
             )

barrios <- st_read("http://cdn.buenosaires.gob.ar/datosabiertos/datasets/barrios/barrios.geojson")
## Reading layer `barrios_badata' from data source `http://cdn.buenosaires.gob.ar/datosabiertos/datasets/barrios/barrios.geojson' using driver `GeoJSON'
## Simple feature collection with 48 features and 4 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -58.53152 ymin: -34.70529 xmax: -58.33515 ymax: -34.52649
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs

El parámetro remove deja las 2 columnas de latitud, longitud. Presunto que serán usadas en el futuro.

Con esto ya puedo empezar a hacer algunos mapas básicos.

ggplot() +
    geom_sf(data = barrios, aes()) +
    geom_sf(data = locales, aes())

Marche al Instagram de los gráficos que son horribles en el primer intento.

Va de nuevo.

ggplot() +
    geom_sf(data = barrios, aes()) +
    geom_sf(data = locales,
            aes(color = TIPO2_16),
            size = 1,
            alpha = 0.5)

Mejoró rotundamente, pero existe mucha superposición, y los puntos de arriba claramente tapan a los de abajo.

Veamos creando un plano por categoría de comercio.

ggplot() +
    geom_sf(data = barrios, aes()) +
    geom_sf(data = locales,
            aes(color = TIPO2_16),
            size = 1,
            alpha = 0.5) +
    facet_wrap(~TIPO2_16)

Mejoró, pero la leyenda es un poco molesta y los mapas se ven muy tapados, por lo que no puedo apreciar la densidad de comercios.

ggplot() +
    geom_sf(data = barrios, aes()) +
    geom_sf(data = locales,
            aes(color = TIPO2_16),
            size = 1,
            alpha = 0.2) +
    facet_wrap(~TIPO2_16) +
    theme(legend.position = "none")

Si se agranda la imagen ya pueden verse algunos patrones, como el predominio de ciertas categorías en el centro de la ciudad, Recoleta y Palermo. También algunas categorías poseen una distribución más uniforme a lo largo de la ciduad, como kioskos, supermercados y almacenes.

Sin embargo, aún no percibo con claridad la densidad, porque la opacidad luego de cierto nivel de superposición es engañosa.

Veamos cómo sería si usamos geom_point en lugar de geom_sf.

Y la paleta viridis, porque puedo. <3

ggplot() +
    geom_sf(data = barrios, aes()) +
    geom_pointdensity(data = locales,
                      aes(x = X, y = Y),
                      size = 1,
                      adjust = 0.1) +
    scale_color_viridis() +
    facet_wrap(~TIPO2_16)

# vuelvo a dejar la leyenda
# theme(legend.position = "none")

El paquete ggpointdensity tiene la virtud de que uno plotea los puntos, y se colorean por la densidad.

Eso permite ver rápidamente qué puntos de la ciudad tienen mayor concentración de cada rubro.

El parámetro adjust permite modificar qué entorno se considera para obtener la densidad. Al reducirlo pueden aparece otras modas en el espacio que habían sido suavizadas, como en los rubros de indumentaria femenina o talleres mecánicos.

Un punto interesante es que para todos los plots utiliza la misma escala de colores, por lo que cada facetado es comparable con los demás.

Otro paquete que ofrece algo interesante es ggmap. Pero primero preciso un de fondo de la bella CABA.

bbox <- barrios %>%
    st_bbox() %>%
    as.numeric()
# si le paso a get_stamenmap() el tipo bbox no funciona -_-

caba <- get_stamenmap(bbox = bbox,
                      maptype = "toner-lite",
                      zoom = 12)

Logré bajar las tiles, a pesar de que habían pasado cosas con una función que pedía API. Alguien tiene la gorra puesta.

Pruebo stats_density_2d, y saco viridis porque muchos valores bajos me oscurece demasiado.

Quito el geom_sf pues pruebo ggmaps, y utilizo fill en lugar de color para la densidad. Y agrego labs para darle contexto.

ggmap(caba) +
    stat_density_2d(data = locales,
                  aes(x = X, y = Y, fill = stat(level)),
                  alpha = 0.6,
                  geom = "polygon",
                  size = 0.5) +
    scale_fill_distiller(palette = "RdYlGn") +
    facet_wrap(~TIPO2_16) +
    labs(title = "Densidad de comercios por rubro",
       subtitle = "Ciudad Autónoma de Buenos Aires",
       x = "",
       y = "",
       caption = "Fuente: Elaboración propia con datos públicos y abiertos.",
       fill = "Nivel")

Estaba feucho, y logré cambiar la paleta a algo más alegre.

Intento algo similar, con la paleta spectral.

ggplot() +
    geom_sf(data = barrios, aes(), alpha = 0.5) +
    geom_density_2d(data = locales,
                  aes(x = X, y = Y, color = stat(level)),
                  alpha = 0.8,
                  size = 1) +
    scale_color_distiller(palette = "Spectral") +
    facet_wrap(~TIPO2_16) +
    labs(title = "Densidad de comercios por rubro",
       subtitle = "Ciudad Autónoma de Buenos Aires",
       x = "",
       y = "",
       caption = "Fuente: Elaboración propia con datos públicos y abiertos.",
       fill = "Nivel")

También me interesa entender la interacción entre los distintos rubros de locales. Pero, suspendo por hoy, queda para el próximo episodio.

No sin antes ser agradecido.