Comenzamos cargando las librerías a utilizar

Parte A ( Datos de la plataforma de datos abiertos)

library(tidyverse)
## -- Attaching packages ------------------------------------------------------------------------------ tidyverse 1.2.1 --
## v ggplot2 3.2.0       v purrr   0.3.2  
## v tibble  2.1.1       v dplyr   0.8.0.1
## v tidyr   0.8.3       v stringr 1.4.0  
## v readr   1.3.1       v forcats 0.4.0
## -- Conflicts --------------------------------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(ggplot2)
library(sf)
## Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
library(rgdal)
## Loading required package: sp
## 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/Elisa/Documents/R/win-library/3.5/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/Elisa/Documents/R/win-library/3.5/rgdal/proj
##  Linking to sp version: 1.3-1
require(PBSmapping)
## Loading required package: PBSmapping
## 
## -----------------------------------------------------------
## 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/Elisa/Documents/R/win-library/3.5/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
## Checking rgeos availability: TRUE
library(ggrepel)

PARTE B ( Datos de OpenStreet Map)

library(leaflet)
library(osmdata)
## Data (c) OpenStreetMap contributors, ODbL 1.0. http://www.openstreetmap.org/copyright
library(tidyverse)
  1. MAPEANDO LOS DATOS OFICIALES

Cargamos los shapes a unir.

barrios <- st_read("Boston_Neighborhoods.geojson")
## Reading layer `Boston_Neighborhoods' from data source `C:\Users\Elisa\Google Drive\[B]ACADEMICA\FORMACION\[02]MEU\CIENCIA DE DATOS II - A BRUST\C2\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 ("Public_Schools.geojson")
## Reading layer `Public_Schools' from data source `C:\Users\Elisa\Google Drive\[B]ACADEMICA\FORMACION\[02]MEU\CIENCIA DE DATOS II - A BRUST\C2\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(escuelasXbarrio)%>%
  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)
## Joining, by = "Name"

  1. MAPEANDO DATOS DE OPENSETREET MAPS

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

Obtehemos 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':   20790 obs. of  268 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 ...
##  $ 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_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.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 ...
##  $ goods                                : chr  NA NA NA NA ...
##  $ handrail                             : chr  NA NA NA NA ...
##  $ handrail.center                      : chr  NA NA NA NA ...
##  $ hazmat                               : chr  NA NA NA NA ...
##  $ height                               : chr  NA NA NA NA ...
##  $ hgv                                  : chr  NA NA NA NA ...
##  $ hgv.conditional                      : chr  NA NA NA NA ...
##  $ highway                              : chr  "footway" "footway" "primary" "secondary" ...
##   [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 “highway”, 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.

colores <-c("good"="light grey","fair" ="grey", "fair;deficient"="#fecc5c","deficient" = "#fd8d3c", "poor" = "#f03b20","intolerable" = "#bd0026" )

ggplot()+
  geom_sf(data=calles, aes(color=condition), alpha=1)+
  scale_color_manual(values=colores, aesthetics = c("color"))+
  theme_void()+
  labs(title= "Boston",
       subtitle="Estado de las calles",
       caption= "Fuente: OpenStreet Map")

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)

II.a. Al mapear las escuelas parececia ser que su distribucion es homogenea a lo largo de todo el shape de barrios. Sin embargo, es dudosa la calidad de la informacion, por lo que habria que identificar con mayor exacittud que representan los puntos( que se agrupan como formando poligonos en vez de ser aislados)

Para hacer un conteo de las escuelas por barrio, utilizaremos 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 = "Name"
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)

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 a la mala calidad de los datos, que deberan ser filtrados para nuevamente,hacer una comparativa mas fidedigna y poder sacar conclusiones.