TP FINAL: Instrumentos de Análisis Territorial II

El trabajo consiste en leer, explorar y estudiar la relación entre dos o más variables de datos que tengan una representación espacial como puntos o polígonos. El trabajo deberá tener una hipótesis o pregunta a responderse, y estructurarse entorno a ella.

Más allá de los datos específicos que se estudien, en el trabajo deben utilizarse obligatoriamente todas las siguientes herramientas vistas en clase: 1) Unión espacial (spatial join usando la función st_join()) 2) Elaboración de mapas (usando cualquier paquete de R o GeoDa) 3) Al menos una de las siguientes tres: geocodificación, areal weighted interpolation o cálculo de distancias (st_distance() ) 4) Análisis de autocorrelación espacial (test de Moran). 5) Estimación de regresión lineal 6) Análisis de necesidad de modificar la regresión lineal por alguno de los modelos espaciales vistos en clase 7) Se valorará incluir en el reporte, si fuera correspondiente para el análisis, el uso de otras herramientas vistas en clase (por ejemplo, algunas de las herramientas de aprendizaje automático o de clustering)

Para el presente trabajo se analizará la ciudad de Barcelona y la cobertura del transporte público La pregunta que el trabajo pretende responder es ¿Cómo es la accesibilidad al transporte público y cómo influye en el valor del suelo?

Para tal fin, cargamos las siguientes librerías:

library(plyr)
library(tidyverse)

## ── Attaching packages ────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──

## ✓ ggplot2 3.3.2     ✓ purrr   0.3.4
## ✓ tibble  3.0.3     ✓ dplyr   1.0.0
## ✓ tidyr   1.1.0     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.5.0

## ── Conflicts ───────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::arrange()   masks plyr::arrange()
## x purrr::compact()   masks plyr::compact()
## x dplyr::count()     masks plyr::count()
## x dplyr::failwith()  masks plyr::failwith()
## x dplyr::filter()    masks stats::filter()
## x dplyr::id()        masks plyr::id()
## x dplyr::lag()       masks stats::lag()
## x dplyr::mutate()    masks plyr::mutate()
## x dplyr::rename()    masks plyr::rename()
## x dplyr::summarise() masks plyr::summarise()
## x dplyr::summarize() masks plyr::summarize()

library(sf)

## Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0

Cargamos un dataset espacial con los límites geográficos de la ciudad de Barcelona, que provee la página del Ayuntamiento en su sección “open data”:

Barcelona_limites <- read_sf("BCN_UNITATS_ADM_POLIGONS.json")

st_crs(Barcelona_limites)

## Coordinate Reference System:
##   User input: 25831 
##   wkt:
## PROJCS["ETRS89 / UTM zone 31N",
##     GEOGCS["ETRS89",
##         DATUM["European_Terrestrial_Reference_System_1989",
##             SPHEROID["GRS 1980",6378137,298.257222101,
##                 AUTHORITY["EPSG","7019"]],
##             TOWGS84[0,0,0,0,0,0,0],
##             AUTHORITY["EPSG","6258"]],
##         PRIMEM["Greenwich",0,
##             AUTHORITY["EPSG","8901"]],
##         UNIT["degree",0.0174532925199433,
##             AUTHORITY["EPSG","9122"]],
##         AUTHORITY["EPSG","4258"]],
##     PROJECTION["Transverse_Mercator"],
##     PARAMETER["latitude_of_origin",0],
##     PARAMETER["central_meridian",3],
##     PARAMETER["scale_factor",0.9996],
##     PARAMETER["false_easting",500000],
##     PARAMETER["false_northing",0],
##     UNIT["metre",1,
##         AUTHORITY["EPSG","9001"]],
##     AXIS["Easting",EAST],
##     AXIS["Northing",NORTH],
##     AUTHORITY["EPSG","25831"]]

Vemos que el dataset está proyectado en EPSG = 25831. Vamos a cambiarlo a Merkator, ya que trabajaremos más adelante con datasets con Coordinated Reference System (CRS) en Merkator (4326).

Barcelona_limites <- st_transform(Barcelona_limites, crs = 4326)

Nos quedamos con los límites de los distritos, algo similar a las comunas en Buenos Aires:

Barcelona_distritos <- Barcelona_limites %>% 
  filter(!(DISTRICTE =="-")) %>% 
  group_by(DISTRICTE) %>% 
  summarise()

## `summarise()` ungrouping output (override with `.groups` argument)

A su vez, nos falta la información de los nombres de los barrios y distritos. Cargamos la información faltante que se encuentra en otro dataset y le hacemos las modificaciones necesarias para luego unirlo al dataset que contiene la información geográfica:

nombres_limites <- read.csv("BCN_districtes_i_barris.csv")

nombres_limites <- nombres_limites %>% 
    mutate(CODI_DISTRICTE = as.character(CODI_DISTRICTE))

nombres_limites$CODI_DISTRICTE <- revalue(nombres_limites$CODI_DISTRICTE, c("1"= "01", "2"="02", "3"="03", "4"="04", "5"="05", "6"="06", "7"="07", "8"="08", "9"="09"))

Elaboración de mapas

Graficamos los límites de distritos:

Barcelona_distritos <- left_join(Barcelona_distritos, nombres_limites, by = c("DISTRICTE" = "CODI_DISTRICTE")) %>% 
  select(-CODI_BARRI, -NOM_BARRI) %>% 
  mutate(Ref_distrito = paste0(DISTRICTE,": ", NOM_DISTRICTE))

ggplot() +
        geom_sf(data = Barcelona_distritos, aes(fill = Ref_distrito), color = NA) +
        geom_sf_text(data = Barcelona_distritos, aes(label = DISTRICTE), size = 2, colour = "white") +
        scale_fill_brewer(palette="PRGn") +
theme_void() +
         labs(
           title = "Límites de Distritos",
           subtitle = "Ciudad de Barcelona",
           fill = "Distritos",
           caption = "Fuente: Ayuntamiento de Barcelona | 2020")

## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data

Del mismo portal de datos elegimos un dataset con registros geo-referenciados de las ubicaciones de las Estaciones / Paradas del transporte público en la ciudad de Barcelona:

transporte <- read.csv("BCN_TRANSPORTS.csv")

Superponemos estas ubicaciones sobre el mapa de distritos que construimos anteriormente:

ggplot() + 
    geom_sf(data = Barcelona_distritos, aes(fill = Ref_distrito), color = NA) +
    geom_point(data = transporte, aes(x=LONGITUD, y=LATITUD), color = "black") +
    scale_fill_brewer(palette="PRGn") +      
    theme_void() +
         labs(title = "Transporte Público en Barcelona",
              subtitle = "Estaciones de metro, tren y tranvía",
              fill = "Distritos",
              caption = "Fuente: Ayuntamiento de Barcelona | 2020")

Analizamos que categorías de Transporte Público tiene el dataset:

ggplot() + 
    geom_sf(data = Barcelona_distritos, fill = "snow3", color = "white") +
    geom_point(data = transporte, aes(x=LONGITUD, y=LATITUD, color = NOM_CAPA), alpha=.9) +
    guides(color = guide_legend(override.aes = list(size = 3, alpha = 1))) + 
    scale_color_viridis_d(direction = -1) +  
    theme_void() +
         labs(title = "Transporte Público en Barcelona",
              subtitle = "Estaciones de metro, tren y tranvía",
              color = "Líneas",
              caption = "Fuente: Ayuntamiento de Barcelona | 2020")

A su vez, podemos analizar la distribución espacial de los datos a partir de un mapa de densidad que muestre donde se concentra la mayor cantidad de observaciones. Primero obtenemos la bounding box de nuestros datos, y luego los pasamos a get_stamenmap() para hacer el mapa:

library(ggmap)

## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.

## Please cite ggmap if you use it! See citation("ggmap") for details.

bbox <- make_bbox(lon = transporte$LONGITUD, lat = transporte$LATITUD)
MAPA <- get_stamenmap(bbox = bbox, maptype = "toner-lite", zoom = 12)

## Source : http://tile.stamen.com/toner-lite/12/2070/1528.png

## Source : http://tile.stamen.com/toner-lite/12/2071/1528.png

## Source : http://tile.stamen.com/toner-lite/12/2072/1528.png

## Source : http://tile.stamen.com/toner-lite/12/2073/1528.png

## Source : http://tile.stamen.com/toner-lite/12/2070/1529.png

## Source : http://tile.stamen.com/toner-lite/12/2071/1529.png

## Source : http://tile.stamen.com/toner-lite/12/2072/1529.png

## Source : http://tile.stamen.com/toner-lite/12/2073/1529.png

## Source : http://tile.stamen.com/toner-lite/12/2070/1530.png

## Source : http://tile.stamen.com/toner-lite/12/2071/1530.png

## Source : http://tile.stamen.com/toner-lite/12/2072/1530.png

## Source : http://tile.stamen.com/toner-lite/12/2073/1530.png

## Source : http://tile.stamen.com/toner-lite/12/2070/1531.png

## Source : http://tile.stamen.com/toner-lite/12/2071/1531.png

## Source : http://tile.stamen.com/toner-lite/12/2072/1531.png

## Source : http://tile.stamen.com/toner-lite/12/2073/1531.png

ggmap(MAPA) +
    geom_bin2d(data = transporte, 
               aes(x = LONGITUD, y = LATITUD), bins = 30, alpha=.9) +
   scale_fill_viridis_c(direction=-1) +
   theme_void()+
   labs(title = "Densidad de Transporte público",
         subtitle = "Ciudad de Barcelona",
         x = "Longitud",
         y = "Latitud",
         fill = "Cantidad",
         caption = "Fuente: Ayuntamiento de Barcelona")

Join spacial

A continuación realizaremos un join espacial, para unir las ubicaciones de Estaciones de Transporte Público con el dataset que contiene los límites de distritos de Barcelona. Para eso, primero convertimos el dataset con las ubicaciones a un formato espacial:

transporte_esp <- transporte %>%
    st_as_sf(coords = c("LONGITUD", "LATITUD"), crs = 4326) %>% select(EQUIPAMENT, NOM_CAPA)

Barcelona_transporte <- st_join(Barcelona_distritos, transporte_esp) %>% 
  group_by(EQUIPAMENT, NOM_CAPA, DISTRICTE) %>% 
  summarise()

## although coordinates are longitude/latitude, st_intersects assumes that they are planar

## `summarise()` regrouping output by 'EQUIPAMENT', 'NOM_CAPA' (override with `.groups` argument)

Realizamos un conteo de la cantidad de Estaciones y Paradas de Transporte Público por Distrito:

Barcelona_transporte  <- Barcelona_transporte  %>% 
      group_by(DISTRICTE) %>%
      summarise(Frecuencia = n())

## `summarise()` ungrouping output (override with `.groups` argument)

Con un gráfico de barras vemos la cantidad de Estaciones y Paradas de Transporte Público por Distrito:

ggplot(Barcelona_transporte) +
  geom_bar(aes( x = DISTRICTE, weight = Frecuencia, fill=Frecuencia), alpha=.95) +
  scale_fill_viridis_c() +
  theme_minimal() +
  labs(title = "Transporte Público por distrito en Barcelona",
       subtitle = "Estaciones de metro, tren y tranvía",
       fill = "Escala",
       caption= "Fuente: Ayuntamiento de Barcelona | 2020",
       y="Cantidad de Estaciones de Transporte Público",
       x="Distritos")

Y mejor aún en un mapa con los límites de los distritos, cuyo color de relleno indica el conteo de Estaciones y Paradas de Transporte Público en cada uno:

ggplot() + 
  geom_sf(data = Barcelona_transporte, aes(fill=Frecuencia), color = NA, alpha =.95) + 
  geom_sf_text(data = Barcelona_transporte, aes(label = Frecuencia), size = 2, colour = "white") +
  scale_fill_viridis_c() +
  theme_void() +
  labs(title = "Accesibilidad al Transporte por distrito en Barcelona",
         subtitle = "Estaciones de metro, tren y tranvía",
         fill = "Escala",
         caption= "Fuente: Ayuntamiento de Barcelona | 2020")

## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data

Al mismo gráfico le superponemos las ubicaciones de las estaciones y paradas de Transporte Público.

ggplot() + 
  geom_sf(data = Barcelona_transporte, aes(fill=Frecuencia), color = NA, alpha =.95) + 
  geom_sf_text(data = Barcelona_transporte, aes(label = Frecuencia), size = 2, colour = "white") +
  geom_sf(data=transporte_esp, color = "black") +
  scale_fill_viridis_c() +
  theme_void() +
  labs(title = "Accesibilidad al Transporte por distrito en Barcelona",
         subtitle = "Estaciones de metro, tren y tranvía",
         fill = "Escala",
         caption= "Fuente: Ayuntamiento de Barcelona | 2020")

## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data

Como podemos ver, hay observaciones fuera de los límites de la ciudad de Barcelona, nos quedaremos sólo con los que estan dentro de los límites de la ciudad y volvemos a graficar:

Transporte_BCN <- st_intersection(Barcelona_distritos, transporte_esp)

## although coordinates are longitude/latitude, st_intersection assumes that they are planar

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

ggplot() + 
  geom_sf(data = Barcelona_transporte, aes(fill=Frecuencia), color = NA, alpha =.95) + 
  geom_sf_text(data = Barcelona_transporte, aes(label = Frecuencia), size = 2, colour = "white") +
  geom_sf(data=Transporte_BCN, color = "black") +
  scale_fill_viridis_c() +
  theme_void() +
  labs(title = "Accesibilidad al Transporte por distrito en Barcelona",
         subtitle = "Estaciones de metro, tren y tranvía",
         fill = "Escala",
         caption= "Fuente: Ayuntamiento de Barcelona | 2020")

## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data

A su vez, con la ayuda de leaflet podemos agregar un popup que despliegue el nombre de cada estación:

library(leaflet)

icons <- awesomeIcons(icon = 'ios-close', iconColor = 'white', library = 'ion', markerColor = "darkpurple")

leaflet(transporte_esp) %>%      
  addProviderTiles(providers$CartoDB.Positron) %>%
  addAwesomeMarkers(icon=icons, popup = ~EQUIPAMENT)

Aereal weight

Medimos la cobertura de las paradas de transporte público, tanto metro como tren y tranvía, considerando 6 minutos de caminata. Para ello, cruzamos estos datos con las manzanas e identificamos aquellas manzanas que no se encuentran en este espacio de cobertura.

library(hereR)

## 
## Attaching package: 'hereR'

## The following objects are masked from 'package:ggmap':
## 
##     geocode, route

Cargamos un shapefile con las manzanas de Barcelona y nos quedamos con la columna que tiene la información espacial.

Bcn_manzanas <- read_sf("BCN_Parcelari/0601040100_Parcel·lari_POL_v.shp") %>% select(geometry)

st_crs(Bcn_manzanas)

## Coordinate Reference System:
##   User input: 25831 
##   wkt:
## PROJCS["ETRS89 / UTM zone 31N",
##     GEOGCS["ETRS89",
##         DATUM["European_Terrestrial_Reference_System_1989",
##             SPHEROID["GRS 1980",6378137,298.257222101,
##                 AUTHORITY["EPSG","7019"]],
##             TOWGS84[0,0,0,0,0,0,0],
##             AUTHORITY["EPSG","6258"]],
##         PRIMEM["Greenwich",0,
##             AUTHORITY["EPSG","8901"]],
##         UNIT["degree",0.0174532925199433,
##             AUTHORITY["EPSG","9122"]],
##         AUTHORITY["EPSG","4258"]],
##     PROJECTION["Transverse_Mercator"],
##     PARAMETER["latitude_of_origin",0],
##     PARAMETER["central_meridian",3],
##     PARAMETER["scale_factor",0.9996],
##     PARAMETER["false_easting",500000],
##     PARAMETER["false_northing",0],
##     UNIT["metre",1,
##         AUTHORITY["EPSG","9001"]],
##     AXIS["Easting",EAST],
##     AXIS["Northing",NORTH],
##     AUTHORITY["EPSG","25831"]]

Usamos set_key() para que here maps sea capaz de asociar las consultas con nuestra cuenta:

set_key("ooP8RLmyLNsZLelE8REZiGDXSLMd2YLiFpak_PqRtAc")

Ahora vamos a usar las isocronas para cada uno de los puntos:

isocronos_esttransporte <- map(1:nrow(transporte_esp), function(x) {isoline(transporte_esp[x,], mode = "pedestrian", range_type = "time",range = 60*6)})

Los juntamos en el mismo data frame:

isocronos_esttransporteJuntos <- bind_rows(isocronos_esttransporte)

Hacemos un gran poligono:

isocronos_esttransporteJuntosUnion <-  st_union(isocronos_esttransporteJuntos)

Vemos la cobertura que estimamos en un mapa de leaflet:

leaflet(isocronos_esttransporteJuntosUnion) %>% 
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons(color = "#762a83")

Vamos a proyectar en 2 dimensiones los datos que cargamos, usando la proyección en metros 25831, que es la que usa la ciudad de Barcelona:

isocronos_esttransporteJuntosUnionproy <- st_transform(isocronos_esttransporteJuntosUnion , crs= 25831)

st_crs(isocronos_esttransporteJuntosUnionproy)

## Coordinate Reference System:
##   User input: EPSG:25831 
##   wkt:
## PROJCS["ETRS89 / UTM zone 31N",
##     GEOGCS["ETRS89",
##         DATUM["European_Terrestrial_Reference_System_1989",
##             SPHEROID["GRS 1980",6378137,298.257222101,
##                 AUTHORITY["EPSG","7019"]],
##             TOWGS84[0,0,0,0,0,0,0],
##             AUTHORITY["EPSG","6258"]],
##         PRIMEM["Greenwich",0,
##             AUTHORITY["EPSG","8901"]],
##         UNIT["degree",0.0174532925199433,
##             AUTHORITY["EPSG","9122"]],
##         AUTHORITY["EPSG","4258"]],
##     PROJECTION["Transverse_Mercator"],
##     PARAMETER["latitude_of_origin",0],
##     PARAMETER["central_meridian",3],
##     PARAMETER["scale_factor",0.9996],
##     PARAMETER["false_easting",500000],
##     PARAMETER["false_northing",0],
##     UNIT["metre",1,
##         AUTHORITY["EPSG","9001"]],
##     AXIS["Easting",EAST],
##     AXIS["Northing",NORTH],
##     AUTHORITY["EPSG","25831"]]

class(isocronos_esttransporteJuntosUnionproy)

## [1] "sfc_MULTIPOLYGON" "sfc"

Lo convertimos en un objeto sf y creamos una columna, cobertura, a la cual le asignamos el valor TRUE:

isocronos_esttransporteJuntosUnionproy <- st_as_sf(isocronos_esttransporteJuntosUnionproy) %>%                                    mutate(cobertura = TRUE)

Unimos espacialmente el dataset de isocronos con el de las manzanas de la ciudad.

Bcn_manzanas <- st_join(Bcn_manzanas, isocronos_esttransporteJuntosUnionproy)

Completamos los datos para los casos en los cuales no hubo ningún resultado en el match:

Bcn_manzanas <- Bcn_manzanas %>% mutate(cobertura=ifelse(is.na(cobertura), FALSE, TRUE))

ggplot(Bcn_manzanas) +
  geom_sf(aes(fill=cobertura), alpha = .85, color=NA) +
  theme_minimal() +
  coord_sf(datum=NA) +
  scale_fill_manual(values = c("#762a83","#a6dba0"),
                    breaks = c(TRUE,FALSE),
                    labels=c("Menos de 6 minutos","Más de 6 minutos")) + 
  labs(title="Accesibilidad al transporte público en Barcelona",
       subtitle = "Estaciones de metro, tren y tranvía",
       fill = "Cobertura",
       caption = "Fuente: Ayuntamiento de Barcelona | 2020")

Autocorrelación espacial y Regresión lineal

Esto lo haremos en GeoDA y para eso modificamos la variable Cobertura true =1 y false = 0, del dataset anterior:

Bcn_manzanas$cobertura [Bcn_manzanas$cobertura == "TRUE"] <- 1
Bcn_manzanas$cobertura [Bcn_manzanas$cobertura == "FALSE"] <- 0
Bcn_manzanas$cobertura <- as.integer(Bcn_manzanas$cobertura)

Cargamos un dataset que brinda la página de airbnb con la oferta de alquileres en Barcelona. Además filtramos para quedarnos solo con 5 tipos de propiedad:

Bcn_airbnb <- read_sf("Bcn_airbnb.geojson") %>% 
  filter(city=="Barcelona") %>% 
  filter(property_type == "Bed & Breakfast" | property_type == "Apartment" | property_type == "Condominium" | property_type == "Boutique hotel" | property_type == "House") %>% 
  filter(price < 200) %>% 
  mutate(property_type_ponderado = property_type) %>% 
  select(price, property_type, property_type_ponderado, bedrooms, bathrooms)

Ponderamos el tipo de propiedad para que nos de un valor numérico y nos sirva para la regresión:

Bcn_airbnb$property_type_ponderado <- revalue(Bcn_airbnb$property_type_ponderado, c("Bed & Breakfast"= "1", "Apartment"="2", "Condominium" ="3", "Boutique hotel"="4", "House"="5"))

Bcn_airbnb$property_type_ponderado <- as.integer(Bcn_airbnb$property_type_ponderado)

st_crs(Bcn_airbnb)

## Coordinate Reference System:
##   User input: 4326 
##   wkt:
## GEOGCS["WGS 84",
##     DATUM["WGS_1984",
##         SPHEROID["WGS 84",6378137,298.257223563,
##             AUTHORITY["EPSG","7030"]],
##         AUTHORITY["EPSG","6326"]],
##     PRIMEM["Greenwich",0,
##         AUTHORITY["EPSG","8901"]],
##     UNIT["degree",0.0174532925199433,
##         AUTHORITY["EPSG","9122"]],
##     AUTHORITY["EPSG","4326"]]

Proyectamos en 2 dimensiones a los datos que cargamos, usando la proyección 25831:

Bcn_airbnb <- st_transform(Bcn_airbnb , crs= 25831)

Y unimos en un mismo dataset la cobertura del transporte público y la oferta de alojamientos publicadas por Airbnb.

Bcn_airbnb_transporte <- st_intersection(Bcn_airbnb, Bcn_manzanas)

## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

ggplot() +
  geom_sf(data=Bcn_manzanas, fill = "snow3", color = NA) +
  geom_sf(data=Bcn_airbnb_transporte %>% mutate(cobertura = as.character(cobertura)), aes(color = cobertura), size = .2, alpha = .7) + 
  guides(color = guide_legend(override.aes = list(size = 3, alpha = .85))) + 
  theme_void() +
  scale_color_manual(values = c("#762a83","#a6dba0"),
                      breaks = c(1,0),
                      labels=c("Menos de 6 minutos","Más de 6 minutos")) + 
  labs(title = "Cobertura de cada airbnb de transporte público",
       color = "Cobertura",
       caption = "Fuente: Ayuntamiento de Barcelona | 2020")

Exportamos este dataset para continuar trabajándolo en GeoDa:

st_write(Bcn_airbnb_transporte, "Bcn_airbnb_transporte.shp", driver = "ESRI Shapefile", delete_layer = TRUE)

## Warning in abbreviate_shapefile_names(obj): Field names abbreviated for ESRI
## Shapefile driver

## Deleting layer `Bcn_airbnb_transporte' using driver `ESRI Shapefile'
## Writing layer `Bcn_airbnb_transporte' to data source `Bcn_airbnb_transporte.shp' using driver `ESRI Shapefile'
## Writing 10724 features with 6 fields and geometry type Point.

Clustering

library(data.table)

## 
## Attaching package: 'data.table'

## The following objects are masked from 'package:dplyr':
## 
##     between, first, last

## The following object is masked from 'package:purrr':
## 
##     transpose

library(dbscan)

Utilizamos los datos de airbnb que preparamos anteriormente. Vemos lo que tenemos con leaflet.

paleta <- colorFactor(
  palette = c('#762a83','#af8dc3','#e7d4e8','#7fbf7b','#1b7837'),
  domain = Bcn_airbnb_transporte$property_type)

leaflet(Bcn_airbnb_transporte %>% st_transform(4326)) %>%
  addProviderTiles(providers$CartoDB.Positron)  %>% 
  addCircles(color = ~ paleta(property_type),
             opacity = .4,
             radius = ~ifelse(property_type == "Apartment", 10, 20),
             label = ~ price) %>% 
addLegend(title = "Oferta de Airbnb",
          position =  "bottomright",
          pal = paleta,
          values = ~ property_type)

En primera instancia, nos quedamos sólo con aquellas propiedades con buena cobertura al transporte público, que están catalogadas como “Apartamento”, y que posean menos de 3 dormitorios:

Cluster_airbnb_P <- Bcn_airbnb_transporte %>% 
  filter(property_type == "Apartment" & bedrooms < 3) %>% 
  filter(cobertura == "1")

paleta_bedrooms <- colorFactor(
  palette = c("#762a83", "#c2a5cf", "#5aae61"),
  domain = as.character(Cluster_airbnb_P$bedrooms))

leaflet(Cluster_airbnb_P %>% st_transform(4326)) %>%
  addProviderTiles(providers$CartoDB.Positron)  %>% 
  addCircles(color = ~ paleta_bedrooms(bedrooms),
             opacity = .4,
             radius = 20,
             label = ~ bedrooms) %>% 
addLegend(title = "Cantidad de Dormitorios",
          position =  "bottomright",
          pal = paleta_bedrooms,
          values = ~ bedrooms)

dbscan_Cluster_airbnb_P <- dbscan(x = st_distance(Cluster_airbnb_P),
                         eps = 3000,
                         minPts = 5)

dbscan_Cluster_airbnb_P

## DBSCAN clustering for 6490 objects.
## Parameters: eps = 3000, minPts = 5
## The clustering contains 167 cluster(s) and 2195 noise points.
## 
##    0    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15 
## 2195 1774    7    9   15    9    5    5    4    5   18    5  100   34   22    7 
##   16   17   18   19   20   21   22   23   24   25   26   27   28   29   30   31 
##    6   10   11    6    9    7    7    5   11    5   10    8   43    7    5    5 
##   32   33   34   35   36   37   38   39   40   41   42   43   44   45   46   47 
##   10    7    8   10   20   29   11   27    7    9    5    5   11    6    9    7 
##   48   49   50   51   52   53   54   55   56   57   58   59   60   61   62   63 
##    7   11    7   10   19   49    9    6   14   15    8   35    5  171    6    8 
##   64   65   66   67   68   69   70   71   72   73   74   75   76   77   78   79 
##    9    7   27    5    6    5    9    9    6   10   10    5    8   20   72    9 
##   80   81   82   83   84   85   86   87   88   89   90   91   92   93   94   95 
##    7    5    5   10   13    5    8    6   10    6  351  127    6   13   16    6 
##   96   97   98   99  100  101  102  103  104  105  106  107  108  109  110  111 
##    5   15   13   25    8   21    5    6   10    9    5    5   29    6    7   41 
##  112  113  114  115  116  117  118  119  120  121  122  123  124  125  126  127 
##    7    6   13    6    5    5    5    5    7    5    8    6   10    6    8    8 
##  128  129  130  131  132  133  134  135  136  137  138  139  140  141  142  143 
##   10    9    8   37   11   51    8   62   11    5    6    6   13   25    5   13 
##  144  145  146  147  148  149  150  151  152  153  154  155  156  157  158  159 
##    7    6    6   10    5    5    5    5    6    9    6    8    5    5    9    5 
##  160  161  162  163  164  165  166  167 
##    9    5    6    6    5    6    5    5 
## 
## Available fields: cluster, eps, minPts

Agregamos los 167 clusters que encontró y los 2195 puntos de ruido a nuestro dataset de oferta de Airbnb y graficamos con leaflet:

Cluster_airbnb_P <- Cluster_airbnb_P %>% 
  mutate(cluster = dbscan_Cluster_airbnb_P$cluster)

paleta_cluster <- colorFactor(palette = c("#bababa", '#40004b','#762a83','#9970ab','#c2a5cf','#e7d4e8','#a6dba0','#5aae61','#1b7837','#00441b'),
  domain = Cluster_airbnb_P$cluster)

leaflet(Cluster_airbnb_P %>% st_transform(4326)) %>%
  addProviderTiles(providers$CartoDB.Positron) %>% 
  addCircles(color = ~ paleta_cluster(cluster), 
             opacity = .5,
             radius = ~ifelse(cluster == "0", 10, 40),
             label = ~ price)

En segunda instancia, nos quedamos sólo con aquellas propiedades con buena cobertura al transporte público, que están catalogadas como “Apartamento”, pero que poseen 3 dormitorios o más:

Cluster_airbnb_G <- Bcn_airbnb_transporte %>% 
  filter(property_type == "Apartment" & bedrooms > 2) %>% 
  filter(cobertura == "1")

paleta_bedrooms <- colorFactor(
  palette = c('#762a83','#af8dc3','#7fbf7b','#1b7837'),
  domain = as.character(Cluster_airbnb_G$bedrooms))

leaflet(Cluster_airbnb_G %>% st_transform(4326)) %>%
  addProviderTiles(providers$CartoDB.Positron)  %>% 
  addCircles(color = ~ paleta_bedrooms(bedrooms),
             opacity = .4,
             radius = 20,
             label = ~ bedrooms) %>% 
addLegend(title = "Cantidad de Dormitorios",
          position =  "bottomright",
          pal = paleta_bedrooms,
          values = ~ bedrooms)

Con la función del paquete dbscan, asignamos la matriz de distancias con un valor de eps y un valor de minPts. Con st_distance() calculamos la distancia en metros y con minPts = 4 pedimos solo tres apartamentos que tengan más de 3 dormitorios, para determinar un núcleo.

dbscan_Cluster_airbnb_G <- dbscan(x = st_distance(Cluster_airbnb_G),
                         eps = 3000,
                         minPts = 5)

dbscan_Cluster_airbnb_G

## DBSCAN clustering for 861 objects.
## Parameters: eps = 3000, minPts = 5
## The clustering contains 24 cluster(s) and 276 noise points.
## 
##   0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19 
## 276 174  59  11 112  13   6  16  38  11  20   5  29   8   6  16   6   5  12   5 
##  20  21  22  23  24 
##   9   5   5   7   7 
## 
## Available fields: cluster, eps, minPts

Agregamos los 24 clusters que encontró y los 276 puntos de ruido a nuestro dataset de oferta de Airbnb y graficamos con leaflet:

Cluster_airbnb_G <- Cluster_airbnb_G %>% 
  mutate(cluster = dbscan_Cluster_airbnb_G$cluster)

paleta_cluster <- colorFactor(palette = c("#bababa", '#40004b','#762a83','#9970ab','#c2a5cf','#e7d4e8','#a6dba0','#5aae61','#1b7837','#00441b'),
  domain = Cluster_airbnb_G$cluster)

leaflet(Cluster_airbnb_G %>% st_transform(4326)) %>%
  addProviderTiles(providers$CartoDB.Positron) %>% 
  addCircles(color = ~ paleta_cluster(cluster), 
             opacity = .5,
             radius = ~ifelse(cluster == "0", 20, 40),
             label = ~ price)

También podemos ver la oferta de transporte según el dataset que utilizamos al principio. Graficamos con leaflet la oferta de transporte de la ciudad de Barcelona:

paleta_transporte <- colorFactor(palette = c('#40004b','#762a83','#9970ab','#c2a5cf','#a6dba0','#5aae61','#1b7837','#00441b'),
  domain = transporte_esp$NOM_CAPA)

leaflet(transporte_esp) %>%
  addProviderTiles(providers$CartoDB.Positron) %>% 
  addCircles(color = ~ paleta_transporte(NOM_CAPA), 
             opacity = .5,
             radius = 20,
             label = ~ EQUIPAMENT) %>% 
addLegend(title = "Oferta de Transporte",
          position =  "bottomright",
          pal = paleta_transporte,
          values = ~ NOM_CAPA)

Nos quedamos sólo con la oferta de metro y le pedimos a R que busque los cluster:

Cluster_transporte <- st_transform(transporte_esp, crs = 25831) %>%
  filter(NOM_CAPA == "Metro i línies urbanes FGC")

dbscan_transporte <- dbscan(x = st_distance(Cluster_transporte),
                        eps = 5000,
                        minPts = 5)

dbscan_transporte

## DBSCAN clustering for 463 objects.
## Parameters: eps = 5000, minPts = 5
## The clustering contains 24 cluster(s) and 235 noise points.
## 
##   0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19 
## 235   6  41   5   8  16  39   6   6   6   6   5  11   6   5  12   5   5   5   8 
##  20  21  22  23  24 
##   6   5   6   5   5 
## 
## Available fields: cluster, eps, minPts

Agregamos los 24 clusters que encontró y los 235 puntos de ruido a nuestro dataset de oferta de Airbnb y graficamos con leaflet:

Cluster_transporte <- Cluster_transporte %>% 
  mutate(cluster = dbscan_transporte$cluster)

paleta_cluster <- colorFactor(palette = c("#bababa", '#40004b','#762a83','#9970ab','#c2a5cf','#e7d4e8','#a6dba0','#5aae61','#1b7837','#00441b'),
  domain = Cluster_transporte$cluster)

leaflet(Cluster_transporte %>% st_transform(4326)) %>%
  addProviderTiles(providers$CartoDB.Positron) %>% 
  addCircles(color = ~ paleta_cluster(cluster), 
             opacity = .5,
             radius = ~ifelse(cluster == "0", 20, 40),
             label = ~ EQUIPAMENT)

Ahora obtenemos los polígonos formados por los clusters de la oferta de Airbnb de apartamentos con más de 3 dormitorios:

Cluster_pol_airbnbG <- lapply(Cluster_airbnb_G %>% filter(!cluster == 0) %>%  pull(cluster) %>% unique(),
                        function(clusterNumber){
  st_convex_hull(x = Cluster_airbnb_G %>%
                   filter(cluster %in% clusterNumber) %>%
                   st_union()) %>% 
                            st_as_sf() %>%
                            mutate(cluster=clusterNumber)}) %>% bind_rows()

Cluster_pol_airbnbG <- Cluster_pol_airbnbG %>% mutate(categoria="Oferta Airbnb: Apartamento > 3 Dormitorios")

leaflet(Cluster_pol_airbnbG %>% st_transform(4326)) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons(color = "#762a83")

Hacemos lo mismo con la oferta de Airbnb de apartamentos con menos de 3 Dormitorios:

Cluster_pol_airbnbP <- lapply(Cluster_airbnb_P %>% filter(!cluster == 0) %>%  pull(cluster) %>% unique(),
                        function(clusterNumber){
  st_convex_hull(x = Cluster_airbnb_P %>%
                   filter(cluster %in% clusterNumber) %>%
                   st_union()) %>% 
                            st_as_sf() %>%
                            mutate(cluster=clusterNumber)}) %>% bind_rows()

Cluster_pol_airbnbP <- Cluster_pol_airbnbP %>% mutate(categoria="Oferta Airbnb: Apartamento < 3 Dormitorios")

leaflet(Cluster_pol_airbnbP %>% st_transform(4326)) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons(color = "#c2a5cf")

Hacemos lo mismo con la oferta de Metro:

Cluster_pol_transporte <- lapply(Cluster_transporte %>% filter(!cluster == 0) %>%  pull(cluster) %>% unique(),
                                 function(clusterNumber){
                                   st_convex_hull(x = Cluster_transporte %>%
                                                    filter(cluster %in% clusterNumber) %>%
                                                    st_union()) %>% 
                                     st_as_sf() %>%
                                     mutate(cluster=clusterNumber)}) %>% bind_rows()

Cluster_pol_transporte <- Cluster_pol_transporte %>% mutate(categoria="Oferta Transporte: Metro")

leaflet(Cluster_pol_transporte %>% st_transform(4326)) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addPolygons(color = "#a6dba0")

Por último unimos los polígonos y graficamos:

Cluster_airbnb_metro <- rbind(Cluster_pol_airbnbG, Cluster_pol_airbnbP, Cluster_pol_transporte)