M1/ Geoprocesamiento

Comenzamos cargando las librerías a utilizar

library(tidyverse)
## Warning: package 'tidyverse' was built under R version 3.6.1
## -- Attaching packages --------------------------------------------------------------------------------------------------- tidyverse 1.2.1 --
## v ggplot2 3.2.1     v purrr   0.3.2
## v tibble  2.1.1     v dplyr   0.8.1
## v tidyr   0.8.3     v stringr 1.4.0
## v readr   1.3.1     v forcats 0.4.0
## Warning: package 'ggplot2' was built under R version 3.6.1
## -- Conflicts ------------------------------------------------------------------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(ggplot2)
library(sf)
## Warning: package 'sf' was built under R version 3.6.1
## Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
library(rgdal)
## Warning: package 'rgdal' was built under R version 3.6.1
## Loading required package: sp
## Warning: package 'sp' was built under R version 3.6.1
## rgdal: version: 1.4-4, (SVN revision 833)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
##  Path to GDAL shared files: C:/Users/fsciutto/Documents/R/R-3.6.0/library/rgdal/gdal
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: C:/Users/fsciutto/Documents/R/R-3.6.0/library/rgdal/proj
##  Linking to sp version: 1.3-1
require(PBSmapping)
## Loading required package: PBSmapping
## Warning: package 'PBSmapping' was built under R version 3.6.1
## 
## -----------------------------------------------------------
## PBS Mapping 2.72.1 -- Copyright (C) 2003-2019 Fisheries and Oceans Canada
## 
## PBS Mapping comes with ABSOLUTELY NO WARRANTY;
## for details see the file COPYING.
## This is free software, and you are welcome to redistribute
## it under certain conditions, as outlined in the above file.
## 
## A complete user guide 'PBSmapping-UG.pdf' is located at 
## C:/Users/fsciutto/Documents/R/R-3.6.0/library/PBSmapping/doc/PBSmapping-UG.pdf
## 
## Packaged on 2019-03-14
## Pacific Biological Station, Nanaimo
## 
## All available PBS packages can be found at
## https://github.com/pbs-software
## 
## To see demos, type '.PBSfigs()'.
## -----------------------------------------------------------
require(maptools)
## Loading required package: maptools
## Warning: package 'maptools' was built under R version 3.6.1
## Checking rgeos availability: FALSE
##      Note: when rgeos is not available, polygon geometry     computations in maptools depend on gpclib,
##      which has a restricted licence. It is disabled by default;
##      to enable gpclib, type gpclibPermit()
library(ggrepel)
## Warning: package 'ggrepel' was built under R version 3.6.1
library(extrafont)
## Registering fonts with R
library(leaflet)
## Warning: package 'leaflet' was built under R version 3.6.1
library(osmdata)
## Warning: package 'osmdata' was built under R version 3.6.1
## Data (c) OpenStreetMap contributors, ODbL 1.0. http://www.openstreetmap.org/copyright
library(tidyverse)

Cargamos los shapes a unir.

barrios <- st_read("./00BASES/Boston_Neighborhoods.geojson")
## Reading layer `Boston_Neighborhoods' from data source `D:\ELISA\[B]ACADEMICA\FORMACION\[02]MEU\CIENCIA DE DATOS II - A BRUST\00BASES\Boston_Neighborhoods.geojson' using driver `GeoJSON'
## Simple feature collection with 26 features and 7 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -71.19125 ymin: 42.22792 xmax: -70.92278 ymax: 42.39699
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
escuelas <- st_read ("./00BASES/Public_Schools.geojson")
## Reading layer `Public_Schools' from data source `D:\ELISA\[B]ACADEMICA\FORMACION\[02]MEU\CIENCIA DE DATOS II - A BRUST\00BASES\Public_Schools.geojson' using driver `GeoJSON'
## Simple feature collection with 131 features and 16 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -71.17434 ymin: 42.2339 xmax: -71.00412 ymax: 42.39162
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
ggplot()+
  geom_sf(data= barrios)+
  geom_sf(data=escuelas)

A los efectos de futura visualizacion vamos a hallar los centroides de los poligonos de barrios.

barriosCent<- st_centroid(barrios)
## Warning in st_centroid.sf(barrios): st_centroid assumes attributes are
## constant over geometries of x
## Warning in st_centroid.sfc(st_geometry(x), of_largest_polygon =
## of_largest_polygon): st_centroid does not give correct centroids for
## longitude/latitude data
ggplot()+
    geom_sf(data=barriosCent)

Ahora vamos a hacer el join espacial entre la data de barrios y escuelas. Primero, asignamos la informacion de barrios a los puntos de las escuelas.

escuelas <- st_join(escuelas, barrios)
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
## although coordinates are longitude/latitude, st_intersects assumes that they are planar

Luego agrupamos las escuelas por barrio y contabilizamos. Nota: anulamos en este paso la informacion geometrica del shape, que traiamos residualmente para que no interfiera en el paso siguiente.

escuelasXbarrio <- escuelas %>% 
  group_by(Name)%>%
  summarise(total=n()) %>% 
  st_set_geometry(NULL)

Acto seguido, hacemos un left join para asignarle al shape de barrios la data de escuelas que le corresponde a cada poligono. Visualizamos a partir de un “Choropleth map” de escuelas por barrio. Nota: Para mapear las etiquetas reordenamos las cordenadas en x, y.

barriosCent_xy <- as.data.frame(sf::st_coordinates(barriosCent))
barriosCent_xy$Name <- barrios$Name
barrios <- left_join(barrios,escuelasXbarrio)
## Joining, by = "Name"
barrios %>% 
  ggplot() + 
  geom_sf(aes(fill = total), color = NA)+
  geom_sf(data=barriosCent, size=.5)+
  geom_text_repel(data = barriosCent_xy, aes(X, Y, label = Name), colour = "black", size=3)+
  theme_void()+
  labs(title= "Escuelas por barrio", subtitle = "Boston, 2019", x=NULL, y=NULL, caption = "Fuente: Analize Boston")+
  theme(text=element_text(size=12, family="Century"))+
  scale_fill_distiller(palette= "YlGn", direction= 1)

A partir del grafico observamos que hay una mayor concentracion de escuelas en el barrio de Dorchester y que hay barrios que no poseen ninguna.

M2 / Acceso a información urbana geo-referenciada

Descargamos la boundign box de la ciudad de Boston

bbox <- getbb("Boston, US")
bbox
##         min       max
## x -71.19126 -70.80449
## y  42.22791  42.39698

Obtenemos un rectangulo, de una superficie mayor a la ciudad. Para eliminar la data innecesaria, la funcion format_out me permite extraer el poligono perimetral, es decir, los limites politicos.

bbox_poly <-getbb("Boston , US", format_out = "sf_polygon")
bbox_poly
## Simple feature collection with 1 feature and 0 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -71.19126 ymin: 42.22791 xmax: -70.80449 ymax: 42.39698
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
##                         geometry
## 1 POLYGON ((-71.19126 42.2828...

Hacemos en primera instancia un mapeo interactivo con el fin de verificar si el poligono se corresponde con el area deseada.

leaflet()%>%
  addTiles()%>%
  addPolygons(data=bbox_poly)
Boston <- opq(bbox) %>%
                add_osm_feature(key= "highway")%>%
                osmdata_sf

Obtenemos el shape del layer de calles a partir del recorte original.

calles <- Boston$osm_lines
View(calles)
ggplot()+
  geom_sf(data=calles)

Clipeamos las calles con el perimetro de la ciudad, es decir sus limites politicos.

calles <- st_intersection (calles, bbox_poly)
## 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()+
  theme_void()+
  geom_sf(data=calles)


#consultar el blog r-spatial.github.io/sf/articles/sf3.html

Ahora vamos a graficar el shape de calles obtenido según alguno de sus atributos. Para ello, comenzamos explorando el dataframe.

str(calles)
## Classes 'sf' and 'data.frame':   21362 obs. of  275 variables:
##  $ osm_id                               : chr  "4790718" "4790735" "8614591" "8614837" ...
##  $ name                                 : chr  "Skybridge to Hilton" NA "Charles River Dam Road" "Charlestown Avenue" ...
##  $ FIXME                                : chr  NA NA NA NA ...
##  $ FIXME.hgv                            : chr  NA NA NA NA ...
##  $ NHS                                  : chr  NA NA NA NA ...
##  $ abutters                             : chr  NA NA NA NA ...
##  $ access                               : chr  NA NA NA NA ...
##  $ access.conditional                   : chr  NA NA NA NA ...
##  $ addr.city                            : chr  NA NA NA NA ...
##  $ addr.housename                       : chr  NA NA NA NA ...
##  $ addr.housenumber                     : chr  NA NA NA NA ...
##  $ addr.postcode                        : chr  NA NA NA NA ...
##  $ addr.state                           : chr  NA NA NA NA ...
##  $ addr.street                          : chr  NA NA NA NA ...
##  $ alt_name                             : chr  NA NA NA NA ...
##  $ amenity                              : chr  NA NA NA NA ...
##  $ architect                            : chr  NA NA NA NA ...
##  $ area                                 : chr  NA NA NA NA ...
##  $ artwork_type                         : chr  NA NA NA NA ...
##  $ attribution                          : chr  NA NA "Office of Geographic and Environmental Information (MassGIS)" "Office of Geographic and Environmental Information (MassGIS)" ...
##  $ barrier                              : chr  NA NA NA NA ...
##  $ bic                                  : chr  NA NA NA NA ...
##  $ bicycle                              : chr  NA NA "yes" NA ...
##  $ bicycle.lanes.backward               : chr  NA NA NA NA ...
##  $ bridge                               : chr  "yes" NA NA "yes" ...
##  $ bridge.material                      : chr  NA NA NA NA ...
##  $ bridge.movable                       : chr  NA NA NA NA ...
##  $ bridge.name                          : chr  NA NA NA NA ...
##  $ bridge.structure                     : chr  NA NA NA NA ...
##  $ building                             : chr  NA NA NA NA ...
##  $ building.architecture                : chr  NA NA NA NA ...
##  $ building.colour                      : chr  NA NA NA NA ...
##  $ building.height                      : chr  NA NA NA NA ...
##  $ building.levels                      : chr  NA NA NA NA ...
##  $ building.material                    : chr  NA NA NA NA ...
##  $ building.part                        : chr  NA NA NA NA ...
##  $ bus                                  : chr  NA NA NA NA ...
##  $ busway                               : chr  NA NA NA NA ...
##  $ busway.right                         : chr  NA NA NA NA ...
##  $ campusbuilding                       : chr  NA NA NA NA ...
##  $ check_date                           : chr  NA NA NA NA ...
##  $ class.bicycle                        : chr  NA NA NA NA ...
##  $ comment                              : chr  NA NA NA NA ...
##  $ condition                            : chr  NA NA "fair" "good" ...
##  $ construction                         : chr  NA NA NA NA ...
##  $ conveying                            : chr  NA NA NA NA ...
##  $ covered                              : chr  NA NA NA NA ...
##  $ created_by                           : chr  NA NA NA NA ...
##  $ crossing                             : chr  NA NA NA NA ...
##  $ crossing.island                      : chr  NA NA NA NA ...
##  $ crossing_1                           : chr  NA NA NA NA ...
##  $ crossing_ref                         : chr  NA NA NA NA ...
##  $ cutting                              : chr  NA NA NA NA ...
##  $ cycleway                             : chr  NA NA NA NA ...
##  $ cycleway.both                        : chr  NA NA NA NA ...
##  $ cycleway.both.lane                   : chr  NA NA NA NA ...
##  $ cycleway.lanes.backward              : chr  NA NA NA NA ...
##  $ cycleway.left                        : chr  NA NA NA NA ...
##  $ cycleway.left.buffer                 : chr  NA NA NA NA ...
##  $ cycleway.left.lane                   : chr  NA NA NA NA ...
##  $ cycleway.right                       : chr  NA NA "track" NA ...
##  $ cycleway.right.buffer                : chr  NA NA NA NA ...
##  $ cycleway.right.lane                  : chr  NA NA NA NA ...
##  $ cycleway.right_1                     : chr  NA NA NA NA ...
##  $ denomination                         : chr  NA NA NA NA ...
##  $ description                          : chr  NA NA NA NA ...
##  $ designation                          : chr  NA NA NA NA ...
##  $ destination                          : chr  NA NA NA NA ...
##  $ destination.backward                 : chr  NA NA NA NA ...
##  $ destination.forward                  : chr  NA NA NA NA ...
##  $ destination.lanes                    : chr  NA NA NA NA ...
##  $ destination.ref                      : chr  NA NA NA NA ...
##  $ destination.ref.lanes                : chr  NA NA NA NA ...
##  $ destination.ref.to                   : chr  NA NA NA NA ...
##  $ destination.street                   : chr  NA NA NA NA ...
##  $ destination.symbol                   : chr  NA NA NA NA ...
##  $ direction                            : chr  NA NA NA NA ...
##  $ driveway                             : chr  NA NA NA NA ...
##  $ ele                                  : chr  NA NA NA NA ...
##  $ embankment                           : chr  NA NA NA NA ...
##  $ emergency                            : chr  NA NA NA NA ...
##  $ fixme                                : chr  NA NA NA NA ...
##  $ floating                             : chr  NA NA NA NA ...
##  $ foot                                 : chr  "yes" "yes" NA NA ...
##  $ footway                              : chr  NA NA NA NA ...
##  $ ford                                 : chr  NA NA NA NA ...
##  $ former.name                          : chr  NA NA NA NA ...
##  $ gnis.county_id                       : chr  NA NA NA NA ...
##  $ gnis.county_name                     : chr  NA NA NA NA ...
##  $ gnis.created                         : chr  NA NA NA NA ...
##  $ gnis.feature_id                      : chr  NA NA NA NA ...
##  $ gnis.import_uuid                     : chr  NA NA NA NA ...
##  $ gnis.reviewed                        : chr  NA NA NA NA ...
##  $ gnis.state_id                        : chr  NA NA NA NA ...
##  $ golf                                 : chr  NA NA NA NA ...
##  $ golf_cart                            : chr  NA NA NA NA ...
##  $ goods                                : chr  NA NA NA NA ...
##  $ handrail                             : chr  NA NA NA NA ...
##  $ handrail.center                      : chr  NA NA NA NA ...
##   [list output truncated]
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
##   ..- attr(*, "names")= chr  "osm_id" "name" "FIXME" "FIXME.hgv" ...

Nos interesa graficar la variable “condition”, donde se expresa el caracter de la avenida según su uso.

calles$condition <- as.factor(calles$condition)
levels(calles$condition)
## [1] "deficient"      "fair"           "fair;deficient" "good"          
## [5] "intolerable"    "poor"
calles$condition <- factor(calles$condition, levels=c( "intolerable", "poor", "deficient" , "fair;deficient",  "fair" , "good" ))

I. Al graficar las calles segun su condicion, observamos que, en general, las calles se encuentran en buen estado pero que, llamativamente, las calles en mal estado se agrupan en determinadas zonas de la ciudad.

calles <- calles %>% filter(!is.na(condition))
colores <-c("good"="light grey","fair" ="grey", "fair;deficient"="#fecc5c","deficient" = "#fd8d3c", "poor" = "#f03b20","intolerable" = "#bd0026" )
 
ggplot()+
  geom_sf(data=calles, aes(color=condition), fill=NA, size=.01)+
  scale_color_manual(values=colores, aesthetics = c("color"))+
  theme_void()+
  labs(title= "Estado de las calles",
       subtitle="Boston 2019",
       caption= "Fuente: OpenStreet Map")+
  theme(text=element_text(size=12, family="Century"))

II. A partir de las categorias explicadas en wiki.openstreet map, nos descargamos la variable “schools” del openStreet map

escuelasOst<-opq(bbox) %>%
  add_osm_feature(key="amenity", value="school") %>%
  osmdata_sf
escuelasOst <- escuelasOst$osm_points
escuelasOst <- st_intersection (escuelasOst, bbox_poly)
## 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=barrios)+
  geom_sf(data=escuelasOst, size=0.5)+
  theme_void()+
  labs(title= "Escuelas por barrio",
       subtitle="Boston 2019",
       caption= "Fuente: OpenStreet Map")+
  theme(text=element_text(size=12, family="Century"))

II.a. Al mapear las escuelas pareceria ser que su distribucion es homogenea a lo largo de todo el shape de barrios. Sin embargo, es dudosa la geometria obtenida es dudosa, y no queda claro si cada punto corresponde o no a un establecimiento. Para verificar cuantitativamente vamos a hacer un conteo de las escuelas por barrio, utilizando nuevamente el shape de barrios oficiales y los cruzamos con la data de Open StreetMap.

escuelasOstXbarrio <- st_join(escuelasOst, barrios)
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
## although coordinates are longitude/latitude, st_intersects assumes that they are planar

Luego agrupamos las escuelas por barrio y contabilizamos. Nota: anulamos en este paso la informacion geometrica del shape, que traiamos residualmente para que no interfiera en el paso siguiente.

escuelasOstXbarrio <- escuelasOstXbarrio %>% 
  group_by(Name)%>%
  summarise(total=n()) %>% 
  st_set_geometry(NULL)

Ahora, asignamos a cada barrio el conteo de escuelas y graficamos

barriosOst <- barrios %>% 
  left_join(escuelasOstXbarrio)
## Joining, by = c("Name", "total")
ggplot(barriosOst) + 
  geom_sf(aes(fill = total), color = NA)+
  geom_sf(data=barriosCent, size=.5)+
  geom_text_repel(data = barriosCent_xy, aes(X, Y, label = Name), colour = "black", size=3)+
  theme_void()+
  labs(title= "Escuelas por barrio",
       subtitle="Boston 2019",
       caption= "Fuente: OpenStreet Map")+
  theme(text=element_text(size=12, family="Century"))+
  scale_fill_distiller(palette= "YlGn", direction = 1)

II. b. Al hacer el conteo observamos que el numero de escuelas es llamativamene mayor en los datos de open street maps. esto se debe al dudoso origen de los datos, que deberan ser filtrados para nuevamente,hacer una comparativa mas fidedigna y poder sacar conclusiones.

M3 / CAPTURANDO Y EXPLORANDO DATOS DE TWEETER

Cargamos las librerias

library(tidyverse)
library(rtweet)
## Warning: package 'rtweet' was built under R version 3.6.1
## 
## Attaching package: 'rtweet'
## The following object is masked from 'package:purrr':
## 
##     flatten
library(ggmap)
## Warning: package 'ggmap' was built under R version 3.6.1
## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.
## Please cite ggmap if you use it! See citation("ggmap") for details.
tweets <- read.csv("./00BASES/tweets_rmba.csv")

Exploramos las variables del dataset

names(tweets)
##  [1] "X"                     "Y"                    
##  [3] "id"                    "in_reply_to_status_id"
##  [5] "in_reply_to_user_id"   "text"                 
##  [7] "created"               "lang"                 
##  [9] "source"                "user_name"            
## [11] "user_id"               "user_created"         
## [13] "user_description"      "user_location"        
## [15] "user_followers"        "user_followed"

Visualizamos el top 5 de usuarios con mas seguidores, asumiendo que sus tweets son los que tienen mayor repercusión.

tweets %>% 
    top_n(5, user_followers) %>% 
    arrange(desc(user_followers)) %>% 
    select(user_name, user_followers, text)
##   user_name user_followers
## 1    maluma        4924516
## 2    maluma        4919851
## 3    maluma        4912535
## 4    maluma        4912499
## 5    maluma        4912441
##                                                                                                                                                                   text
## 1                                              Gracias hermosa por tenerme en tu show.. la pasamos INCREÃ\215BLE!! @Su_Gimenez @ Buenos Aires,â\200¦ https://t.co/UiAzGZ7JRM
## 2                                                         #â\235Œ THE FILM. Ya lo vieron ?? Vayan a @vevo / go to @vevo @ Buenos Aires, Argentina https://t.co/ymm7BQVQfx
## 3                                                                   "Si con #23..." LOS AMO ARGENTINA 🇦🇷â\235¤ï¸\217 @ Buenos Aires, Argentina https://t.co/d5JDYoXIBX
## 4 15.000 almas unidas ðŸ\231\217ðŸ\217»ðŸ\231\217ðŸ\217»ðŸ\231\217ðŸ\217»  \nEstás imágenes hablan por sí solas!! \n#MalumaArgentinaTour 🇦🇷🇦🇷🇦🇷 @â\200¦ https://t.co/BfUAqE5zKF
## 5                                YESSSS!!! ✊ðŸ\217»âœŠðŸ\217»âœŠðŸ\217»\nPrimera noche en el hipódromo de Buenos Aires. Una noche que va a quedarâ\200¦ https://t.co/vUpUMl0yCL

Graficamos la distribución de los usuarios según su popularidad ( cantidad de seguidores).

options(scipen = 20)
ggplot(tweets) +
    geom_histogram(aes(x = user_followers), fill= "blue")+
  theme_minimal()+
   labs(title= "Top 5 de usuarios",
       subtitle="Twits RMBA",
       caption= "Fuente: Twitter",
       x= "Seguidores",
       y= NULL)+
  theme(text=element_text(size=12, family="Century"))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

#write.csv (tweets, "~/tweets2.csv) para guardar como con un nuevo nombre
#db scan.
tweets <-  tweets %>% 
    filter(user_followers>5000)
bboxBA <- c(min(tweets$X),
          min(tweets$Y),
          max(tweets$X),
          max(tweets$Y))
mapaBa <- get_stamenmap(bboxBA, maptype = "toner-lite")
## Source : http://tile.stamen.com/toner-lite/10/343/614.png
## Source : http://tile.stamen.com/toner-lite/10/344/614.png
## Source : http://tile.stamen.com/toner-lite/10/345/614.png
## Source : http://tile.stamen.com/toner-lite/10/346/614.png
## Source : http://tile.stamen.com/toner-lite/10/347/614.png
## 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
ggmap(mapaBa)

ggmap(mapaBa)+
  geom_point(data=tweets, aes(x=X,y=Y, size=user_followers, color=user_followers))+
  scale_color_distiller(palette = "Spectral", direction = -1)+
  theme_void()+
  labs(title= "Distribución de tweets populares",
       subtitle="Tweets RMBA",
       caption= "Fuente: Twitter")+
  theme(text=element_text(size=12, family="Century"))

top_tweets <-  tweets %>% 
    top_n(20, user_followers)
ggmap(mapaBa)+
  geom_point(data=top_tweets, aes(x=X,y=Y, color=user_followers))+
  scale_color_distiller(palette = "Spectral", direction = -1, name= "Seguidores")+
  geom_label_repel(data =top_tweets, aes(X, Y, label = user_name))+
  theme_void()+
  labs(title= "Distribución de usuarios populares",
       subtitle="Twits RMBA",
       caption= "Fuente: Twitter")+
  theme(text=element_text(size=12, family="Century"))
## Warning in min(x): ningún argumento finito para min; retornando Inf
## Warning in max(x): ningun argumento finito para max; retornando -Inf
## Warning in min(x): ningún argumento finito para min; retornando Inf
## Warning in max(x): ningun argumento finito para max; retornando -Inf

M4 / ANALIZANDO PATRONES TEMPORALES

Vamos a trabajar con la base de datos de delitos de la ciudad de Boston, obtenida del portal de datos abiertos.

delitos <- read.csv ("./00BASES/tmpcl2rr915.csv")

Exploramos el dataset

summary(delitos)
##       INCIDENT_NUMBER    OFFENSE_CODE 
##  I152071596   :    20   Min.   : 111  
##  I172053750   :    18   1st Qu.:1102  
##  I192025403   :    15   Median :3001  
##  I162067346   :    14   Mean   :2328  
##  I182051210   :    14   3rd Qu.:3201  
##  I130041200-00:    13   Max.   :3831  
##  (Other)      :414855                 
##                        OFFENSE_CODE_GROUP
##  Motor Vehicle Accident Response: 48221  
##  Larceny                        : 33721  
##  Medical Assistance             : 31699  
##  Investigate Person             : 24220  
##  Other                          : 23358  
##  Drug Violation                 : 21494  
##  (Other)                        :232236  
##                             OFFENSE_DESCRIPTION    DISTRICT     
##  SICK/INJURED/MEDICAL - PERSON        : 25480   B2     : 66216  
##  INVESTIGATE PERSON                   : 24225   C11    : 55883  
##  M/V - LEAVING SCENE - PROPERTY DAMAGE: 21166   D4     : 53432  
##  ASSAULT SIMPLE - BATTERY             : 19292   B3     : 46934  
##  VANDALISM                            : 19173   A1     : 46402  
##  VERBAL DISPUTE                       : 17286   C6     : 30185  
##  (Other)                              :288327   (Other):115897  
##  REPORTING_AREA  SHOOTING              OCCURRED_ON_DATE       YEAR     
##  Min.   :  0.0    :413237   2016-08-01 00:00:00:    33   Min.   :2015  
##  1st Qu.:178.0   Y:  1712   2017-06-01 00:00:00:    33   1st Qu.:2016  
##  Median :345.0              2015-07-01 00:00:00:    27   Median :2017  
##  Mean   :384.2              2017-08-01 00:00:00:    26   Mean   :2017  
##  3rd Qu.:544.0              2015-12-07 11:38:00:    25   3rd Qu.:2018  
##  Max.   :962.0              2017-01-01 00:00:00:    23   Max.   :2019  
##  NA's   :26484              (Other)            :414782                 
##      MONTH          DAY_OF_WEEK         HOUR             UCR_PART     
##  Min.   : 1.00   Friday   :63049   Min.   : 0.00             :   109  
##  1st Qu.: 4.00   Monday   :59426   1st Qu.: 9.00   Other     :  1600  
##  Median : 7.00   Saturday :58393   Median :14.00   Part One  : 78522  
##  Mean   : 6.59   Sunday   :52357   Mean   :13.11   Part Three:208137  
##  3rd Qu.: 9.00   Thursday :60535   3rd Qu.:18.00   Part Two  :126581  
##  Max.   :12.00   Tuesday  :60192   Max.   :23.00                      
##                  Wednesday:60997                                      
##             STREET            Lat             Long       
##  WASHINGTON ST : 18798   Min.   :-1.00   Min.   :-71.18  
##                : 12209   1st Qu.:42.30   1st Qu.:-71.10  
##  BLUE HILL AVE : 10279   Median :42.33   Median :-71.08  
##  BOYLSTON ST   :  9296   Mean   :42.22   Mean   :-70.92  
##  DORCHESTER AVE:  6546   3rd Qu.:42.35   3rd Qu.:-71.06  
##  TREMONT ST    :  6422   Max.   :42.40   Max.   : -1.00  
##  (Other)       :351399   NA's   :26853   NA's   :26853   
##                         Location     
##  (0.00000000, 0.00000000)   : 26853  
##  (42.34862382, -71.08277637):  1652  
##  (42.36183857, -71.05976489):  1620  
##  (42.28482577, -71.09137369):  1425  
##  (42.32866284, -71.08563401):  1310  
##  (42.25621592, -71.12401947):  1210  
##  (Other)                    :380879
head(delitos)
##   INCIDENT_NUMBER OFFENSE_CODE      OFFENSE_CODE_GROUP
## 1      I192066149         3125         Warrant Arrests
## 2      I192066145         3106 Property Related Damage
## 3      I192066142         3301         Verbal Disputes
## 4      I192066140         3115      Investigate Person
## 5      I192066139         2907              Violations
## 6      I192066138         3301         Verbal Disputes
##               OFFENSE_DESCRIPTION DISTRICT REPORTING_AREA SHOOTING
## 1                  WARRANT ARREST                      NA         
## 2    PROPERTY - ACCIDENTAL DAMAGE      E18            556         
## 3                  VERBAL DISPUTE      C11            347         
## 4              INVESTIGATE PERSON       B2            311         
## 5 VAL - OPERATING AFTER REV/SUSP.      D14            763         
## 6                  VERBAL DISPUTE      C11            356         
##      OCCURRED_ON_DATE YEAR MONTH DAY_OF_WEEK HOUR   UCR_PART        STREET
## 1 2019-08-21 21:31:00 2019     8   Wednesday   21 Part Three              
## 2 2019-08-21 20:09:00 2019     8   Wednesday   20 Part Three   PINEDALE RD
## 3 2019-08-21 21:08:00 2019     8   Wednesday   21 Part Three       EAST ST
## 4 2019-08-21 20:52:00 2019     8   Wednesday   20 Part Three     WARREN ST
## 5 2019-08-21 20:28:00 2019     8   Wednesday   20   Part Two WASHINGTON ST
## 6 2019-08-21 20:47:00 2019     8   Wednesday   20 Part Three    CHARLES ST
##        Lat      Long                    Location
## 1       NA        NA    (0.00000000, 0.00000000)
## 2 42.27685 -71.12171 (42.27684989, -71.12170659)
## 3 42.30761 -71.05971 (42.30761110, -71.05970561)
## 4 42.31941 -71.08141 (42.31940816, -71.08140505)
## 5 42.34904 -71.15748 (42.34903684, -71.15747989)
## 6 42.29991 -71.06425 (42.29990887, -71.06424768)

II. Realizamos un primer grafico que ilustre los delitos en funcion año en que sucedieron. Observamos un aumento de los reportes entre los años 2016 a 2018, y un descenso en 2019. El descenso en 2019 se podria deber - probablemente- a que se trata del corriente año y aun no ha transcurrido en su totalidad. Nos falta informacion para dilucidar el incremento respecto de 2015.

options(scipen = 20)

ggplot(delitos) + 
  geom_bar(aes(x = YEAR), fill= "red")+
   theme_minimal()+
   labs(title= "Delitos reportados",
       subtitle="Boston 2015-2019",
       caption= "Fuente: Analyze Boston",
       x= "Year",
       y= NULL)+
  theme(text=element_text(size=12, family="Century"))

Para indagar un poco mas sobre el tema vamos a graficar los delitos en funcion de su categoria, a ver si esta informacion nos aporta alguna pista.

delitos %>% 
    count(OFFENSE_CODE_GROUP) %>% 
    top_n(5) %>% 
    arrange(desc(n))
## Selecting by n
## # A tibble: 5 x 2
##   OFFENSE_CODE_GROUP                  n
##   <fct>                           <int>
## 1 Motor Vehicle Accident Response 48221
## 2 Larceny                         33721
## 3 Medical Assistance              31699
## 4 Investigate Person              24220
## 5 Other                           23358
delitos_frecuentes <- delitos %>% 
    count(OFFENSE_CODE_GROUP) %>% 
    top_n(5) %>% 
    pull(OFFENSE_CODE_GROUP)
## Selecting by n

Al graficar observamos que el incremento de los reportes esta asociado en gran medida a los accidentes automobilisticos. Para darnos cuenta si es que efectivamente ha habido un incremento de todos modos, deberiamos contar con informacion historica de un rango mayor de años, ya que un pequeño salto entre dos años podria no tener que ver con la tendencia general.

delitos %>% 
    filter(OFFENSE_CODE_GROUP %in% delitos_frecuentes)%>% 
    ggplot() +
        geom_bar(aes(x = YEAR, fill = OFFENSE_CODE_GROUP))+
   theme_minimal()+
   labs(title= "Delitos frecuentes ",
       subtitle="Boston 2015-2019",
       caption= "Fuente: Analyze Boston",
       x= "Año",
       y= NULL)+
  theme(text=element_text(size=12, family="Century"))+
  scale_fill_brewer(palette= "Set1", name= "Tipo de delito")

Observamos que la proporcion en el tiempo se mantiene,y que el delito mas frecuente suele ser el accidente vehicular.

Pasamos entonces a ver si hay algun tipo de fluctuacion relacionada con el mes del año.


# Primero realizamos un conteo de delitos por tipo y por mes del año
conteo <-  delitos %>% 
    filter( OFFENSE_CODE_GROUP %in% delitos_frecuentes) %>% 
    count(OFFENSE_CODE_GROUP, MONTH)

# Y ahora a mostras las cantidades mensuales como líneas 

meses<- seq(1,12, by= 1)
ggplot(conteo) +
   geom_line(aes(x = MONTH, y = n, group = OFFENSE_CODE_GROUP, color = OFFENSE_CODE_GROUP), fill=NA)+
   theme_minimal()+
   labs(title= "Frecuencia mensual de delitos",
       subtitle="Boston 2015-2019",
       caption= "Fuente: Analyze Boston",
       x= "Mes",
       y= NULL)+
  theme(text=element_text(size=12, family="Century"))+
  scale_color_brewer(palette= "Set1", name= "Tipo de delito")+
  scale_x_continuous(breaks= meses)
## Warning: Ignoring unknown parameters: fill

Observando el gráfico, pareceria ser que en septiembre hay un pico en la cantidad de delitos en relacion al resto del año. Ya que estamos analizando el hemisferio norte, esto podria deberse a la mayor actividad que se produce durante los meses calidos.

III.a. Vamos a analizar la distribución espacial de los datos a partir de un mapa de densidad que muestre donde se concentran la mayor cantidad de observaciones.

Para ello comenzamos descargando la base de OpenStreet maps, con la tematica toner-lite del estudio Stamen.

mapaBoston <- get_stamenmap(bbox, maptype = "toner-lite", zoom = 13)
## 63 tiles needed, this may take a while (try a smaller zoom).
## Source : http://tile.stamen.com/toner-lite/13/2476/3028.png
## Source : http://tile.stamen.com/toner-lite/13/2477/3028.png
## Source : http://tile.stamen.com/toner-lite/13/2478/3028.png
## Source : http://tile.stamen.com/toner-lite/13/2479/3028.png
## Source : http://tile.stamen.com/toner-lite/13/2480/3028.png
## Source : http://tile.stamen.com/toner-lite/13/2481/3028.png
## Source : http://tile.stamen.com/toner-lite/13/2482/3028.png
## Source : http://tile.stamen.com/toner-lite/13/2483/3028.png
## Source : http://tile.stamen.com/toner-lite/13/2484/3028.png
## Source : http://tile.stamen.com/toner-lite/13/2476/3029.png
## Source : http://tile.stamen.com/toner-lite/13/2477/3029.png
## Source : http://tile.stamen.com/toner-lite/13/2478/3029.png
## Source : http://tile.stamen.com/toner-lite/13/2479/3029.png
## Source : http://tile.stamen.com/toner-lite/13/2480/3029.png
## Source : http://tile.stamen.com/toner-lite/13/2481/3029.png
## Source : http://tile.stamen.com/toner-lite/13/2482/3029.png
## Source : http://tile.stamen.com/toner-lite/13/2483/3029.png
## Source : http://tile.stamen.com/toner-lite/13/2484/3029.png
## Source : http://tile.stamen.com/toner-lite/13/2476/3030.png
## Source : http://tile.stamen.com/toner-lite/13/2477/3030.png
## Source : http://tile.stamen.com/toner-lite/13/2478/3030.png
## Source : http://tile.stamen.com/toner-lite/13/2479/3030.png
## Source : http://tile.stamen.com/toner-lite/13/2480/3030.png
## Source : http://tile.stamen.com/toner-lite/13/2481/3030.png
## Source : http://tile.stamen.com/toner-lite/13/2482/3030.png
## Source : http://tile.stamen.com/toner-lite/13/2483/3030.png
## Source : http://tile.stamen.com/toner-lite/13/2484/3030.png
## Source : http://tile.stamen.com/toner-lite/13/2476/3031.png
## Source : http://tile.stamen.com/toner-lite/13/2477/3031.png
## Source : http://tile.stamen.com/toner-lite/13/2478/3031.png
## Source : http://tile.stamen.com/toner-lite/13/2479/3031.png
## Source : http://tile.stamen.com/toner-lite/13/2480/3031.png
## Source : http://tile.stamen.com/toner-lite/13/2481/3031.png
## Source : http://tile.stamen.com/toner-lite/13/2482/3031.png
## Source : http://tile.stamen.com/toner-lite/13/2483/3031.png
## Source : http://tile.stamen.com/toner-lite/13/2484/3031.png
## Source : http://tile.stamen.com/toner-lite/13/2476/3032.png
## Source : http://tile.stamen.com/toner-lite/13/2477/3032.png
## Source : http://tile.stamen.com/toner-lite/13/2478/3032.png
## Source : http://tile.stamen.com/toner-lite/13/2479/3032.png
## Source : http://tile.stamen.com/toner-lite/13/2480/3032.png
## Source : http://tile.stamen.com/toner-lite/13/2481/3032.png
## Source : http://tile.stamen.com/toner-lite/13/2482/3032.png
## Source : http://tile.stamen.com/toner-lite/13/2483/3032.png
## Source : http://tile.stamen.com/toner-lite/13/2484/3032.png
## Source : http://tile.stamen.com/toner-lite/13/2476/3033.png
## Source : http://tile.stamen.com/toner-lite/13/2477/3033.png
## Source : http://tile.stamen.com/toner-lite/13/2478/3033.png
## Source : http://tile.stamen.com/toner-lite/13/2479/3033.png
## Source : http://tile.stamen.com/toner-lite/13/2480/3033.png
## Source : http://tile.stamen.com/toner-lite/13/2481/3033.png
## Source : http://tile.stamen.com/toner-lite/13/2482/3033.png
## Source : http://tile.stamen.com/toner-lite/13/2483/3033.png
## Source : http://tile.stamen.com/toner-lite/13/2484/3033.png
## Source : http://tile.stamen.com/toner-lite/13/2476/3034.png
## Source : http://tile.stamen.com/toner-lite/13/2477/3034.png
## Source : http://tile.stamen.com/toner-lite/13/2478/3034.png
## Source : http://tile.stamen.com/toner-lite/13/2479/3034.png
## Source : http://tile.stamen.com/toner-lite/13/2480/3034.png
## Source : http://tile.stamen.com/toner-lite/13/2481/3034.png
## Source : http://tile.stamen.com/toner-lite/13/2482/3034.png
## Source : http://tile.stamen.com/toner-lite/13/2483/3034.png
## Source : http://tile.stamen.com/toner-lite/13/2484/3034.png
ggmap(mapaBoston)

Mapeamos rapidamente la data de delitos descargada de Analize Boston

delitos_geo <- delitos%>%
  filter(Long != -1,Lat !=-1)
 
ggmap(mapaBoston)+
  geom_point(data=delitos_geo, aes(x=Long, y=Lat), size=.1)

Ahora vamos a realizar una serie de mapas mostrando la densidad de delitos a nivel territorial. Comenzamos por estudiar los delitos segun la division barrial.

Vamos a asignar a los delitos la información del barrio al que pertenecen. Para esto primero debemos compatibilizar las coordenadas.

delitos_geo <- delitos_geo %>%
filter(!is.na(Long), !is.na(Lat)) %>%
st_as_sf(coords = c("Long", "Lat"), crs = 4326)

Hacemos el spatial join

delitosXbarrio <- st_join(delitos_geo, barrios)
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
delitosXbarrio <- delitosXbarrio %>% 
  group_by(Name)%>%
  summarise(total=n()) %>%
  st_set_geometry(NULL)
## Warning: Factor `Name` contains implicit NA, consider using
## `forcats::fct_explicit_na`

Acto seguido, hacemos un left join para asignarle al shape de barrios la data de delitos que le corresponde a cada poligono. Visualizamos a partir de un “Choropleth map” de delitos por barrio.

barrios <- barrios %>% left_join(delitosXbarrio)
## Joining, by = c("Name", "total")

barrios %>% 
  ggplot() + 
  geom_sf(aes(fill = total), color = NA)+
  scale_fill_distiller(palette= "YlOrRd", direction= 1)+
  geom_sf(data=barriosCent, size=.5)+
  geom_text_repel(data = barriosCent_xy, aes(X, Y, label = Name), colour = "black", size=3)+
  theme_void()+
  labs(title= "Crimenes por barrio", subtitle = "Boston, 2015-2019", x=NULL, y=NULL, caption = "Analize Boston")+
  theme(text=element_text(size=12, family="Century"))

ggplot(data=barrios, aes(x=Name, y=total))+
  geom_bar(stat="identity", fill="red")+
  coord_flip()+
  labs(title= "Crimenes por barrio", subtitle = "Boston, 2015-2019", x=NULL, y=NULL, caption = "Analize Boston")+
  theme(text=element_text(size=12, family="Century"))+
  theme_minimal()+
  theme(text=element_text(size=12, family="Century"))
## Warning: Removed 6 rows containing missing values (position_stack).

Pero tambien podemos mostrar las densidades territoriales a partir de ciertas funciones, que analizan los puntos detectando patrones de agrupacion. Con estas funciones haremos 2 tipos de mapa: de pixelado por un lado, y otro de isolineas, a partir del calculo de la densidad de Kernel.

ggmap(mapaBoston) +
    geom_bin2d(data = delitos, 
               aes(x =Long, y = Lat), bins=100)+
    scale_fill_distiller(palette="YlOrRd", direction= 1)+
  labs(title= "Crimenes por barrio", subtitle = "Boston, 2015-2019", x=NULL, y=NULL, caption = "Analize Boston")+
  theme_void()
## Warning: Removed 27780 rows containing non-finite values (stat_bin2d).

ggmap(mapaBoston) +
    geom_density2d(data = delitos, aes(x = Long, y = Lat, color = stat(level))) +
        scale_color_distiller(palette="YlOrBr", direction= 1)+
  theme(text=element_text(size=12, family="Century"))+
  labs(title= "Crimenes por barrio", subtitle = "Boston, 2015-2019", x=NULL, y=NULL, caption = "Analize Boston")+
  theme_void()
## Warning: Removed 27780 rows containing non-finite values (stat_density2d).

b. Vamos a estudiar los patrones del delito en Boston segun el tipo de crimen. Los mapas de facetado nos permiten visualizar en simultaneo distintas categorias para una misma variable, y asi compararlas mas claramente.

ggmap(mapaBoston) +
    geom_point(data = filter(delitos, OFFENSE_CODE_GROUP %in% delitos_frecuentes), 
               aes(x = Long, y = Lat, color = OFFENSE_CODE_GROUP),
               size = 0.1, alpha = 0.1) +
    scale_color_brewer(palette = "Set1", name= "Tipo de delito") +
    guides(color = guide_legend(override.aes = list(size=.1, alpha = 1))) +
    facet_wrap(~ OFFENSE_CODE_GROUP)+
    labs(title= "Crimenes por barrio segun tipo", subtitle = "Boston, 2015-2019", x=NULL, y=NULL, caption = "Analize Boston")+
    theme_void()+
    theme(text=element_text(size=12, family="Century"))
## Warning: Removed 13001 rows containing missing values (geom_point).

ggmap(mapaBoston) +
    geom_density2d(data = filter(delitos, OFFENSE_CODE_GROUP %in% delitos_frecuentes), aes(x = Long, y = Lat, color = stat(level))) +
    scale_color_distiller(palette="Reds", direction= 1) +
    facet_wrap(~OFFENSE_CODE_GROUP)+
    theme_void()+
    theme(text=element_text(size=12, family="Century"))+
    guides(color = guide_legend(override.aes = list(size=2, alpha = 1)))+
    labs(title= "Crimenes por barrio segun tipo", subtitle = "Boston, 2015-2019", x=NULL, y=NULL, caption = "Analize Boston")
## Warning: Removed 13001 rows containing non-finite values (stat_density2d).

En conclusion, a partir del facetado, observamos que, por un lado, la distribucion de los delitos es bastante homogenea. Lo que varia entre tipo de delitos es la cantidad, siendo los robos o hurtos el tipo mas usual.

Por ultimo, hacemos un analisis introduciendo la variable temporal

ggmap(mapaBoston) +
   geom_point(data = filter(delitos, 
                             OFFENSE_CODE_GROUP %in% c("Motor Vehicle Accident Response", "Medical Assistance")),
               aes(x = Long, y = Lat, color = OFFENSE_CODE_GROUP), alpha = .5, size = .1, fill=NA) +
        facet_wrap(~DAY_OF_WEEK)+
   labs(title= "Crimenes por barrio segun dia de la semana", subtitle = "Boston, 2015-2019", y=NULL, caption = "Analize Boston")+
   theme_void()+
   theme(text=element_text(size=12, family="Century"))+
   scale_color_brewer(palette="Set1", name= "Tipo de delito")+
   guides(color = guide_legend(override.aes = list(size=2, alpha = 1)))
## Warning: Removed 9570 rows containing missing values (geom_point).

Al facetar temporalmente vemos que en efecto hay una variacion de la distribucion de los reportes segun el dia de la semana.

M5 / Analizando y visualizando flujos de viajes urbanos

bluebikes <- read.csv("./00BASES/201907-bluebikes-tripdata.csv")
head(bluebikes)
##   tripduration                starttime                 stoptime
## 1          581 2019-07-01 00:00:41.0320 2019-07-01 00:10:22.7280
## 2          600 2019-07-01 00:00:49.7470 2019-07-01 00:10:50.5160
## 3          349 2019-07-01 00:05:34.4880 2019-07-01 00:11:24.4720
## 4        29455 2019-07-01 00:05:36.1540 2019-07-01 08:16:31.5770
## 5          363 2019-07-01 00:05:45.4590 2019-07-01 00:11:48.8100
## 6         1769 2019-07-01 00:06:28.5960 2019-07-01 00:35:58.0710
##   start.station.id                               start.station.name
## 1               52                        Newbury St at Hereford St
## 2               52                        Newbury St at Hereford St
## 3               91 One Kendall Square at Hampshire St / Portland St
## 4               33                                   Kenmore Square
## 5               33                                   Kenmore Square
## 6              389          Everett Square (Broadway at Chelsea St)
##   start.station.latitude start.station.longitude end.station.id
## 1               42.34872               -71.08595             58
## 2               42.34872               -71.08595             58
## 3               42.36628               -71.09169             96
## 4               42.34871               -71.09701             33
## 5               42.34871               -71.09701             33
## 6               42.40726               -71.05546             76
##                                                         end.station.name
## 1                                                 Mugar Way at Beacon St
## 2                                                 Mugar Way at Beacon St
## 3                     Cambridge Main Library at Broadway / Trowbridge St
## 4                                                         Kenmore Square
## 5                                                         Kenmore Square
## 6 Central Sq Post Office / Cambridge City Hall at Mass Ave / Pleasant St
##   end.station.latitude end.station.longitude bikeid   usertype birth.year
## 1             42.35554             -71.07287   1377   Customer       1997
## 2             42.35554             -71.07287   1254   Customer       1997
## 3             42.37338             -71.11107   3807 Subscriber       1998
## 4             42.34871             -71.09701   3815   Customer       1969
## 5             42.34871             -71.09701   4364   Customer       1969
## 6             42.36643             -71.10550   1106 Subscriber       1994
##   gender
## 1      2
## 2      2
## 3      1
## 4      0
## 5      0
## 6      1

Ploteamos

    ggplot() +
        geom_bar(data=bluebikes, aes(x=birth.year))

Al ver el grafico observamos que hay un usuario triatlonista con record de viajes que nos esta generando un outlier. Para homogeneizar la data simplemente vamos a filtrarlo, y a su vez, reacomodar el grafico para su mejor lectura.

bluebikes2 <-bluebikes %>% 
  group_by(birth.year)%>%
  summarise(total=n()) 
max(bluebikes2$total)
## [1] 40483

Vamos a la tabla y vemos que ese valor se corresponde con el año 1969

bluebikes <-bluebikes %>% 
  filter(birth.year != 1969)
max(bluebikes$birth.year)
## [1] 2003
min(bluebikes$birth.year)
## [1] 1888
years <- seq(1945,2005, by= 5)
    ggplot() +
        geom_bar(data=bluebikes, aes(x=birth.year), fill="blue")+
   coord_cartesian (xlim=c(1940,2005))+
   scale_x_continuous(breaks=years)+
   theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))

A partir del ploteo vemos que los milennials tienen la delantera en el uso del sistema de bicicletas.

Ahora si, vamos a analizar los viajes según su origen-destino. Primero comenzamos identificando los tramos que se dan con mayor frecuencia.

conteo <- bluebikes %>% 
    group_by(start.station.id, end.station.id) %>% 
    summarise(total = n())

Graficamos el conteo

ggplot() + 
    geom_tile(data = conteo, aes(x = start.station.id, y = end.station.id, fill = total)) +
    scale_fill_distiller(palette = "Spectral")+
  theme_minimal()+
   labs(title= "Crimenes por barrio segun dia de la semana", subtitle = "Boston, 2015-2019", x= "Estación inicial", y="Estación terminal", caption = "Analize Boston")+
    theme(text=element_text(size=12, family="Century"))

A partir del mapa de calor observamos que la frecuencia se concentra entre ciertas estaciones, y que hay algunas que directamente, nunca fueron conectadas por un recorrido. Determinadas estaciones tampoco parecerian ser usadas para viajes circulares. Las franjas vacias indicarian tambien que la numeracion no es continua o hay estaciones en desuso.

Para hacer mas clara la visualizacion vamos a filtrar el top 10 de viajes, a su vez eliminando los “viajes circulares”.

top10 <- conteo %>% 
    ungroup() %>% 
    filter(start.station.id!= end.station.id) %>% 
    top_n(10)
## Selecting by total

top10
## # A tibble: 10 x 3
##    start.station.id end.station.id total
##               <int>          <int> <int>
##  1               67            179   369
##  2               80            178   467
##  3              100            118   414
##  4              118            100   423
##  5              178             80   365
##  6              178            107   423
##  7              179             67   364
##  8              179            107   384
##  9              239            328   450
## 10              328            239   349
ggplot() + 
    geom_tile(data = top10, 
              aes(x = as.factor(start.station.id),
                  y = as.factor(end.station.id),
                  fill = total)) +
    scale_fill_distiller(palette = "Spectral")+
    labs(title= "Heatmap delitos", subtitle = "Boston, 2015-2019", x= "Estación inicial", y="Estación terminal", caption = "Analize Boston")+
    theme(text=element_text(size=12, family="Century"))+
    theme_minimal()

Observamos que el viaje mas popular se da entre las estaciones 328 y 239, seguido de 178 y 328.

Vamos a mapear las rutas Hacemos un left join para asignar las coordinadas al filtro top 10.

top10 <- left_join(top10, bluebikes)
## Joining, by = c("start.station.id", "end.station.id")
str(top10)
## Classes 'tbl_df', 'tbl' and 'data.frame':    4008 obs. of  16 variables:
##  $ start.station.id       : int  67 67 67 67 67 67 67 67 67 67 ...
##  $ end.station.id         : int  179 179 179 179 179 179 179 179 179 179 ...
##  $ total                  : int  369 369 369 369 369 369 369 369 369 369 ...
##  $ tripduration           : int  373 252 333 321 313 265 291 354 675 282 ...
##  $ starttime              : Factor w/ 316898 levels "2019-07-01 00:00:41.0320",..: 4165 4488 4568 4580 4687 5911 6034 7281 7373 7930 ...
##  $ stoptime               : Factor w/ 316899 levels "2019-07-01 00:10:22.7280",..: 4055 4377 4475 4493 4584 5701 5829 7025 7248 7663 ...
##  $ start.station.name     : Factor w/ 302 levels "175 N Harvard St",..: 192 192 192 192 192 192 192 192 192 192 ...
##  $ start.station.latitude : num  42.4 42.4 42.4 42.4 42.4 ...
##  $ start.station.longitude: num  -71.1 -71.1 -71.1 -71.1 -71.1 ...
##  $ end.station.name       : Factor w/ 303 levels "175 N Harvard St",..: 196 196 196 196 196 196 196 196 196 196 ...
##  $ end.station.latitude   : num  42.4 42.4 42.4 42.4 42.4 ...
##  $ end.station.longitude  : num  -71.1 -71.1 -71.1 -71.1 -71.1 ...
##  $ bikeid                 : int  4439 41 2572 3707 3421 172 2851 3520 2497 1573 ...
##  $ usertype               : Factor w/ 2 levels "Customer","Subscriber": 2 2 2 2 2 2 2 2 2 2 ...
##  $ birth.year             : int  1960 1987 1981 1982 1992 1991 1990 1992 1998 1996 ...
##  $ gender                 : int  1 1 1 1 1 1 1 1 1 1 ...

Ruteamos…

library(osrm)
## Warning: package 'osrm' was built under R version 3.6.1
## Data: (c) OpenStreetMap contributors, ODbL 1.0 - http://www.openstreetmap.org/copyright
## Routing: OSRM - http://project-osrm.org/
obtener_recorrido <- function(o_nombre, o_x, o_y, d_nombre, d_x, d_y) {
    
    ruta <- osrmRoute(src = c(o_nombre, o_x, o_y),
                      dst = c(d_nombre, d_x, d_y))
    
    cbind(start.station.id = o_nombre, end.station.id = d_nombre, ruta)
    
}


argumentos <- list(top10$start.station.id, top10$start.station.longitude, top10$start.station.latitude,
                  top10$end.station.id, top10$end.station.longitude, top10$end.station.latitude)
recorridos <- pmap_df(argumentos, obtener_recorrido)
recorridos <- recorridos %>% 
    mutate(ID = paste(start.station.id, "-",end.station.id))
ggmap(mapaBoston) +
    geom_path(data = recorridos, aes(x = lon, y = lat, color = ID)) +
    theme_nothing()
## Warning: Removed 4617 rows containing missing values (geom_path).

Los recorridos populares se producen en el area mas centrica en su mayoria. Vamos a tener que generar un nuevo mapa base para mejorar la visualizacion.

bbox <- make_bbox(recorridos$lon, recorridos$lat)

bbox
##      left    bottom     right       top 
## -71.13345  42.35339 -71.08601  42.40332
mapaRecorridos <- get_stamenmap(bbox, color = "bw", zoom = 13)
## Source : http://tile.stamen.com/terrain/13/2477/3028.png
## Source : http://tile.stamen.com/terrain/13/2478/3028.png
## Source : http://tile.stamen.com/terrain/13/2477/3029.png
## Source : http://tile.stamen.com/terrain/13/2478/3029.png
## Source : http://tile.stamen.com/terrain/13/2477/3030.png
## Source : http://tile.stamen.com/terrain/13/2478/3030.png
recorridosGroup <- recorridos %>% 
    group_by(ID) %>% 
    summarise(total = n())
recorridos_geo <- left_join(recorridosGroup, recorridos)
## Joining, by = "ID"
ggmap(mapaRecorridos) +
    geom_path(data = recorridos_geo, aes(x = lon, y = lat, color = ID, size=total), alpha = 0.7) +
    labs(title= "Top 10 recorridos en bicicleta", subtitle = "Boston, 2019", y=NULL, caption = "Analize Boston")+
    theme_void()+
    theme(text=element_text(size=12, family="Century"))