library("tidyverse")
## Warning: package 'tidyverse' was built under R version 4.2.3
## Warning: package 'ggplot2' was built under R version 4.2.3
## Warning: package 'tibble' was built under R version 4.2.3
## Warning: package 'tidyr' was built under R version 4.2.3
## Warning: package 'readr' was built under R version 4.2.3
## Warning: package 'purrr' was built under R version 4.2.3
## Warning: package 'dplyr' was built under R version 4.2.3
## Warning: package 'stringr' was built under R version 4.2.3
## Warning: package 'forcats' was built under R version 4.2.3
## Warning: package 'lubridate' was built under R version 4.2.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.1     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.2     ✔ tibble    3.2.1
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library("ggplot2")
library("sf")
## Linking to GEOS 3.9.1, GDAL 3.3.2, PROJ 7.2.1; sf_use_s2() is TRUE
library("osmdata")
## Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright
library("rmarkdown")
library("lubridate")
library("ggmap") 
## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.
## Please cite ggmap if you use it! See citation("ggmap") for details.
library("mapview")
library("hrbrthemes")
## NOTE: Either Arial Narrow or Roboto Condensed fonts are required to use these themes.
##       Please use hrbrthemes::import_roboto_condensed() to install Roboto Condensed and
##       if Arial Narrow is not on your system, please see https://bit.ly/arialnarrow
library("RColorBrewer")
library("readr")
library("leaflet")
library("XML")
## Warning: package 'XML' was built under R version 4.2.1
library("osmextract")
## Warning: package 'osmextract' was built under R version 4.2.1
## Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright.
## Check the package website, https://docs.ropensci.org/osmextract/, for more details.
library("forcats")
library("spatialreg")
## Warning: package 'spatialreg' was built under R version 4.2.2
## Loading required package: spData
## Warning: package 'spData' was built under R version 4.2.2
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
## Loading required package: Matrix
## Warning: package 'Matrix' was built under R version 4.2.1
## 
## Attaching package: 'Matrix'
## 
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
library("spdep")
## Warning: package 'spdep' was built under R version 4.2.2
## Loading required package: sp
## 
## Attaching package: 'spdep'
## 
## The following objects are masked from 'package:spatialreg':
## 
##     get.ClusterOption, get.coresOption, get.mcOption,
##     get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
##     set.coresOption, set.mcOption, set.VerboseOption,
##     set.ZeroPolicyOption
library("rgdal")
## Warning: package 'rgdal' was built under R version 4.2.1
## Please note that rgdal will be retired by the end of 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## 
## rgdal: version: 1.5-32, (SVN revision 1176)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.4.3, released 2022/04/22
## Path to GDAL shared files: C:/Users/sixto/AppData/Local/R/win-library/4.2/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Users/sixto/AppData/Local/R/win-library/4.2/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.5-0
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
library("readxl")
## Warning: package 'readxl' was built under R version 4.2.3
library("maptools")
## Warning: package 'maptools' was built under R version 4.2.1
## Checking rgeos availability: TRUE
## Please note that 'maptools' will be retired by the end of 2023,
## plan transition at your earliest convenience;
## some functionality will be moved to 'sp'.

Descargamos la base de datos de radios censales de 2010

radios = st_read("https://cdn.buenosaires.gob.ar/datosabiertos/datasets/direccion-general-de-estadisticas-y-censos/informacion-censal-por-radio/caba_radios_censales.geojson")
## Reading layer `informacion_censal_por_radio_wgs84' from data source 
##   `https://cdn.buenosaires.gob.ar/datosabiertos/datasets/direccion-general-de-estadisticas-y-censos/informacion-censal-por-radio/caba_radios_censales.geojson' 
##   using driver `GeoJSON'
## Simple feature collection with 3554 features and 15 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -58.53152 ymin: -34.70529 xmax: -58.33514 ymax: -34.52755
## Geodetic CRS:  WGS 84

Vamos a tener que ajustar varias cosas.

Descargamos la base de datos del subte

estaciones_subte = st_read("https://cdn.buenosaires.gob.ar/datosabiertos/datasets/sbase/subte-estaciones/estaciones-de-subte.geojson")
## Reading layer `estaciones_de_subte_WGS84' from data source 
##   `https://cdn.buenosaires.gob.ar/datosabiertos/datasets/sbase/subte-estaciones/estaciones-de-subte.geojson' 
##   using driver `GeoJSON'
## Simple feature collection with 90 features and 3 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -58.48639 ymin: -34.64331 xmax: -58.36993 ymax: -34.55564
## Geodetic CRS:  WGS 84

Vamos a calcular las NBI de cada radio censal

 radios = radios %>% mutate(nbi_per = as.numeric(radios$H_CON_NBI) / as.numeric(radios$T_HOGAR))

Subte Buffer

estaciones_buffer = st_buffer(estaciones_subte, 300)
mapview(estaciones_buffer)

Descargamos la cantidad de pasajeros del subte.

result <- st_join(estaciones_buffer, radios, join = st_intersects)
result <- result %>%
  group_by(ID.x) %>%
  summarise(total_population = sum(as.numeric(TOTAL_POB)))
mapview(result,  zcol ="total_population" )

Nuevas líneas de subte

linea_f =  st_read("C:/Users/sixto/OneDrive/Documentos/LSE/Urban Economics/subtenuevos/Linea_F.geojson")
## Reading layer `Linea F' from data source 
##   `C:\Users\sixto\OneDrive\Documentos\LSE\Urban Economics\subtenuevos\Linea_F.geojson' 
##   using driver `GeoJSON'
## Simple feature collection with 15 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -58.42135 ymin: -34.64661 xmax: -58.373 ymax: -34.58097
## Geodetic CRS:  WGS 84
linea_f$linea <- "f"
mapview(linea_f)
linea_g =  st_read("C:/Users/sixto/OneDrive/Documentos/LSE/Urban Economics/subtenuevos/Linea_g.geojson")
## Reading layer `Linea G' from data source 
##   `C:\Users\sixto\OneDrive\Documentos\LSE\Urban Economics\subtenuevos\Linea_G.geojson' 
##   using driver `GeoJSON'
## Simple feature collection with 15 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -58.49278 ymin: -34.6075 xmax: -58.37499 ymax: -34.59149
## Geodetic CRS:  WGS 84
linea_g$linea <- "g"
mapview(linea_g)
linea_I =  st_read("C:/Users/sixto/OneDrive/Documentos/LSE/Urban Economics/subtenuevos/Linea_I.geojson")
## Reading layer `Linea I' from data source 
##   `C:\Users\sixto\OneDrive\Documentos\LSE\Urban Economics\subtenuevos\Linea_I.geojson' 
##   using driver `GeoJSON'
## Simple feature collection with 10 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -58.44577 ymin: -34.63104 xmax: -58.42157 ymax: -34.58133
## Geodetic CRS:  WGS 84
linea_I$linea <- "I"
mapview(linea_I)

` Cambiamos el nombre de estaciones que tienen nombres idénticos

linea_I <- linea_I %>% 
  mutate(Name = ifelse(Name == "Rivadavia", "Rivadavia_I", Name))
linea_I <- linea_I %>% 
  mutate(Name = ifelse(Name == "Cid Campeador", "Cid Campeador_I", Name))
linea_f <- linea_f %>% 
  mutate(Name = ifelse(Name == "Corrientes", "Corrientes_f", Name))
linea_I <- linea_I %>% 
  mutate(Name = ifelse(Name == "Plaza Italia", "Plaza Italia_I", Name))
lineas = rbind(linea_f,linea_g, linea_I)
mapview(lineas, zcol = "linea")
lineas_buffer = lineas %>% st_buffer(300)
mapview(lineas_buffer, zcol = "linea")
interseccion_nuevo_subte_buffer <- st_join(lineas_buffer, radios, join = st_intersects)
interseccion_nuevo_subte_buffer <- interseccion_nuevo_subte_buffer %>%
  group_by(Name, linea) %>%
  summarise(total_population = sum(as.numeric(TOTAL_POB)))
## `summarise()` has grouped output by 'Name'. You can override using the
## `.groups` argument.
Lineas_G_F_I = interseccion_nuevo_subte_buffer

mapview(Lineas_G_F_I, zcol = "total_population")
names(Lineas_G_F_I)
## [1] "Name"             "linea"            "total_population" "geometry"
summary_by_linea <- Lineas_G_F_I %>%
  group_by(linea) %>%
  summarise(total_population = sum(total_population))
library("DT")
## Warning: package 'DT' was built under R version 4.2.3
datatable(summary_by_linea, options = list(pageLength = 10))

Ahora quiero comparar con NBI entre todas las líneas.

Líneas clásicas

result <- st_join(estaciones_buffer, radios, join = st_intersects)
result_NBI <- result %>%
  group_by(LINEA) %>%
  summarise(Hogares_NBI = sum(as.numeric(H_CON_NBI)))
print(result_NBI)
## Simple feature collection with 6 features and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -58.48968 ymin: -34.64604 xmax: -58.36661 ymax: -34.5529
## Geodetic CRS:  WGS 84
## # A tibble: 6 × 3
##   LINEA Hogares_NBI                                                     geometry
##   <chr>       <dbl>                                           <MULTIPOLYGON [°]>
## 1 A            8595 (((-58.41444 -34.61342, -58.41439 -34.6134, -58.41436 -34.6…
## 2 B            5284 (((-58.4131 -34.60655, -58.4131 -34.60657, -58.41306 -34.60…
## 3 C            4852 (((-58.38314 -34.61103, -58.38314 -34.61107, -58.38319 -34.…
## 4 D            3483 (((-58.38386 -34.60408, -58.38386 -34.60411, -58.38386 -34.…
## 5 E            8178 (((-58.38087 -34.61275, -58.38087 -34.61283, -58.38087 -34.…
## 6 H            4478 (((-58.39647 -34.63396, -58.39647 -34.6339, -58.39652 -34.6…

Líneas nuevas

interseccion_nuevo_subte_buffer <- st_join(lineas_buffer, radios, join = st_intersects)
result_NBI_nuevas_lineas <- interseccion_nuevo_subte_buffer %>%
  group_by(linea) %>%
  summarise(Hogares_NBI = sum(as.numeric(H_CON_NBI)))
print(result_NBI_nuevas_lineas)
## Simple feature collection with 3 features and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -58.49605 ymin: -34.64934 xmax: -58.36968 ymax: -34.57823
## Geodetic CRS:  WGS 84
## # A tibble: 3 × 3
##   linea Hogares_NBI                                                     geometry
##   <chr>       <dbl>                                           <MULTIPOLYGON [°]>
## 1 I            1721 (((-58.44015 -34.62754, -58.44012 -34.62753, -58.44012 -34.…
## 2 f            7622 (((-58.37575 -34.64813, -58.37575 -34.64817, -58.3757 -34.6…
## 3 g            2974 (((-58.37811 -34.59234, -58.37811 -34.59242, -58.37811 -34.…