TAREA 1: CALCULANDO Y MAPEANDO AGREGADOS POR ÁREA + ADQUIRIENDO OPEN DATA URBANA.

Paso I y II: seleccionamos una ciudad Boston con un dataset con registros geo-referenciados con sus coordenadas: Bibliotecas Públicas.

Cargamos las librerias que necesitamos para trabajar.

library(tidyverse)
## ── Attaching packages ────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.0     ✓ purrr   0.3.3
## ✓ tibble  2.1.3     ✓ dplyr   0.8.5
## ✓ tidyr   1.0.2     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.5.0
## ── Conflicts ───────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(sf)
## Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0

Ahora vamos a leer los dataset seleccionados:

Bibliotecas Públicas utilizando el comando st_read () y lo leeremos directamente de la web:

Bibliotecas_Publicas_OK <- st_read("http://bostonopendata-boston.opendata.arcgis.com/datasets/cb00f9248aa6404ab741071ca3806c0e_6.geojson?outSR={%22latestWkid%22:2249,%22wkid%22:102686}")
## Reading layer `bbfc2420-19c6-49d8-ac42-015431c7cce2202042-1-18eeelv.jwzp' from data source `http://bostonopendata-boston.opendata.arcgis.com/datasets/cb00f9248aa6404ab741071ca3806c0e_6.geojson?outSR={%22latestWkid%22:2249,%22wkid%22:102686}' using driver `GeoJSON'
## Simple feature collection with 26 features and 7 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -71.16788 ymin: 42.2571 xmax: -71.02858 ymax: 42.37777
## CRS:            4326

Filtramos los datos que nos interesan del dataset para que quede más ordenado:

Bibliotecas_Publicas_Ordenadas <- Bibliotecas_Publicas_OK %>% 
  select(LIBRARIES_, DISTRICT, ST_ADDRESS, BRANCH, geometry)

head(Bibliotecas_Publicas_Ordenadas)
## Simple feature collection with 6 features and 4 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -71.12214 ymin: 42.2571 xmax: -71.06501 ymax: 42.36145
## CRS:            4326
##   LIBRARIES_      DISTRICT        ST_ADDRESS        BRANCH
## 1          5                151 Cambridge St      West End
## 2         11               690 Washington St Codman Square
## 3         17     HYDE PARK     35 Harvard Av     Hyde Park
## 4         18 JAMAICA PLAIN     433 Centre St      Connolly
## 5         22       ROXBURY     149 Dudley St        Dudley
## 6          1 COPLEY SQUARE   700 Boylston St       Central
##                     geometry
## 1 POINT (-71.06501 42.36145)
## 2 POINT (-71.07097 42.28752)
## 3  POINT (-71.12214 42.2571)
## 4  POINT (-71.1111 42.32065)
## 5 POINT (-71.08385 42.32849)
## 6 POINT (-71.07887 42.34944)

Ahora leemos el dataset de Barrios de Boston, utilizando el comando st_read ya que es un elemento geográfico.

Barrios_Boston_OK <- st_read("Boston_Neighborhoods.shp")
## Reading layer `Boston_Neighborhoods' from data source `/Users/ani/Desktop/UTDT - materias/04 - 0118 - MU119 Ciencia De Datos Para Ciudades II/Boston/Ciencia de Datos de Ciudades II/Boston_Neighborhoods.shp' using driver `ESRI Shapefile'
## Simple feature collection with 26 features and 7 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 739715.3 ymin: 2908294 xmax: 812255.1 ymax: 2970086
## CRS:            2249
head(Barrios_Boston_OK)
## Simple feature collection with 6 features and 7 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 751472.1 ymin: 2922750 xmax: 776215.4 ymax: 2953742
## CRS:            2249
##   OBJECTID             Name      Acres Neighborho SqMiles  ShapeSTAre
## 1       27       Roslindale 1605.56824         15    2.51  69938272.9
## 2       28    Jamaica Plain 2519.24539         11    3.94 109737890.8
## 3       29     Mission Hill  350.85356         13    0.55  15283120.0
## 4       30         Longwood  188.61195         28    0.29   8215903.5
## 5       31      Bay Village   26.53984         33    0.04   1156070.8
## 6       32 Leather District   15.63991         27    0.02    681271.7
##   ShapeSTLen                       geometry
## 1  53563.913 MULTIPOLYGON (((757409.1 29...
## 2  56349.937 MULTIPOLYGON (((762983.8 29...
## 3  17918.724 MULTIPOLYGON (((766903.6 29...
## 4  11908.757 MULTIPOLYGON (((764826.9 29...
## 5   4650.635 MULTIPOLYGON (((773315.7 29...
## 6   3237.141 MULTIPOLYGON (((775544.1 29...
names(Barrios_Boston_OK)
## [1] "OBJECTID"   "Name"       "Acres"      "Neighborho" "SqMiles"   
## [6] "ShapeSTAre" "ShapeSTLen" "geometry"

Vamos a explorar el dataset Barrios_Boston, utilizando el comando plot, para tener una rápida observación de las 7 variables del dataset:

plot(Barrios_Boston_OK)

Observamos las 7 variables de nuestro dataset de forma correcta, algunas más relevantes que otras.

Ahora cruzamos nuestros datos de Barrios y Bibliotecas Públicas en un mapa a ver que observamos:

ggplot() +
  geom_sf(data = Barrios_Boston_OK) +
  geom_sf(data = Bibliotecas_Publicas_Ordenadas,
          color = "red")

Observamos las bibliotecas asignadas a los distintos barrios.

Para poder continnuar debemos verificar si los sistemas de coordenadas de ambos datasets coinciden, en caso contrario debemos transformarlo a EPSG:4326 con el comando st_transform().

st_crs(Barrios_Boston_OK)
## Coordinate Reference System:
##   User input: 2249 
##   wkt:
## PROJCS["NAD83 / Massachusetts Mainland (ftUS)",
##     GEOGCS["NAD83",
##         DATUM["North_American_Datum_1983",
##             SPHEROID["GRS 1980",6378137,298.257222101,
##                 AUTHORITY["EPSG","7019"]],
##             TOWGS84[0,0,0,0,0,0,0],
##             AUTHORITY["EPSG","6269"]],
##         PRIMEM["Greenwich",0,
##             AUTHORITY["EPSG","8901"]],
##         UNIT["degree",0.0174532925199433,
##             AUTHORITY["EPSG","9122"]],
##         AUTHORITY["EPSG","4269"]],
##     PROJECTION["Lambert_Conformal_Conic_2SP"],
##     PARAMETER["standard_parallel_1",42.68333333333333],
##     PARAMETER["standard_parallel_2",41.71666666666667],
##     PARAMETER["latitude_of_origin",41],
##     PARAMETER["central_meridian",-71.5],
##     PARAMETER["false_easting",656166.667],
##     PARAMETER["false_northing",2460625],
##     UNIT["US survey foot",0.3048006096012192,
##         AUTHORITY["EPSG","9003"]],
##     AXIS["X",EAST],
##     AXIS["Y",NORTH],
##     AUTHORITY["EPSG","2249"]]
st_crs(Bibliotecas_Publicas_Ordenadas)
## 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"]]

El dataset Barrios_Boston_OK deber ser modificado:

Barrios_Boston_OK <- st_transform(Barrios_Boston_OK, 4326)

Para verificar si dicho dataset se transformo correctamente al sistema de coordenadas de Proyeccion Mercator (4326) volvemos a observarlo:

st_crs(Barrios_Boston_OK)
## Coordinate Reference System:
##   User input: EPSG: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"]]

Paso III: Ahora mediante un join espacial, utilizando el comando st_join(). A cada biblioteca le asignaremos el barrio que le corresponda.

Bibliotecas_Publicas_Ordenadas_con_Barrios <- st_join(Bibliotecas_Publicas_Ordenadas, Barrios_Boston_OK)
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
head(Bibliotecas_Publicas_Ordenadas_con_Barrios)
## Simple feature collection with 6 features and 11 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -71.12214 ymin: 42.2571 xmax: -71.06501 ymax: 42.36145
## CRS:            4326
##   LIBRARIES_      DISTRICT        ST_ADDRESS        BRANCH OBJECTID
## 1          5                151 Cambridge St      West End       40
## 2         11               690 Washington St Codman Square       48
## 3         17     HYDE PARK     35 Harvard Av     Hyde Park       46
## 4         18 JAMAICA PLAIN     433 Centre St      Connolly       28
## 5         22       ROXBURY     149 Dudley St        Dudley       35
## 6          1 COPLEY SQUARE   700 Boylston St       Central       37
##            Name     Acres Neighborho SqMiles ShapeSTAre ShapeSTLen
## 1      West End  190.4907         31    0.30    8297743   17728.59
## 2    Dorchester 4662.8795          6    7.29  203114217  104344.03
## 3     Hyde Park 2927.2212         10    4.57  127509243   66861.24
## 4 Jamaica Plain 2519.2454         11    3.94  109737891   56349.94
## 5       Roxbury 2108.4691         16    3.29   91844546   49488.80
## 6      Back Bay  399.3144          2    0.62   17394066   19455.67
##                     geometry
## 1 POINT (-71.06501 42.36145)
## 2 POINT (-71.07097 42.28752)
## 3  POINT (-71.12214 42.2571)
## 4  POINT (-71.1111 42.32065)
## 5 POINT (-71.08385 42.32849)
## 6 POINT (-71.07887 42.34944)
names(Bibliotecas_Publicas_Ordenadas_con_Barrios)
##  [1] "LIBRARIES_" "DISTRICT"   "ST_ADDRESS" "BRANCH"     "OBJECTID"  
##  [6] "Name"       "Acres"      "Neighborho" "SqMiles"    "ShapeSTAre"
## [11] "ShapeSTLen" "geometry"

Con este nuevo dataset, estamos en condiciones de saber, cuantas Bibliotecas Públias hay por Barrio = Neighborho. Asimismo, podremos graficar mapas de valores máximos, mínimos, promedios, entre otros, de la cantidad de Bibliotecas que haya por Barrio.

Paso IV:

Vamos a graficar un gráfico de puntos para visualizar la cantidad de bilioteas por barrio. Primero debemos agrupar los bibliotecas por barrio = Neighborhood_ID.

Total_Bibliotecas_Publicas_Ordenadas_con_Barrios <- Bibliotecas_Publicas_Ordenadas_con_Barrios %>% 
  group_by(Neighborho) %>% 
  summarise(total=n())

Total_Bibliotecas_Publicas_Ordenadas_con_Barrios
## Simple feature collection with 18 features and 2 fields
## geometry type:  GEOMETRY
## dimension:      XY
## bbox:           xmin: -71.16788 ymin: 42.2571 xmax: -71.02858 ymax: 42.37777
## CRS:            4326
## # A tibble: 18 x 3
##    Neighborho total                                                     geometry
##  * <fct>      <int>                                               <GEOMETRY [°]>
##  1 10             1                                    POINT (-71.12214 42.2571)
##  2 11             2       MULTIPOINT ((-71.11488 42.30863), (-71.1111 42.32065))
##  3 12             2      MULTIPOINT ((-71.09314 42.27752), (-71.06851 42.27355))
##  4 13             1                                   POINT (-71.09908 42.33246)
##  5 14             1                                   POINT (-71.05502 42.36408)
##  6 15             1                                   POINT (-71.12865 42.28565)
##  7 16             2      MULTIPOINT ((-71.09563 42.31406), (-71.08385 42.32849))
##  8 17             1                                   POINT (-71.03872 42.33584)
##  9 19             1                                    POINT (-71.15745 42.2834)
## 10 2              1                                   POINT (-71.07887 42.34944)
## 11 24             1                                   POINT (-71.12811 42.36011)
## 12 25             2      MULTIPOINT ((-71.16788 42.35132), (-71.15291 42.34769))
## 13 31             1                                   POINT (-71.06501 42.36145)
## 14 32             1                                   POINT (-71.07703 42.34133)
## 15 4              1                                   POINT (-71.06444 42.37582)
## 16 6              5 MULTIPOINT ((-71.08119 42.30799), (-71.07097 42.28752), (-7…
## 17 7              1                                   POINT (-71.06306 42.35214)
## 18 8              1                                   POINT (-71.02858 42.37777)
ggplot(Total_Bibliotecas_Publicas_Ordenadas_con_Barrios) +
  geom_point(aes(x=Neighborho, y=total), color="red") + 
  labs(title="Bibliotecas por Barrio",
       subtitle="Boston",
       caption="Fuente:https://data.boston.gov/dataset/public-libraries/resource/ed710c4a-689a-47af-9dd0-6e4215003c24",
       x="Barrio ID",
       y="Total")

Se observa que en los barrios prevalecen 1 Biblioteca Pública por Barrio aunque en 4 Barrios se destacan 2 Bibliotecas Púbicas. Es llamativo en el Barrio ID: 6 la concentración de 5 Bibliotecas Públicas, habrá que estudiar a que se debe esto.

Mapeamos las biliotecas por barrio distinguiendo por color:

ggplot() +
  geom_sf(data = Barrios_Boston_OK) +
  geom_sf(data = Total_Bibliotecas_Publicas_Ordenadas_con_Barrios, aes(color = Neighborho), size=1) +
  theme(legend.position = "bottom")

Ahora vamos a pintar cada barrio = Neighborho con la cantidad máxima de Bibliotecas Públicas, utilizando un gráfico de coropletas.

Para ello utilizamos al dataset Total_Bibliotecas_Publicas_Ordenadas_con_Barrios y le quitamos la geometría. Luego lo unimos al dataset Barrios_Boston_OK:

Total_Bibliotecas_Publicas_Ordenadas_con_Barrios <- Total_Bibliotecas_Publicas_Ordenadas_con_Barrios %>% 
  st_set_geometry(NULL)

head(Total_Bibliotecas_Publicas_Ordenadas_con_Barrios)
## # A tibble: 6 x 2
##   Neighborho total
##   <fct>      <int>
## 1 10             1
## 2 11             2
## 3 12             2
## 4 13             1
## 5 14             1
## 6 15             1

Modificamos el nombre de la columa Total = Total Bibliotecas Publicas

colnames(Total_Bibliotecas_Publicas_Ordenadas_con_Barrios)[colnames(Total_Bibliotecas_Publicas_Ordenadas_con_Barrios) %in% c('total')] <- paste0(c('Total_Bibliotecas_Publicas'))

head(Total_Bibliotecas_Publicas_Ordenadas_con_Barrios)
## # A tibble: 6 x 2
##   Neighborho Total_Bibliotecas_Publicas
##   <fct>                           <int>
## 1 10                                  1
## 2 11                                  2
## 3 12                                  2
## 4 13                                  1
## 5 14                                  1
## 6 15                                  1

Mediante el comando left_join() unificamos el datasat Total_Bibliotecas_Publicas_Ordenadas_con_Barrios con Barrios_Boston_OK:

Bibliotecas_Barrios_Total <- Barrios_Boston_OK %>% 
  left_join(Total_Bibliotecas_Publicas_Ordenadas_con_Barrios, by="Neighborho")

head(Bibliotecas_Barrios_Total)
## Simple feature collection with 6 features and 8 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -71.14779 ymin: 42.26761 xmax: -71.05588 ymax: 42.35237
## CRS:            EPSG:4326
##   OBJECTID             Name      Acres Neighborho SqMiles  ShapeSTAre
## 1       27       Roslindale 1605.56824         15    2.51  69938272.9
## 2       28    Jamaica Plain 2519.24539         11    3.94 109737890.8
## 3       29     Mission Hill  350.85356         13    0.55  15283120.0
## 4       30         Longwood  188.61195         28    0.29   8215903.5
## 5       31      Bay Village   26.53984         33    0.04   1156070.8
## 6       32 Leather District   15.63991         27    0.02    681271.7
##   ShapeSTLen Total_Bibliotecas_Publicas                       geometry
## 1  53563.913                          1 MULTIPOLYGON (((-71.12593 4...
## 2  56349.937                          2 MULTIPOLYGON (((-71.10499 4...
## 3  17918.724                          1 MULTIPOLYGON (((-71.09043 4...
## 4  11908.757                         NA MULTIPOLYGON (((-71.09811 4...
## 5   4650.635                         NA MULTIPOLYGON (((-71.06663 4...
## 6   3237.141                         NA MULTIPOLYGON (((-71.05838 4...
names(Bibliotecas_Barrios_Total)
## [1] "OBJECTID"                   "Name"                      
## [3] "Acres"                      "Neighborho"                
## [5] "SqMiles"                    "ShapeSTAre"                
## [7] "ShapeSTLen"                 "Total_Bibliotecas_Publicas"
## [9] "geometry"

Finalmente mapeamos:

ggplot() +
  geom_sf(data = Bibliotecas_Barrios_Total, aes(fill = Total_Bibliotecas_Publicas)) + 
  geom_sf_text(data = Bibliotecas_Barrios_Total, aes(label = Name), size = 2, colour = "black") +
  labs(title = "Total Bibliotecas Públicas por Barrio",
         subtitle = "Boston",
         fill = " Total Bibliotecas Públicas",
         caption= "Fuente: https://data.boston.gov/dataset/boston-neighborhoods y
                   https://data.boston.gov/dataset/public-libraries/resource/ed710c4a-689a-47af-9dd0-6e4215003c24",
         y="",
         x="") +
  scale_fill_gradient(low="yellow", high="red")
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data
## Warning: Computation failed in `stat_sf_coordinates()`:
## Evaluation error: TopologyException: Input geom 1 is invalid: Ring Self-intersection at or near point -71.125747341600402 42.27233853935423 at -71.125747341600402 42.27233853935423.

Ahora realizamos un mapa de coropletas contemplando la densidad de bibliotecas por barrio, para verificar si los barrios más grande tienen la mayor cantidad de bibliotecas:

ggplot() +
  geom_sf(data = Bibliotecas_Barrios_Total, aes(fill = Total_Bibliotecas_Publicas/Acres)) + 
  geom_sf_text(data = Bibliotecas_Barrios_Total, aes(label = Name), size = 2, colour = "black") +
  labs(title = "Total Bibliotecas Públicas por Barrio",
         subtitle = "Boston",
         fill = " Total Bibliotecas Públicas",
         caption= "Fuente: https://data.boston.gov/dataset/boston-neighborhoods y
                   https://data.boston.gov/dataset/public-libraries/resource/ed710c4a-689a-47af-9dd0-6e4215003c24",
         y="",
         x="") +
  scale_fill_gradient(low="yellow", high="red")
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data
## Warning: Computation failed in `stat_sf_coordinates()`:
## Evaluation error: TopologyException: Input geom 1 is invalid: Ring Self-intersection at or near point -71.125747341600402 42.27233853935423 at -71.125747341600402 42.27233853935423.

Contemplamos que la distribución de bibliotecas es prácticamente unifome. Llama la atención el Barrio ID 14, de menor tamaño, posee la mayor cantidad, al igual que sus barrios circundantes. Esto se debe a que al poseer 1 sola Biblioteca Pública y tener una superficie tan pequeña supera ampliamente en densidad al Barrio ID 6 de mayor tamaño que queda coloreado con poca densidad.

TAREA 2: EXPLORANDO Y MAPEANDO INFORMACIÓN GERREFERENCIADA DE OpenStreetMap

En primer lugar debemos instalar el paquete osmdata que nos permitira acceder a los datos de OpenStreetMap y el paquete leaflet que nos posiilitará generar napas intercativos de forma rápida. Luego cargamos sus librerias correspondientes.

library(osmdata)
## Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright
library(leaflet)

Paso I: Debemos descargar de OpenStreetMap la grilla de las calles de la Ciudad de Boston. Para ello primero debemos delimiar la ciudad.

bbox_Boston <- getbb("Boston, United States of America")
bbox_poly_Boston <- getbb("Boston, United States of America", format_out = "sf_polygon")

Para verificar que estamos viendo utilizamos el comando leaflet():

leaflet(bbox_poly_Boston) %>% 
  addTiles() %>% 
  addPolygons()

Efectivamente se corresponde con la Ciudad de Boston que veníamos trabajando.

Ahora debemos decargar la información, utilizando el comando osmdata_sf().

Vamos a descargar las Calles = Highway y luego vamos a mapearla según su atributo de Velocidad.

Boston_Calles <- opq(bbox_Boston) %>% 
    add_osm_feature(key = "highway") %>% 
    osmdata_sf()

Recortamos las calles que exceden nuestra ciudad, utilizando el comando st_intersection

Boston_Calles <- st_intersection(Boston_Calles$osm_lines, bbox_poly_Boston)
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
Boston_Calles <- Boston_Calles %>%
  mutate(maxspeed = as.numeric(str_remove(maxspeed, "mph")),
         lanes = ifelse(is.na(lanes), 1, as.numeric(lanes)))
## Warning: NAs introduced by coercion

Ahora graficamos:

ggplot(Boston_Calles) +
    geom_sf(aes(color = maxspeed), alpha = 0.5) +
    scale_color_viridis_c() +
      theme_void() +
    labs(title = "Boston",
         subtitle = "Vías de Circulación - Velocidad",
         caption = "Fuente: OpenStreetMap",
         color = "velocidad máxima")

Paso II: Ahora descargaremos una capa de datos de tipo puntos, serán las Bibliotecas = Library, que no son solo las públicas, sino que incluye todas las bibliotecas donde se pueden tomar prestados libros en la ciudad.

Boston_Biliotecas_OSM <- opq(bbox_Boston) %>% 
  add_osm_feature(key = "amenity", value = "library") %>% 
  osmdata_sf()

Y también eliminamos las bibliotecas que estan afuera de la ciudad:

Boston_Biliotecas_OSM <- st_intersection(Boston_Biliotecas_OSM$osm_points, bbox_poly_Boston)
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries

Paso II-A: Mapeamos para ver la distibución general de Bibliotecas en la Ciudad:

ggplot() +
  geom_sf(data = Boston_Calles, 
            color = "gray48") +
  geom_sf(data = Boston_Biliotecas_OSM, 
            color = "red", alpha = .1) +
  labs(title = "Bibliotecas",
       subtitle = "Boston",
       caption = "Fuente: OpenStreetMap")

Paso II-A: Ahora agruparemos las biliotecas por barrio y utilizando un mapa de coropletas veremos como se distibuyen la mayor cantidad de biliotecas en la ciudad.

Para ello primero realizamos un join espacial de los Barrios de Boston con las Bibliotecas OSM utilizando el comando st_join:

Bibliotecas_OSM_en_Barrios <- st_join(Boston_Biliotecas_OSM, Barrios_Boston_OK)
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
names(Bibliotecas_OSM_en_Barrios)
##  [1] "osm_id"                 "name"                   "addr.city"             
##  [4] "addr.housenumber"       "addr.place"             "addr.postcode"         
##  [7] "addr.state"             "addr.street"            "addr.unit"             
## [10] "alt_name"               "amenity"                "attribution"           
## [13] "building"               "created_by"             "designation"           
## [16] "ele"                    "email"                  "entrance"              
## [19] "gnis.county_name"       "gnis.feature_id"        "gnis.import_uuid"      
## [22] "gnis.reviewed"          "is_in"                  "level"                 
## [25] "library"                "note"                   "office"                
## [28] "opening_hours"          "operator"               "operator.type"         
## [31] "phone"                  "postal_code"            "ref.isil"              
## [34] "region"                 "source"                 "source.url"            
## [37] "source_url"             "toilets"                "toilets.wheelchair"    
## [40] "type"                   "website"                "wheelchair"            
## [43] "wheelchair.description" "wifi"                   "wikidata"              
## [46] "wikipedia"              "OBJECTID"               "Name"                  
## [49] "Acres"                  "Neighborho"             "SqMiles"               
## [52] "ShapeSTAre"             "ShapeSTLen"             "geometry"

Ahora agrupamos por barrio:

Total_Bibliotecas_OSM_en_Barrios <- Bibliotecas_OSM_en_Barrios %>% 
  group_by(Neighborho) %>% 
  summarise(total=n())

Total_Bibliotecas_OSM_en_Barrios
## Simple feature collection with 23 features and 2 fields
## geometry type:  GEOMETRY
## dimension:      XY
## bbox:           xmin: -71.17893 ymin: 42.23389 xmax: -71.02771 ymax: 42.38127
## CRS:            EPSG:4326
## # A tibble: 23 x 3
##    Neighborho total                                                     geometry
##  * <fct>      <int>                                               <GEOMETRY [°]>
##  1 10             2      MULTIPOINT ((-71.13317 42.23389), (-71.11774 42.26296))
##  2 11            18 MULTIPOINT ((-71.12824 42.30121), (-71.12094 42.30745), (-7…
##  3 12            20 MULTIPOINT ((-71.09324 42.27703), (-71.09323 42.27745), (-7…
##  4 13             2       MULTIPOINT ((-71.10667 42.33007), (-71.09556 42.3356))
##  5 14             3 MULTIPOINT ((-71.05499 42.36401), (-71.05477 42.36649), (-7…
##  6 16            28 MULTIPOINT ((-71.0919 42.33246), (-71.08593 42.3168), (-71.…
##  7 17             1                                   POINT (-71.07094 42.33292)
##  8 19             4 MULTIPOINT ((-71.17893 42.27802), (-71.17454 42.28218), (-7…
##  9 2             53 MULTIPOINT ((-71.08714 42.34664), (-71.08616 42.34512), (-7…
## 10 24            52 MULTIPOINT ((-71.12839 42.36023), (-71.12825 42.36033), (-7…
## # … with 13 more rows

Le quitamos la geometría al dataset generado:

Total_Bibliotecas_OSM_en_Barrios <- Total_Bibliotecas_OSM_en_Barrios %>% 
  st_set_geometry(NULL)

head(Total_Bibliotecas_OSM_en_Barrios)
## # A tibble: 6 x 2
##   Neighborho total
##   <fct>      <int>
## 1 10             2
## 2 11            18
## 3 12            20
## 4 13             2
## 5 14             3
## 6 16            28

Utilizando el comando left_join() unificamos el datasat Total_Bibliotecas_OSM_en_Barrios con Barrios_Boston_OK:

Total_Bibliotecas_OSM_en_Barrios <- Barrios_Boston_OK %>% 
  left_join(Total_Bibliotecas_OSM_en_Barrios, by="Neighborho")

head(Total_Bibliotecas_OSM_en_Barrios)
## Simple feature collection with 6 features and 8 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -71.14779 ymin: 42.26761 xmax: -71.05588 ymax: 42.35237
## CRS:            EPSG:4326
##   OBJECTID             Name      Acres Neighborho SqMiles  ShapeSTAre
## 1       27       Roslindale 1605.56824         15    2.51  69938272.9
## 2       28    Jamaica Plain 2519.24539         11    3.94 109737890.8
## 3       29     Mission Hill  350.85356         13    0.55  15283120.0
## 4       30         Longwood  188.61195         28    0.29   8215903.5
## 5       31      Bay Village   26.53984         33    0.04   1156070.8
## 6       32 Leather District   15.63991         27    0.02    681271.7
##   ShapeSTLen total                       geometry
## 1  53563.913    NA MULTIPOLYGON (((-71.12593 4...
## 2  56349.937    18 MULTIPOLYGON (((-71.10499 4...
## 3  17918.724     2 MULTIPOLYGON (((-71.09043 4...
## 4  11908.757    72 MULTIPOLYGON (((-71.09811 4...
## 5   4650.635     2 MULTIPOLYGON (((-71.06663 4...
## 6   3237.141    NA MULTIPOLYGON (((-71.05838 4...

Ahora mapeamos:

ggplot() +
  geom_sf(data = Total_Bibliotecas_OSM_en_Barrios, aes(fill = total)) + 
  geom_sf_text(data = Total_Bibliotecas_OSM_en_Barrios, aes(label = Name), size = 2, colour = "black") +
  labs(title = "Total Bibliotecas por Barrio",
         subtitle = "Boston",
         fill = "Total Bibliotecas",
         caption= "Fuente: https://data.boston.gov/dataset/boston-neighborhoods y
                   OpenStreetMap",
         y="",
         x="") +
  scale_fill_gradient(low="yellow", high="red")
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data
## Warning: Computation failed in `stat_sf_coordinates()`:
## Evaluation error: TopologyException: Input geom 1 is invalid: Ring Self-intersection at or near point -71.125747341600402 42.27233853935423 at -71.125747341600402 42.27233853935423.

Para obtener una mejor observación de la distribución de las Bibliotecas por Barrio realizaremos un gráfico de puntos, de este modo tendremos la información completa para comprar nuestra información actual con la Tarea 1:

ggplot(Total_Bibliotecas_OSM_en_Barrios) +
  geom_point(aes(x=Neighborho, y=total), color="red") + 
  labs(title="Bibliotecas por Barrio",
       subtitle="Boston",
       caption="OpenStreetMap",
       x="Barrio ID",
       y="Total")
## Warning: Removed 3 rows containing missing values (geom_point).

Se advierte que las Biliotecas de OSM alcanzan 23 Barrios de la ciudad, mientras que Bibliotecas Públicas estan presentes en 18 Barrios. Asimismo, la cantidad de entidades que pueden observarse en cada uno de ellos es muy distinto. En el dataset de Biliotecas de OSM los barrios tienen un mínimo de 10 establecimientos, mientras que el dataset de Bibliotecas Públicas tiene 1 establecimiento. Además, en el primer dataset mencionado, se distingue un rango entre 10 y 75 Bibliotecas, primando el Barrio ID 6 con mas de 200 entidades. Por otro lado, en los barrios donde se encuentran presentes las Bibliotecas Públicas, éstas oscilan entre 1 y 2 entidades, pero se asemeja con la concentración de bibliotecas en el Barrio ID 6. Podemos deducir que esta es el área central de la ciudad, donde se concentran actividades, servicios e infraestructura. Los Barrios ID 11, 12, 16 y 25 contaban con 2 Bibliotecas Públicas y considerando las Bibliotecas de OSM llegan a un alcance aproximado de 25 establecimientos. Al considerar las Bibliotecas OSM se observan nuevos barrios que no contaban con instituciones, alcanzando altos números como los Barrios ID 2, 24, 28 y 34 que tienen un promedio aproximado entre 50 y 75 bibliotecas. Podría considerarse a dichos distritos como núcleos urbanos con alta densidad, necesarios de contar con el mencionado servicio.

Por otro lado realizaremos un mapa de coropletas, para contemplar las densidad de Bibliotecas por Barrio:

ggplot() +
  geom_sf(data = Total_Bibliotecas_OSM_en_Barrios, aes(fill = total/Acres)) + 
  geom_sf_text(data = Total_Bibliotecas_OSM_en_Barrios, aes(label = Name), size = 2, colour = "black") +
  labs(title = "Total Bibliotecas por Barrio",
         subtitle = "Boston",
         fill = "Total Bibliotecas",
         caption= "Fuente: https://data.boston.gov/dataset/boston-neighborhoods y
                   OpenStreetMap",
         y="",
         x="") +
  scale_fill_gradient(low="yellow", high="red")
## Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
## give correct results for longitude/latitude data
## Warning: Computation failed in `stat_sf_coordinates()`:
## Evaluation error: TopologyException: Input geom 1 is invalid: Ring Self-intersection at or near point -71.125747341600402 42.27233853935423 at -71.125747341600402 42.27233853935423.

Verificamos al igual que en la Tarea 1, que la distribución de Bibliotecas no es coincidente con la superficie de los Barrios. Los Barrios con menor superfice, suelen tener la mayor cantidad de densidad de Bibliotecas. A diferencia que en el dataset de Bibliotecas Públicas, donde la mayor densidad se disntiguía en el Barrio ID 14, aquí la mayor densidad se da en el Barrio ID 13. La similitud entre ambos dataset Bibliotecas Públicas y Bibliotecas OMS se da en que los barrios más pequeños concentran la mayor densidad de bibliotecas.