Subte población y q pasajeros

Descargamos la base de datos de radios censales de 2010

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, 600)

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")
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: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -58.49298 ymin: -34.64879 xmax: -58.36331 ymax: -34.55017
## Geodetic CRS:  WGS 84
## # A tibble: 6 × 3
##   LINEA Hogares_NBI                                                     geometry
##   <chr>       <dbl>                                                <POLYGON [°]>
## 1 A           27903 ((-58.40396 -34.60491, -58.40396 -34.60484, -58.40405 -34.6…
## 2 B           16042 ((-58.41689 -34.60746, -58.41689 -34.60754, -58.41682 -34.6…
## 3 C           14868 ((-58.38675 -34.6111, -58.38675 -34.61114, -58.3868 -34.611…
## 4 D           10216 ((-58.41974 -34.57612, -58.41974 -34.57604, -58.41979 -34.5…
## 5 E           22259 ((-58.46264 -34.64867, -58.46264 -34.64871, -58.46255 -34.6…
## 6 H           15953 ((-58.3945 -34.6303, -58.39442 -34.63027, -58.39442 -34.630…

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.43872 -34.62203, -58.43869 -34.62202, -58.43869 -34.…
## 2 f            7622 (((-58.39532 -34.60146, -58.39532 -34.60146, -58.39531 -34.…
## 3 g            2974 (((-58.44467 -34.60999, -58.44445 -34.60993, -58.44436 -34.…

# Para las líneas clásicas (estaciones_subte)
# Población 2010
pop_2010_classical <- st_join(estaciones_buffer, radios, join = st_intersects) %>%
  group_by(LINEA) %>%
  summarise(pop_2010 = sum(as.numeric(TOTAL_POB), na.rm = TRUE))

# Transformar ambas capas a CRS 4326
estaciones_buffer <- st_transform(estaciones_buffer, crs = 4326)
poblacion_2022 <- st_transform(poblacion_2022, crs = 4326)

# Población 2022 para líneas clásicas
pop_2022_classical <- st_join(estaciones_buffer, poblacion_2022, join = st_intersects) %>%
  group_by(LINEA) %>%
  summarise(pop_2022 = sum(as.numeric(POB_TOT), na.rm = TRUE))

# Para las líneas nuevas (lineas)
# Población 2010
pop_2010_new <- st_join(lineas_buffer, radios, join = st_intersects) %>%
  group_by(linea) %>%
  summarise(pop_2010 = sum(as.numeric(TOTAL_POB), na.rm = TRUE))

# Población 2022
pop_2022_new <- st_join(lineas_buffer, poblacion_2022, join = st_intersects) %>%
  group_by(linea) %>%
  summarise(pop_2022 = sum(as.numeric(POB_TOT), na.rm = TRUE))

# Combinar ambas series para cada conjunto
# Líneas clásicas: se renombra para unificar el nombre de la columna que identifica la línea

pop_classical <- st_drop_geometry(pop_2010_classical) %>%
  rename(linea = LINEA) %>%
  full_join(pop_2022_classical %>% rename(linea = LINEA), by = "linea")

# Líneas nuevas ya tienen el campo 'linea'
pop_new <- full_join(pop_2010_new, st_drop_geometry(pop_2022_new), by = "linea")

# Unir ambas tablas en una sola
pop_comparison <- bind_rows(pop_classical, pop_new)

pop_comparison <- pop_comparison %>%
  mutate(
    pct_change = paste0(round((pop_2022 - pop_2010) / pop_2010 * 100, 2), "%"),
    annual_growth_rate = paste0(round(((pop_2022 / pop_2010)^(1/(2022 - 2010)) - 1) * 100, 2), "%") 
  )


pop_comparison_no_geometry <- subset(pop_comparison, select = -geometry)


colnames(pop_comparison_no_geometry) <- c("Linea", "Poblacion_2010", "Poblacion_2022", "Cambio_Porcentual", "Tasa_Crecimiento_Anual")
# Visualizar la tabla con las nuevas columnas


datatable(pop_comparison_no_geometry, options = list(pageLength = 10))

# Visualizar la tabla

# Calculate area in km² and population density for radios (2010 data)
radios <- radios %>% 
  mutate(area_km2 = as.numeric(st_area(geometry)) /1e6,
         pop_density = as.numeric(TOTAL_POB) / area_km2)

library(ggplot2)
library(plotly)

## Warning: package 'plotly' was built under R version 4.3.2

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggmap':
## 
##     wind

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

library(sf)
library(viridis)

## Loading required package: viridisLite

radios <- radios %>%
  mutate(across(where(is.character), ~ iconv(., from = "", to = "UTF-8")))


# Create a ggplot map for the radios layer
p_radios <- ggplot(radios) +
  geom_sf(aes(fill = pop_density), color = NA) +
  scale_fill_viridis_c(option = "plasma") +
  theme_minimal() +
  labs(
    title = "Population Density Map (Radios 2010)",
    fill = "People/km2"
  )

# Convert the ggplot to an interactive Plotly map
p_radios_plotly <- ggplotly(p_radios)

# Display the interactive map
p_radios_plotly

mean(radios$pop_density)

## [1] 27470.43

poblacion_2022 = poblacion_2022 %>% mutate(DEN_POB = DEN_POB *100)

poblacion_2022 <- poblacion_2022 %>%
  mutate(across(where(is.character), ~ iconv(., from = "", to = "UTF-8")))

# Create a ggplot map for the radios layer
p_radios_2022 <- ggplot(poblacion_2022) +
  geom_sf(aes(fill = DEN_POB), color = NA) +
  scale_fill_viridis_c(option = "plasma") +
  theme_minimal() +
  labs(
    title = "Population Density Map (Radios 2022)",
    fill = "People/km2"
  )

# Convert the ggplot to an interactive Plotly map
p_radios_plotly_2022 <- ggplotly(p_radios_2022)

# Display the interactive map

p_radios_plotly_2022