Cartografiando la Equidad Urbana

Especialización: Big Data e Inteligencia Territorial

1.Datos abiertos utilizados:

Límites de barrios de Valecia:

Datos de población por barrio

Zonas verdes urbanas (parques, jardines, plazas)

2.Planteamiento del problema

La calidad de vida urbana está estrechamente relacionada con el acceso a espacios verdes que benefician la salud física, mental y ambiental de la población y permiten mitigando la alta densidad poblacional que tienen las ciudades. Según la OMS, se recomienda un mínimo de 9 m² de área verde por habitante.

Este trabajo analizó la densidad poblacional de Valencia, España en relación con las áreas verdes como medio de mitigación de problemas socioambiental a nivel de barrios, con datos poblacionales hasta 2024. Se calculó el Índice Verde Urbano (IVU) y la Densidad Poblacional para evaluar si es necesario crear nuevos espacios que mejoren la calidad de vida local.

3.Formulación del problema

¿Los barrios más densamente poblados de Valencia tienen zonas verdes proporcionales para aumentar la calidad de vida local?

4.Metodología utilizada

5. Resultados de la investigación

Inciamos cargando librerías

#Librerías importantes para resolver la investigación.
library(tidyverse) #manipulación de dataframes.
library(sf) #manipulción de datos geográficos.
library(leaflet) # mapas dinámicos.
library(tidygeocoder) #ruteo de distancia.
library(osmdata) #datos libres de OSM.
library(ggspatial)
library(ggmap)

Cargamos datasets geográfico y dataframe de excel.

#Usamos librería sf para datos geográficos.
barrios <- st_read("Geografico/barrios.geojson")
## Reading layer `barrios' from data source 
##   `C:\R_directory\2. Datascience_II - Geografico\Trabajo final\Geografico\barrios.geojson' 
##   using driver `GeoJSON'
## Simple feature collection with 88 features and 6 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -0.432535 ymin: 39.27893 xmax: -0.2753685 ymax: 39.56659
## Geodetic CRS:  WGS 84
distritos <- st_read("Geografico/distritos.geojson")
## Reading layer `distritos' from data source 
##   `C:\R_directory\2. Datascience_II - Geografico\Trabajo final\Geografico\distritos.geojson' 
##   using driver `GeoJSON'
## Simple feature collection with 22 features and 6 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -0.432535 ymin: 39.27893 xmax: -0.2753685 ymax: 39.56659
## Geodetic CRS:  WGS 84
#Librería Fread para dataframes
poblacion <- read.csv("dataframes/Poblacion.csv",
                   encoding = "UTF-8",
                   sep = ";")

Revisamos el tipo de objeto cargado.

class(barrios)
## [1] "sf"         "data.frame"
class(distritos)
## [1] "sf"         "data.frame"
class(poblacion)
## [1] "data.frame"

Vemos que barrios y distritos es de tipo sf, mientras que, población es de tipo data.frame. Aplicamos un resumen para conocer su contenido.

summary(barrios)
##   coddistbar           nombre           coddistrit         codbarrio        
##  Length:88          Length:88          Length:88          Length:88         
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##                                                                             
##  gis_gis_barrios_area geo_point_2d                geometry 
##  Min.   :  279952     Length:88          POLYGON      :88  
##  1st Qu.: 2470073     Class :character   epsg:4326    : 0  
##  Median : 3854929     Mode  :character   +proj=long...: 0  
##  Mean   : 5027684                                          
##  3rd Qu.: 6020852                                          
##  Max.   :25073008                                          
##  NA's   :32

Vemos que el Municipio de Valencia tiene 88 barrios. El campo gis_gis_barrios_area es donde se almacenada el área de cada barrio, el cual contiene 32 NAs, los cuales debemos corregir para calcular la densidad población. La geometría es de tipo polígono y proyectado 4326. No tiene las variables X-Y con las coordenadas para graficar el texto de cada atributo. Vemos también que tiene la variable coddistbar de tipo character, muy importante para vincular con población y poder obtener la población por cada barrio.

Aplicamos un resumen a población

summary(poblacion)
##    Barrios            coddistbar        Total          Varones     
##  Length:88          Min.   : 11.0   Min.   :    0   Min.   :    0  
##  Class :character   1st Qu.: 59.5   1st Qu.: 3936   1st Qu.: 1953  
##  Mode  :character   Median :103.5   Median : 7142   Median : 3364  
##                     Mean   :104.3   Mean   : 9439   Mean   : 4481  
##                     3rd Qu.:151.2   3rd Qu.:12239   3rd Qu.: 5689  
##                     Max.   :200.0   Max.   :43916   Max.   :21237  
##     Mujeres     
##  Min.   :    0  
##  1st Qu.: 1982  
##  Median : 3770  
##  Mean   : 4958  
##  3rd Qu.: 6420  
##  Max.   :22679

Para este caso, vemos que tenemos 88 barrios con su población. La población mínima es de 0 habitantes, media de 9439 de los cuales 4481 son Varones y 4958 son Mujeres. Por otro lado, la población máxima es de 43916 personas de las cuales 21237 son Varores y 22679 son Mujeres. Tenemos el coddistbar para vincularnos con barrios pero su estructura es de tipo entero porque nos muestra estadísticos. Esto debemos cambiarlo a character antes de poder utilizarlo con barrios para que la función lef_join no arroje inconvenientes.

Graficamos los barrios con un mapa base de stadia_map. Primero calculamos el bounding box para descargar el mapa.

#calculamos el bounding box.
bbox <- as.numeric(st_bbox(barrios)) #tipo numérico.

#creamos el mapa que vamos a usar.
mapa <- get_stadiamap(bbox = bbox,
                      maptype = "stamen_terrain_background",
                      zoom = 10)
## ℹ © Stadia Maps © Stamen Design © OpenMapTiles © OpenStreetMap contributors.

Para mostrar las etiquetas de cada barrio necesitamos calcular los centroides de cada barrio, las coordenadas X-Y para usarlo posteriormente como etiqueta.

#calculamos las coordenadas x-y para graficar etiquetas.
barrios_coords <- barrios |>
  st_centroid() |> 
  st_coordinates() |> 
  cbind(barrios)

#graficamos mapa de barrios
ggmap(mapa) +
  geom_sf(data = barrios, aes(fill = nombre), color = "black", alpha = 0.7, inherit.aes = FALSE) +
  geom_text(data = barrios_coords, aes(x = X, y = Y, label = coddistbar),
            size = 1, color = "black", fontface = "bold") +
  labs(
    title = "Distribución de barrios",
    subtitle = "Valencia, España",
    caption = "Fuente: Geoportal del Ayuntamiento de Valencia") +
  theme_void() +
  theme(legend.position = "none")

Para calcular la densidad población debemos tener el área de los barrios en km2. Por lo cual, primero recalculamos el valor del área en km2 y revisamos su distribución.

#calculamos area nuevamente
barrios <- barrios |> 
  mutate(gis_gis_barrios_area =  as.numeric(st_area(geometry)/ 1e6))

#revisamos el área de los barrios.
summary(barrios$gis_gis_barrios_area)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##  0.09355  0.32813  0.51113  1.55493  1.23631 36.91034

Podemos observar que tenemos barrios que tiene un área menor a 1 km2 (0,09355), mientras otros barrios tienen una extensión de 36,91 km2 con una media de 1,55 km2. Este primer análisis exploratorio de datos nos está indicando que la gran mayoría son barrios pequeños con poca superficie.

Se procede con la unión de población y barrios a partir del campo en común coddistbar. Para ello, debemos cambiar el tipo de dato de la variable poblacion$coddistbar para que sea de tipo character.

#cambiamos tipo de dato a poblacion
poblacion <- poblacion |> 
  mutate(coddistbar = as.character(coddistbar))

#join con poblacion
barrios <- barrios |> 
  left_join(poblacion, by = "coddistbar")

Ahora calculamos la densidad población con el área y la población y revisamos su resultado

barrios <-  barrios |> 
  mutate(densidad_hab_km2 = Total / gis_gis_barrios_area)

summary(barrios$densidad_hab_km2)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       0    6849   17780   18709   31050   55530

Dentro del barrio RAFALELL-VISTABELLA no tenemos regitrados datos por ello nos da un mínimo de 0 hab/km2. La media es de 18.709 hab/km2 y un valor máximo de 55.530 hab/km2. Graficamos el resultado para conocer su distribución en el territorio.

#graficamos el resultado numérico 
ggplot () +
  geom_sf(data = barrios, aes(fill = densidad_hab_km2), inherit.aes = FALSE) +
  scale_fill_distiller(palette="Oranges", direction = 1) +
  labs(
    title = "Densidad poblacional a nivel de barrios",
    subtitle = "Valencia, España",
    fill = "hab/km2",
    caption = "Fuente: Geoportal del Ayuntamiento de Valencia") +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, size = 13))

Generamos una gráfica faceteado de barra para conocer los 10 barrios con mayor y menor densidad poblacional.

df_topB <- barrios %>% 
  arrange(desc(densidad_hab_km2)) %>% 
  slice(1:10) %>% 
  mutate(grupo = "Top 10 más densos")

df_bottomB <- barrios %>%
  arrange(densidad_hab_km2) %>%
  slice(1:10) %>%
  mutate(grupo = "Top 10 menos densos")

# 2) Unir y transformar en un solo df
df_facet_barrios <- bind_rows(df_topB, df_bottomB)

# 3) Graficar faceteado
ggplot(df_facet_barrios, aes(x = reorder(nombre, densidad_hab_km2),
                     y = densidad_hab_km2,
                     fill = densidad_hab_km2)) +
  geom_col() +
  coord_flip() +
  scale_fill_viridis_c(option = "plasma") +
  facet_wrap(~ grupo, ncol = 1, scales = "free_y") +
  labs(
    title = "Densidad poblacional (Top 10 más y menos densos)",
    x = "Barrio",
    y = "Densidad (hab/km²)",
    fill = "hab/km²"
  ) +
  theme_minimal() +
  theme(
    strip.text = element_text(face = "bold", size = 12)
  )

Podemos ver que el barrio más densamente poblado es El Calvari con 55.530 hab/km2, siguiéndole Els Orriols con 45575 hab/km2. En el otro extremo, tenemos barrios con densidades muy bajas. Por ejemplo, El Palmar tiene una densidad de 20 hab/km2 mientras que Borboto tiene 521 hab/km2. Rafalell-Vistabella no tiene registros de población por lo que no se pudo calcular la densidad poblacional.

Se realizó una categorización para mostrar el resultado por niveles: Muy bajo, Bajo, Medio, Alto y Muy Alto y revisamos el conteo de cada categoría para ver cuantos barrios tenemos dentro de cada categoría.

#Categorizamo por niveles
barrios <-  barrios |> 
  mutate( densidad_clase = case_when(densidad_hab_km2 >= 0 & densidad_hab_km2 <10000 ~ " Muy Baja",
                                     densidad_hab_km2 >= 10000 & densidad_hab_km2 < 20000 ~"Baja",
                                     densidad_hab_km2 >= 20000 & densidad_hab_km2 < 30000 ~ "Media",
                                     densidad_hab_km2 >= 30000 & densidad_hab_km2 < 40000 ~ "Alta",
                                     densidad_hab_km2 >= 40000 ~ "Muy Alto"))

# Agrupamos para un conteo de las clases
barrios |> 
  group_by(densidad_clase) |> 
  summarise(recuento = n())
## Simple feature collection with 5 features and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -0.432535 ymin: 39.27893 xmax: -0.2753685 ymax: 39.56659
## Geodetic CRS:  WGS 84
## # A tibble: 5 × 3
##   densidad_clase recuento                                               geometry
##   <chr>             <int>                                     <MULTIPOLYGON [°]>
## 1 " Muy Baja"          28 (((-0.412439 39.49773, -0.4124427 39.49774, -0.412444…
## 2 "Alta"               18 (((-0.3397535 39.47156, -0.3394047 39.47227, -0.34039…
## 3 "Baja"               22 (((-0.3538934 39.47585, -0.3527422 39.47799, -0.35467…
## 4 "Media"              15 (((-0.3502911 39.47823, -0.3512561 39.4782, -0.351506…
## 5 "Muy Alto"            5 (((-0.3905235 39.48823, -0.3920742 39.48755, -0.39327…

Podemos observar que tenemos 28 barrios que tienen un densidad poblacional muy baja, seguido de 22 barrios con densidad baja, 15 con densidad media, 18 con alta y 5 con una densidad muy alta.

Graficamos el resultado por categorías.

#graficamos el resultado por categoría
ggplot () +
  geom_sf(data = barrios, aes(fill = densidad_clase), inherit.aes = FALSE) +
  scale_fill_brewer(palette = "Pastel1", direction = -1) +
  labs(
    title = "Densidad poblacional a nivel de barrios",
    subtitle = "Valencia, España",
    fill = "Categoría",
    caption = "Fuente: Geoportal del Ayuntamiento de Valencia") +
  theme_void() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, size = 13))

Una vez categorizados, podemos identificar un patrón bastante claro es que la densidad poblacional en la periferia de la ciudad es muy baja porque existe mayor extensión de terreno. Además, podemos ver que en el centro de la ciudad la densidad es baja. En los alrededores vemos que ya tenemos densidades medias, alta y muy altas.

Para calcular el índice de verde urbano, necesitamos obtener las áreas verdes por barrio. Como no tenemos una capa de información que detalle las áreas verdes, vamos a obtenerla a partir de Open Street Map y luego hacer un spatial join para obtener el barrio que le corresponde.

#Paso 1 obtenemos el bbox de los elementos
areas_verdes <- opq(bbox)
#Paso 2 obtenemos los features de los servicios financieros
areas_verdes <- areas_verdes |>  
  add_osm_feature(key = "leisure", value = c("park", "garden"))
#Paso 3 descargamos los feature como sf
areas_verdes <-  osmdata_sf(areas_verdes)
#Paso 4 dejamos solo los puntos.
areas_verdes <- areas_verdes$osm_polygons
#Paso 5 dejamos solo las áreas que quedan dentro de los barrios.
areas_verdes <- st_intersection(areas_verdes, barrios)

Con la descarga de OSM tenemos 1803 elementos que cumplen con la condición parques y jardines dentro de barrios. Vamos a graficar el resultado para concoer su ubicación.

ggmap(mapa) +
  geom_sf(data = barrios, fill = NA, color = "black", lwd = 0.75, inherit.aes = FALSE) +
  geom_sf(data = areas_verdes, aes(color = leisure), inherit.aes = FALSE) +
  labs(title = "Mapa de áreas verdes a nivel de barrios",
       subtitle = "Valencia, España",
       caption = "Fuente: OpenStreetMap",
       color = "leisure") +
  scale_color_manual(values = c("firebrick3", "darkorange1", "blue4", "grey")) +
  theme_void() +
  theme(title = element_text(size = 10, face = "bold"), #tamaño de titulo del mapa
        plot.caption = element_text(face = "italic", colour = "gray35",size = 7)) #tamaño de nota al pie

El resultado arrojo también valores que no deberían estar dentro de la selección como playground, swimming_pool y NA. Vamos a eliminar estos valores para continuar.

areas_verdes <-  areas_verdes |> 
  filter(!leisure %in% c("playground", "swimming_pool") & !is.na(leisure))

Ahora si nos quedamos con 1624 elemetos que cumplen con la condición de parques y jardines.

Para continuar con el cálculo del índice verde urbano como m2/hab, necesitamos calcular la sumatoria de las áreas de parques y jardines para luego calcular el índice de verde urbano.

Para ello, debemos realizar ciertas transformaciones para pasar datos de un objeto a otro, realizar agregaciones y demás procesamiento espacial.

#Calculamos las áreas verdes.
areas_verdes <-  areas_verdes |> 
  select(osm_id, name, leisure) |>
  mutate(area_m2 = as.numeric(st_area(geometry)))

#Spatial Join para conservar campos específicos  
areas_verdes_barrios <- st_join(areas_verdes, barrios |> select(coddistbar, Total)) # solo atributos necesarios.

#Agrupamos por barrios y calculamos área verdes totales y obtenemos el IVU
barrios_IVU <- areas_verdes_barrios |> 
  group_by (coddistbar, Total) |> 
  summarise (area_m2_Total = sum(area_m2, na.rm = TRUE)) |> 
  ungroup() |> 
  mutate(IVU = area_m2_Total / Total)
## `summarise()` has grouped output by 'coddistbar'. You can override using the
## `.groups` argument.

Hacemos un join sencillo para unir la información calculada del IVU al objeto de barrios y representar luego su resultado.

#join de atributos
barrios_IVU_f <- barrios |> 
  left_join(
    barrios_IVU |> st_drop_geometry(), 
    by = "coddistbar")

Procedemos a graficar el resultado para conocer el Índice de Verde Urbano a nivel de barrio con la población y las áreas verdes totales por barrio.

ggplot() +
  geom_sf(data = barrios_IVU_f, aes(fill = IVU), color = "black") +
  labs (title = "Maapa con Índice de Verde Urbano a nivel de Barrios",
        subtitle = "Valencia, España",
        caption = "Ayundamiento de Valencia, Anuario Estadístico 2024") +
  scale_fill_distiller(name = "IVU", palette = "YlOrRd", direction = 1) +
  theme_void() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, size = 13))

Realizamos una gráfica faceteada para conocer cuales son los barrios con mayor y menor índice de verde urbano en Valencia.

df_TOPIVU<- barrios_IVU_f %>% 
  arrange(desc(IVU)) %>% 
  slice(1:10) %>% 
  mutate(grupo = "Top 10 más Verde")

df_bottomIVU <- barrios_IVU_f %>%
  arrange(IVU) %>%
  slice(1:10) %>%
  mutate(grupo = "Top 10 menos Verde")

# 2) Unir y transformar en un solo df
df_facetIVU <- bind_rows(df_TOPIVU, df_bottomIVU)

# 3) Graficar faceteado
ggplot(df_facetIVU, aes(x = reorder(nombre, IVU),
                     y = IVU,
                     fill = IVU)) +
  geom_col() +
  coord_flip() +
  scale_fill_viridis_c(option = "plasma") +
  facet_wrap(~ grupo, ncol = 1, scales = "free_y") +
  labs(
    title = "Índice de Verde Urbano (Top 10 más y menos Verde)",
    x = "Barrio",
    y = "Densidad (m²/hab)",
    fill = "m²/hab"
  ) +
  theme_minimal() +
  theme(
    strip.text = element_text(face = "bold", size = 12)
  )

Como se observa en la gráfica, tenemos los 10 barrios con mayor y menor m2/hab de áreas verdes. El barrio que mayor Índice de Verde Urbano obtuvo fue El Botanic con 39.89 m2/hab, seguido de Trinitat con 37.65 m2/hab; mientras que La Xerea obtuvo 27.02 m2/hab siendo el último del top 10.

Por otro lado, dentro del top 10 con menor áreas verdes por habitante, se encuentra con 0.011 m2/hab, seguido de El Perellonet con 0.09 m2/hab; sieno el último del 10 con menor área verde fue En Corts que obtuvo 0.58 m2/hab.

Se aplicó una regresión lineal entre la Densidad Poblacional y el Índice de Verde Urbano para conocer su relación.

ggplot(barrios_IVU_f, aes(x = densidad_hab_km2, y = IVU)) +
  geom_point(alpha = 0.7, size = 3) +
  geom_smooth(method = "lm", se = FALSE, linetype = "dashed") +
  geom_hline(yintercept = 9, color = "darkgreen", linetype = "solid", size = 1) +
  labs( title = "Relación entre Densidad Poblacional y Verde Urbano",
        x = "Densidad Poblacional (hab/km²)",
        y = "m² de Áreas Verdes por Habitante",
        subtitle = "Tendencia: a mayor densidad, menor disponibilidad verde") +
  theme_minimal() +
  theme (plot.title = element_text(hjust = 0.5, size = 12, face = "bold"),
         plot.subtitle = element_text(hjust = 0.5, size = 10, face = "bold"))

En líneas generales, podemos ver que, existe una relación negativa entre la densidad poblacional y disponibilidad de verde urbano. A medida que la densidad poblacional aumenta, los m2 de verde urbano por habitante tienden a disminuir. Esto significa que, muchos barrios con densidades altas y muy altas tienen muy poca disponibilidad de áreas verdes.

Existen casos particulares (puntos aislados por encima de línea de tendencia), donde algunos barrios con densidad media y alta cumplen con el criterio de más de 9 m2/hab que menciona la Organización Mundial de la Salud.

Como último punto, los barrios con muy baja y baja densidad poblacional tienen buena oferta de áreas verdes por habitante. Sin embargo, la mayoría de los puntos se encuentran por debajo de la línea de tendencia (azul) y la línea que marca el criterio de la OMS (verde), indicando déficit de áreas verdes.

Para graficar esta relación en el mapa, se realizó un proceso de normalizado de valores y se crearon 3 rangos de relaciones.

#reescalado de valores
dens_verde <- barrios_IVU_f |> 
  mutate(dens_norm = scales::rescale(densidad_hab_km2),
         verde_norm = scales::rescale(IVU)
  ) |>
  #rangos por categorías
  mutate(
    dens = cut(dens_norm, breaks = c(0, .33, .66, 1), labels = c('Baja', 'Media', 'Alta')),
    verde = cut(verde_norm, breaks = c(0, .33, .66, 1), labels = c('Bajo', 'Medio', 'Alto')),
    combinacion = paste(dens, verde, sep = "-"))

#Paleta manual de colores bivariado 3x3 (9 colores)
paleta_bivariada <- c(
  "Baja-Bajo"   = "#c3c3be",
  "Media-Bajo"  = "#ace4e4",
  "Alta-Bajo"   = "#5ac8c8",
  "Baja-Medio"  = "#dfb0d6",
  "Media-Medio" = "#a5add3",
  "Alta-Medio"  = "#5698b9",
  "Baja-Alto"   = "#be64ac",
  "Media-Alto"  = "#8c62aa",
  "Alta-Alto"   = "#3b4994",
  "NA-NA"       = "#10100f"
)

#grafica del resultado
ggmap(mapa) +
  geom_sf(data = dens_verde, aes(fill = combinacion), color = "white",
          inherit.aes = FALSE) +
  scale_fill_manual(values = paleta_bivariada, drop = FALSE) +
  annotation_scale(location = "bl", width_hint = 0.3, style = "bar", text_col = "black",
                   unit_category = "metric", unit_value = 2500) + 
  annotation_north_arrow(location = "tr", which_north = "true", 
                         style =   north_arrow_nautical) +
  theme_void() +
  theme( legend.position = "right",
         plot.title = element_text(hjust = 0.5, size = 10, face = "bold"),
         plot.subtitle = element_text(hjust = 0.5, size = 10, face = "bold")) +
  labs(title = "Mapa Bivariado: Densidad poblacional vs Verde Urbano",
       subtitle = "Valencia, España",
       caption = "Fuente: OSM, GAD Valencia",
       fill = "nDensidad / - Verde")

Podemos observar que la mayoría de barrios que tienen una densidad baja en ambos casos se encuetran en la periferia de la urbe. Mientras los que tienen una densidad poblacional media y alta y un bajo índice de verde urbano tieden agruparse en el centro de la ciudad.

6. Conclusiones de la investigación

Como puntos importantes de la investigación, se puede concluir en que: