Natural alero

Análisis espacial del arbolado público en la ciudad de Montevideo

“El mejor momento para plantar un árbol fue hace 20 años. El segundo mejor momento es ahora.” - Proverbio chino

Introducción

Montevideo, con 530 km², es el departamento menos extenso de Uruguay y aún así el más poblado, con casi la mitad de su población -1,3 millones de personas- habitando la capital. Montevideo destaca por su rambla, su cerro, el candombe resonando por la mayoría de sus calles, y hay hasta quienes dicen que se lo reconoce por su calma y silencio.

Más allá de la idiosincracia citadina, la Intendencia de Montevideo se encuentra en un camino de recuperación de los espacios públicos como puntos de cohesión social, equidad y accesibilidad, al mismo tiempo que apuesta a la construcción de cultura que este tipo de espacios genera. Pero somos y nos construimos también en las calles.

Sobran motivos para dar cuenta de la importancia de los árboles en las vías públicas: quizás la mejora de la calidad de aire surja en primer lugar, la reducción de la temperatura ambiente y la contribución a la proliferación del ecosistema natural que hemos devastado con la irrupción del cemento. Sus beneficios también incluyen controlar la escorrentía del agua de lluvia al absorber parte de ella y reducir la velocidad a la que llega al suelo, lo que ayuda a prevenir inundaciones y a proteger la infraestructura urbana. Y después está la mejora estética y paisajística, y el beneficio psicológico y emocional que tiene la presencia de árboles en las personas. Numerosos estudios demuestran su impacto en la reducción del estrés, la mejora del estado de ánimo y la promoción de la salud mental. Y sobre esto, las autoridades públicas tienen injerencia.

Objetivo

El presente trabajo busca analizar la asociación entre la cantidad de árboles en espacios públicos y su distribución en los barrios de la zona urbana, teniendo en cuenta el área del barrio y la cantidad de personas que lo habitan. En definitiva, analizar la distribución del arbolado en la ciudad y si existe una correlación espacial en relación los barrios, su área y su cantidad de habitantes.

En Montevideo hay más de 200.000 árboles en el espacio común. Desde 2001 cuando se contabilizaron 210.717 ejemplares, el último censo del arbolado público, de 2008 -viejo, quizás, pero vamos, son árboles- recogió información de poco más del 65% de estos, unos 137.219. Los restantes no pudieron ser georreferenciados (32 %) o fueron extraídos o cortados. Usando este conjunto de datos de arbolado público se investigará si su disposición geoespacial puede proporcionar alguna pista sobre los lugares donde sería beneficioso llevar a cabo una intervención. Es importante mencionar que este trabajo es un ejercicio académico y que tiene sus limitaciones de complejidad, técnicas y de extensión.

Dataset

Se trabajará con tres dataset, dos de ellos espaciales:

  • El correspondiente al mapa digital que contiene la representación de los árboles del ornato público (no incluye los de parques y plazas ni los del área rural del departamento). Los datos provienen del sistema de arbolado que mantiene el Servicio de Arbolado Público de la Intendencia.

  • El correspondiente al mapa digital que contiene los límites correspondientes a los barrios de la ciudad de Montevideo según definición del Instituto Nacional de Estadística (INE). El mismo se actualiza coincidentemente con los censos nacionales y el actual corresponde al Censo 2011, disponible para descargar desde: https://www.gub.uy/instituto-nacional-estadistica/datos-y-estadisticas/estadisticas/mapas-vectoriales-ano-2011

Ambos dataset fueron extraídos del Sistema de Información Geográfica (SIG) de la Intendencia de Montevideo (IM): https://sig.montevideo.gub.uy/

Librerías

Usaremos las librerías tidyverse, ggpubr, RColorBrewer, sf, nngeo, gstat, skimr, spdep, tmap, units y spatstat.

library(tidyverse)
library(ggpubr)
library(RColorBrewer)
library(sf)
library(nngeo)
library(gstat) 
library(skimr)
library(spdep)
library(tmap)
library(units)
library(spatstat)
Ingreso de datos

Vamos a utilizar el dataset de arbolado público publicado por el SIG.

arbol <- st_read("MVD/v_sig_arboles/v_sig_arboles.shp")
## Reading layer `v_sig_arboles' from data source 
##   `/Users/virginia/Documents/BDIT/5. Geoestadística/Geoestadística/MVD/v_sig_arboles/v_sig_arboles.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 137219 features and 13 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 559099.3 ymin: 6134600 xmax: 588049.8 ymax: 6156332
## Projected CRS: WGS 84 / UTM zone 21S

Y el de barrios de Montevideo:

barrios <- st_read("MVD/v_sig_barrios/v_sig_barrios/v_sig_barrios.shp") 
## Reading layer `v_sig_barrios' from data source 
##   `/Users/virginia/Documents/BDIT/5. Geoestadística/Geoestadística/MVD/v_sig_barrios/v_sig_barrios/v_sig_barrios.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 62 features and 4 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 551982.7 ymin: 6133499 xmax: 589227.2 ymax: 6159810
## Projected CRS: WGS 84 / UTM zone 21S
barrios_2 <- st_read("MVD/v_sig_barrios/ine_barrios_mvd_nbi85.shp")
## Reading layer `ine_barrios_mvd_nbi85' from data source 
##   `/Users/virginia/Documents/BDIT/5. Geoestadística/Geoestadística/MVD/v_sig_barrios/ine_barrios_mvd_nbi85.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 62 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 551982.7 ymin: 6133499 xmax: 589227.2 ymax: 6159810
## CRS:           NA
# Extraer la columna "area" del dataset barrios_2
barrios_2 <- as.data.frame(barrios_2) %>% 
  select(NROBARRIO, AREA_KM)

# Realizar la unión basada en el número de barrio
barrios <- merge(barrios, barrios_2[, c("NROBARRIO", "AREA_KM")], by = "NROBARRIO", all.x = TRUE)

También el de población, para posteriormente unir el dato al de barrios:

poblacion <- st_read("MVD/VivHogPer/Marco2011_ZONA_Montevideo_VivHogPer.shp")
## Reading layer `Marco2011_ZONA_Montevideo_VivHogPer' from data source 
##   `/Users/virginia/Documents/BDIT/5. Geoestadística/Geoestadística/MVD/VivHogPer/Marco2011_ZONA_Montevideo_VivHogPer.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 13621 features and 56 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 551982.7 ymin: 6133499 xmax: 589227.2 ymax: 6159810
## Projected CRS: WGS 84 / UTM zone 21S
Exploración básica de los datos

Con la función summary() de rbase, vamos a ver las primeras filas del dataset.

summary(arbol) 
##      ARBOL          ORDINAL            DISTANCIA            CAP         
##  Min.   :     1   Length:137219      Min.   :  0.000   Min.   :  0.000  
##  1st Qu.: 49022   Class :character   1st Qu.:  2.000   1st Qu.:  0.630  
##  Median :119738   Mode  :character   Median :  2.800   Median :  1.270  
##  Mean   :113816                      Mean   :  2.938   Mean   :  1.216  
##  3rd Qu.:174262                      3rd Qu.:  3.400   3rd Qu.:  1.630  
##  Max.   :226249                      Max.   :751.000   Max.   :761.000  
##    ID_GENERO      ID_TIPOVEG            ALTURA           TRONCOS      
##  Min.   : 10.0   Length:137219      Min.   :  0.000   Min.   : 0.000  
##  1st Qu.: 35.0   Class :character   1st Qu.:  5.000   1st Qu.: 1.000  
##  Median : 51.0   Mode  :character   Median :  8.000   Median : 1.000  
##  Mean   : 59.5                      Mean   :  8.462   Mean   : 1.145  
##  3rd Qu.: 63.0                      3rd Qu.: 12.000   3rd Qu.: 1.000  
##  Max.   :999.0                      Max.   :122.000   Max.   :45.000  
##    INCLINACIO       NOM_CIENTI         NOM_COMUN           ID_ESPECIE   
##  Min.   : 0.0000   Length:137219      Length:137219      Min.   :  1.0  
##  1st Qu.: 0.0000   Class :character   Class :character   1st Qu.:  1.0  
##  Median : 0.0000   Mode  :character   Mode  :character   Median :  1.0  
##  Mean   : 0.9537                                         Mean   : 12.1  
##  3rd Qu.: 0.0000                                         3rd Qu.:  3.0  
##  Max.   :45.0000                                         Max.   :888.0  
##   DESC_TIPO                  geometry     
##  Length:137219      POINT        :137219  
##  Class :character   epsg:32721   :     0  
##  Mode  :character   +proj=utm ...:     0  
##                                           
##                                           
##

Vamos a extraer los datos de arbol sobre “ejemplares secos”.

arbol <- arbol %>% 
  select(ARBOL, ALTURA, NOM_CIENTI, NOM_COMUN, geometry) %>% 
  filter(!is.na(geometry), 
         !NOM_COMUN=="Ejemplar seco") #Eliminamos los 1.580 ejemplares secos

Ahora, vemos que arbol es un dataset con 135.639 registros y 5 variables:

names(arbol)
## [1] "ARBOL"      "ALTURA"     "NOM_CIENTI" "NOM_COMUN"  "geometry"
Variable Descripción
ARBOL Código de identificación de cada árbol
ALTURA Altura del árbol
NOM_CIENTI Nombre científico del árbol
NOM_COMUN Nombre “común” del árbol
geometry Georreferenciación del árbol

Con la función glimpse() de dplyr podemos ver una breve descripción de las variables del dataset, ver qué tipo de dato forman parte del dataset, etc.

glimpse(arbol)
## Rows: 135,639
## Columns: 5
## $ ARBOL      <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 11, 12, 13, 14, 15, 16, 17, 18, …
## $ ALTURA     <dbl> 14.0, 12.0, 12.0, 13.5, 12.0, 3.0, 0.6, 14.0, 13.0, 15.0, 1…
## $ NOM_CIENTI <chr> "Platanus x acerifolia", "Platanus x acerifolia", "Melia az…
## $ NOM_COMUN  <chr> "Platano", "Platano", "Paraiso", "Paraiso", "Fresno america…
## $ geometry   <POINT [m]> POINT (574038.2 6139000), POINT (574039.1 6138990), P…

barrioses un dataset que contiene los 62 barrios de Montevideo y 6 columnas, una de ellas correspondiente a la geometría. Nos vamos a quedar solo con las variables que nos interesan.

barrios <- barrios %>% 
  select(AREA_KM, BARRIO, NROBARRIO, geometry) %>% 
  filter(!is.na(geometry))

Ahora veamos el contenido de poblacion:

skim(poblacion)
## Warning: Couldn't find skimmers for class: sfc_MULTIPOLYGON, sfc; No
## user-defined `sfl` provided. Falling back to `character`.
Data summary
Name poblacion
Number of rows 13621
Number of columns 57
_______________________
Column type frequency:
character 10
numeric 47
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
CODCOMP_A 0 1 10 10 0 13608 0
NOMBDEPTO 0 1 10 10 0 1 0
NOMBLOC 0 1 5 10 0 2 0
CDEPTO_ISO 0 1 4 4 0 1 0
CLOC_ISO 0 1 7 7 0 2 0
BARRIO 15 1 2 25 0 64 0
MUNICIPIO 15 1 1 2 0 8 0
ZONA_LEGAL 15 1 5 5 0 18 0
ZONA_BASIC 15 1 1 2 0 4 0
geometry 0 1 152 45557 0 13621 0

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
AREA 0 1 38765.70 200251.87 1.508e+01 4.32969e+03 10192.18 17830.30 7298824.10 ▇▁▁▁▁
PERIMETER 0 1 597.86 783.95 1.627e+01 3.12710e+02 419.40 600.28 15216.27 ▇▁▁▁▁
DEPTO 0 1 1.00 0.00 1.000e+00 1.00000e+00 1.00 1.00 1.00 ▁▁▇▁▁
SECCION 0 1 19.92 20.99 0.000e+00 1.10000e+01 15.00 20.00 99.00 ▇▂▁▁▁
SEGMENTO 0 1 111.83 98.12 0.000e+00 2.80000e+01 101.00 204.00 633.00 ▇▃▁▁▁
ZONA 0 1 185.38 266.65 0.000e+00 6.00000e+00 14.00 402.00 990.00 ▇▁▂▁▁
LOCALIDAD 0 1 51.40 163.24 2.000e+01 2.00000e+01 20.00 20.00 900.00 ▇▁▁▁▁
CODSEC 0 1 119.92 20.99 1.000e+02 1.11000e+02 115.00 120.00 199.00 ▇▂▁▁▁
CODSEG 0 1 120030.70 20981.52 1.000e+05 1.11008e+05 115110.00 120039.00 199267.00 ▇▂▁▁▁
CODCOMP 0 1 120030887.24 20981524.61 1.000e+08 1.11008e+08 115110600.00 120039600.00 199267600.00 ▇▂▁▁▁
CODLOC 0 1 1051.40 163.24 1.020e+03 1.02000e+03 1020.00 1020.00 1900.00 ▇▁▁▁▁
codcomp_1 0 1 119920762.70 21345668.66 0.000e+00 1.11008e+08 115110600.00 120039600.00 199267600.00 ▁▁▇▂▁
CCZ 0 1 10.01 4.75 0.000e+00 7.00000e+00 10.00 13.00 18.00 ▃▃▆▇▅
AREA_1 0 1 38764.71 200252.04 0.000e+00 4.32909e+03 10191.82 17830.30 7298824.00 ▇▁▁▁▁
NROBARRIO 0 1 32.52 17.82 0.000e+00 1.70000e+01 32.00 48.00 63.00 ▅▆▇▃▇
V_TOT 0 1 38.22 46.84 0.000e+00 2.00000e+00 26.00 53.00 636.00 ▇▁▁▁▁
V_TOT_OC 0 1 34.65 41.62 0.000e+00 2.00000e+00 24.00 48.00 501.00 ▇▁▁▁▁
V_TOT_DES 0 1 3.56 6.82 0.000e+00 0.00000e+00 1.00 5.00 248.00 ▇▁▁▁▁
V_PAR 0 1 38.14 46.76 0.000e+00 2.00000e+00 26.00 53.00 635.00 ▇▁▁▁▁
V_PAR_OC 0 1 34.59 41.56 0.000e+00 2.00000e+00 24.00 48.00 501.00 ▇▁▁▁▁
V_PAR_MA 0 1 0.80 2.03 0.000e+00 0.00000e+00 0.00 1.00 43.00 ▇▁▁▁▁
V_PAR_DES 0 1 3.55 6.80 0.000e+00 0.00000e+00 1.00 5.00 247.00 ▇▁▁▁▁
V_COL 0 1 0.08 0.36 0.000e+00 0.00000e+00 0.00 0.00 7.00 ▇▁▁▁▁
V_COL_OC 0 1 0.07 0.32 0.000e+00 0.00000e+00 0.00 0.00 6.00 ▇▁▁▁▁
V_COL_DES 0 1 0.01 0.13 0.000e+00 0.00000e+00 0.00 0.00 3.00 ▇▁▁▁▁
H_TOT 0 1 35.83 43.07 0.000e+00 2.00000e+00 25.00 50.00 505.00 ▇▁▁▁▁
H_PAR 0 1 35.76 43.00 0.000e+00 2.00000e+00 25.00 50.00 505.00 ▇▁▁▁▁
H_PAR_MA 0 1 0.80 2.03 0.000e+00 0.00000e+00 0.00 1.00 43.00 ▇▁▁▁▁
H_COL 0 1 0.06 0.31 0.000e+00 0.00000e+00 0.00 0.00 5.00 ▇▁▁▁▁
P_TOT 0 1 96.82 106.92 0.000e+00 6.00000e+00 76.00 140.00 3381.00 ▇▁▁▁▁
P_TOT_HOM 0 1 45.06 55.02 0.000e+00 3.00000e+00 36.00 66.00 3381.00 ▇▁▁▁▁
P_TOT_MUJ 0 1 51.76 56.84 0.000e+00 2.00000e+00 40.00 75.00 590.00 ▇▁▁▁▁
P_TOT_PAR 0 1 95.42 101.99 0.000e+00 5.00000e+00 75.00 139.00 1042.00 ▇▁▁▁▁
P_HOM_PAR 0 1 44.33 46.53 0.000e+00 2.00000e+00 36.00 65.00 478.00 ▇▁▁▁▁
P_MUJ_PAR 0 1 51.09 55.98 0.000e+00 2.00000e+00 40.00 74.00 576.00 ▇▁▁▁▁
P_TOT_MA 0 1 1.95 5.53 0.000e+00 0.00000e+00 0.00 2.00 145.00 ▇▁▁▁▁
P_HOM_MA 0 1 0.92 2.72 0.000e+00 0.00000e+00 0.00 1.00 80.00 ▇▁▁▁▁
P_MUJ_MA 0 1 1.03 2.92 0.000e+00 0.00000e+00 0.00 1.00 65.00 ▇▁▁▁▁
P_TOT_COL 0 1 1.40 30.13 0.000e+00 0.00000e+00 0.00 0.00 3381.00 ▇▁▁▁▁
P_HOM_COL 0 1 0.73 29.36 0.000e+00 0.00000e+00 0.00 0.00 3381.00 ▇▁▁▁▁
P_MUJ_COL 0 1 0.67 5.55 0.000e+00 0.00000e+00 0.00 0.00 389.00 ▇▁▁▁▁
T_P_0_14 0 1 18.57 20.64 0.000e+00 0.00000e+00 14.00 27.00 284.00 ▇▁▁▁▁
T_P_15_64 0 1 63.33 74.12 0.000e+00 4.00000e+00 49.00 91.00 3334.00 ▇▁▁▁▁
T_P_65mas 0 1 14.92 20.50 0.000e+00 0.00000e+00 8.00 22.00 328.00 ▇▁▁▁▁
DENS_V 0 1 32.81 40.93 0.000e+00 2.90000e-01 25.58 43.70 616.44 ▇▁▁▁▁
DENS_H 0 1 31.00 37.98 0.000e+00 2.10000e-01 24.49 41.81 570.71 ▇▁▁▁▁
DENS_P 0 1 84.32 91.78 0.000e+00 6.20000e-01 71.48 121.33 1616.42 ▇▁▁▁▁

Vamos a tomar la población del dataset poblacion para asociarla a cada barrio.

# Calculamos la suma de población por barrio
poblacion_por_barrio <- poblacion %>%
  group_by(NROBARRIO) %>%
  summarize(P_TOT_Sum = sum(P_TOT)) 
  
# Sacamos la geometría
poblacion_por_barrio <- as.data.frame(poblacion_por_barrio) %>% 
  select(NROBARRIO, P_TOT_Sum)

# Realizamos la unión basada en el ID de barrio
barrios <- merge(barrios, poblacion_por_barrio[, c("NROBARRIO", "P_TOT_Sum")], by = "NROBARRIO", all.x = TRUE)
Visualización de distribuciones y relaciones

En primer lugar, vamos a visualizar los tipos de árboles que existen en Montevideo.

# Obtenemos una lista de los tipos de árboles únicos
tipos_arboles <- unique(arbol$NOM_COMUN)

# Ordenamos la lista alfabéticamente
tipos_arboles_ordenados <- sort(tipos_arboles)

Según el dataset, existen 309 tipos de árboles plantados en espacios públicos de Montevideo.

Ahora vamos a convertir arbol a dataframe.

arbol_df <- as.data.frame(arbol)
# Contamos la cantidad de árboles por tipo de árbol
arbol_df %>%
  group_by(NOM_COMUN) %>%
  summarise(cantidad_arboles = n()) %>%
  arrange(desc(cantidad_arboles)) %>%
  top_n(10) %>%  # Seleccionamos los 10 tipos de árboles más comunes
# Creamos el gráfico de barras
ggplot(aes(x = reorder(NOM_COMUN, cantidad_arboles), 
           y = cantidad_arboles)) +
  geom_bar(stat = "identity", fill = "forestgreen") +
  labs(x = "", 
       y = "Cantidad", 
       title = "10 tipos de árboles más comunes en Montevideo", 
       caption= "Fuente: Sistema de Información Geográfica (SIG)") +
  guides(fill=guide_legend(title.position = "top", ncol=1))+
  theme(plot.margin = margin(0.25, 1, 0.25, 0.1, "cm"), #ajustar los margenes del mapa
        panel.background = element_rect(fill = "gray100", colour = "gray100", linewidth = 2, linetype = "solid"), 
        panel.grid.major = element_line(linewidth = 0.5, linetype = "dashed", colour = "gray80"), 
        panel.grid.minor = element_line(linewidth = 0.25, linetype = "dashed", colour = "gray90"), 
        title=element_text(size=12, face = "bold"),
        plot.caption=element_text(face = "italic", colour = "gray35",size=6), 
        axis.text.x = element_text(hjust = 1)) +
  coord_flip()

Veamos cómo se distribuyen espacialmente.

# Crear el gráfico de dispersión
ggplot() +
  geom_sf(data = barrios, color="grey100")+
  geom_sf(data = arbol, color="forestgreen", alpha=0.1) +
  labs(title = "Distribución espacial del arbolado en Montevideo", 
       caption= "Fuente: Sistema de Información Geográfica (SIG)") +
  guides(fill=guide_legend(title.position = "top", ncol=1))+
  theme(plot.margin = margin(0.25, 1, 0.25, 0.1, "cm"), #ajustar los margenes del mapa
        panel.background = element_rect(fill = "gray100", colour = "gray100", linewidth = 2, linetype = "solid"), 
        panel.grid.major = element_line(linewidth = 0.5, linetype = "dashed", colour = "gray80"), 
        panel.grid.minor = element_line(linewidth = 0.25, linetype = "dashed", colour = "gray90"), 
        title=element_text(size=12, face = "bold"),
        plot.caption=element_text(face = "italic", colour = "gray35",size=6), 
        axis.text.x = element_text(hjust = 1)) +
  guides(color = FALSE)  # Desactivar la guía de color
## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Vemos que hay una concentración sobre el centro y sur este del departamento, la llamada área urbana del departamento.

Y ahora vamos a construir dos indicadores de interés: cantidad de árboles por área de cada barrio (arbol_km2) y proporción de árboles por habitante (arbol_hab), también al interior del barrio.

arbol_barrio <- st_intersection(barrios,arbol) %>% 
  group_by(BARRIO) %>% 
    summarise(cant_arboles=sum(n())) 
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
arbol_barrio <- st_join(barrios,arbol_barrio)

arbol_barrio <- arbol_barrio %>% 
  mutate(arbol_km2 = cant_arboles/AREA_KM, # Cantidad de árboles por km_2 
         arbol_hab = cant_arboles/P_TOT_Sum) # Proporción de árboles por habitante

Y veamos ahora los barrios con mayor cantidad de árboles.

ggplot() +
  geom_sf(data = arbol_barrio, aes(fill = arbol_km2)) +
  geom_sf(data = barrios, fill=NA)+
  labs(title = "Proporción de árboles por km2, por barrio", 
       caption= "Fuente: Sistema de Información Geográfica (SIG)") +
  guides(fill = guide_colorbar(title = "Árboles por km²", fill=guide_legend(title.position = "top", ncol=1)))+
  theme(plot.margin = margin(0.25, 1, 0.25, 0.1, "cm"), #ajustar los margenes del mapa
        panel.background = element_rect(fill = "gray100", colour = "gray100", linewidth = 2, linetype = "solid"), 
        panel.grid.major = element_line(linewidth = 0.5, linetype = "dashed", colour = "gray80"), 
        panel.grid.minor = element_line(linewidth = 0.25, linetype = "dashed", colour = "gray90"), 
        title=element_text(size=12, face = "bold"),
        plot.caption=element_text(face = "italic", colour = "gray35",size=6), 
        axis.text.x = element_text(hjust = 1)) +
  scale_fill_gradient(low="#edf8e9", high="forestgreen")

Parece haber una concentración de densidad de árboles al centro de la ciudad.

ggplot() +
  geom_sf(data = arbol_barrio, aes(fill = arbol_hab)) +
  geom_sf(data = barrios, fill=NA)+
  labs(title = "Proporción de árboles por habitantes, por barrio", 
       caption= "Fuente: Sistema de Información Geográfica (SIG)") +
  guides(fill = guide_colorbar(title = "Árboles por hab.", fill=guide_legend(title.position = "top", ncol=1)))+
  theme(plot.margin = margin(0.25, 1, 0.25, 0.1, "cm"), #ajustar los margenes del mapa
        panel.background = element_rect(fill = "gray100", colour = "gray100", linewidth = 2, linetype = "solid"), 
        panel.grid.major = element_line(linewidth = 0.5, linetype = "dashed", colour = "gray80"), 
        panel.grid.minor = element_line(linewidth = 0.25, linetype = "dashed", colour = "gray90"), 
        title=element_text(size=12, face = "bold"),
        plot.caption=element_text(face = "italic", colour = "gray35",size=6), 
        axis.text.x = element_text(hjust = 1)) +
  scale_fill_gradient(low="#edf8e9", high="forestgreen")

Destaca, por sobre todo, la concentración de árboles que existe en el barrio Carrasco.

Exploración de datos espaciales, de manera interactiva.

#Activamos modo "view"
tmap_mode("view")

#Chequeamos que no haya errores
tmap_options(check.and.fix = TRUE)

#Visualizamos
tm_shape(arbol_barrio) + 
  tm_fill("arbol_km2", 
          style = "quantile", 
          n = 10, 
          palette = "Greens",
          popup.vars = c("BARRIO.x",
                         "cant_arboles",
                         "AREA_KM",
                         "arbol_km2")) +
  tm_borders(alpha = 0.1) +
  tm_layout(main.title = "Cantidad de árboles por barrio", 
            main.title.size = 0.7,
            legend.position = c("right", "bottom"),
            legend.title.size = 0.8)
#Visualizamos
tm_shape(arbol_barrio) + 
  tm_fill("arbol_hab", 
          style = "quantile", 
          n = 10, 
          palette = "Greens",
          popup.vars = c("BARRIO.x",
                         "cant_arboles",
                         "P_TOT_Sum",
                         "arbol_hab")) +
  tm_borders(alpha = 0.1) +
  tm_layout(main.title = "Cantidad de árboles", 
            main.title.size = 0.7,
            legend.position = c("right", "bottom"),
            legend.title.size = 0.8)

Análisis de asociación espacial

1 - Lista con vecinos

Generamos los datos vecinos (clase nb) usando poly2nb().

w <- poly2nb(arbol_barrio, 
             row.names = "NROBARRIO", 
             queen = FALSE) # Seteo queen = false para que se necesite más de un punto de contacto para configurarse como vecinos.  

class(w) # Clase del componente
## [1] "nb"

Exploramos los vecinos

w
## Neighbour list object:
## Number of regions: 62 
## Number of nonzero links: 322 
## Percentage nonzero weights: 8.376691 
## Average number of links: 5.193548

Mapeamos cómo son estas relaciones entre vecinos.

plot(st_geometry(arbol_barrio), 
     border="grey80", 
     main = paste0("Relaciones entre vecinos"))
plot(w, 
     coords =  st_coordinates(st_centroid(arbol_barrio)),
     add = TRUE,
     col="darkorchid", 
     )

2 - Asignación de pesos a los vecinos

Existen diferentes posibilidades para “ponderar” (asignar pesos) a cada relación con los vecinos. Algunas de las opciones incluye hacer una clasificación “binaria” (con style='B') o “estandarizada” (con style='W').

lwb <-  nb2listw(w, style='B')
lwb
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 62 
## Number of nonzero links: 322 
## Percentage nonzero weights: 8.376691 
## Average number of links: 5.193548 
## 
## Weights style: B 
## Weights constants summary:
##    n   nn  S0  S1   S2
## B 62 3844 322 644 7144

Hay 322 relaciones de primer grado y 644 de segundo grado.

3 - Calcular I Moran (Global)
# Retenemos sólo el Índice de Moran
I_MORAN <- moran(arbol_barrio$arbol_hab, #Variable
                 listw = lwb, #Lista de vecinos
                 n = length(w), #Cantidad de polígonos
                 S0 = Szero(lwb)) [1] #Suma total de los pesos

print(I_MORAN)
## $I
## [1] 0.3533156

El valor del índice de Moran es de 0.35 lo cual no representa un valor alto para este tipo de prueba. Podría existir cierta tendencia de agrupamiento de valores similares en el espacio.

4 - Testear significatividad

A continuación realizamos una prueba para testear si el I de Moran obtenido es significativo.

moran.test(arbol_barrio$arbol_hab, # Data
           lwb) # Pesos
## 
##  Moran I test under randomisation
## 
## data:  arbol_barrio$arbol_hab  
## weights: lwb    
## 
## Moran I statistic standard deviate = 4.9821, p-value = 3.145e-07
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.353315637      -0.016393443       0.005506761

Con un p-valor tan bajo, no podría descartarse la hipótesis de no existe autocorrelación espacial en la cantidad de árboles por habitante.

5 - Correlograma de Moran
I_CORR <- sp.correlogram(neighbours = w, # Vecinos
                         var = arbol_barrio$arbol_hab, # Variable de interés
                         order = 3,  # La cantidad de órdenes que vamos a estudiar
                         method = "I", # Método: I de Moran
                         style = "B", # B corresponde a Binaria
                         zero.policy = TRUE)

I_CORR
## Spatial correlogram for arbol_barrio$arbol_hab 
## method: Moran's I
##          estimate expectation   variance standard deviate Pr(I) two sided    
## 1 (62)  0.3533156  -0.0163934  0.0055068           4.9821        6.29e-07 ***
## 2 (62) -0.0136653  -0.0163934  0.0027357           0.0522       0.9584015    
## 3 (62) -0.1657481  -0.0163934  0.0020396          -3.3071       0.0009427 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(I_CORR, 
     main = paste0("Correlograma de Moran"))

Lo que nos indica el correlograma de Moran es que el único orden que puede ser significativo es el primero, es decir los vecinos directos. Los vecinos de los vecinos no tendría sentido.

6 - Diagrama de dispersión de Moran
mp <- moran.plot(arbol_barrio$arbol_hab,  # Variable de interés
                 listw = lwb,  # Lista de vecinos
                 labels = arbol_barrio$BARRIO, 
                 main = paste0("Diagrama de dispersión")) 

mp
##              x        wx is_inf labels        dfb.1_         dfb.x
## 1  0.007566706 0.4135600   TRUE      1  1.383449e-02 -0.0119479540
## 2  0.073508137 0.2935707  FALSE      2 -1.297112e-01  0.0789483203
## 3  0.061397262 0.1675896  FALSE      3 -2.293867e-01  0.1573001380
## 4  0.087785001 1.0116196   TRUE      4  2.349478e-01 -0.1117467739
## 5  0.086514778 0.2498467  FALSE      5 -1.381163e-01  0.0676914366
## 6  0.100664400 0.3389300  FALSE      6 -7.950224e-02  0.0223887534
## 7  0.081841115 0.1834536  FALSE      7 -1.816450e-01  0.0977881099
## 8  0.082789152 0.5433438  FALSE      8  9.055549e-03 -0.0047917479
## 9  0.128925888 0.5002068  FALSE      9 -1.328803e-02 -0.0166725350
## 10 0.144127371 0.6191785  FALSE     10 -3.550071e-05 -0.0023142820
## 11 0.138410503 0.4243477  FALSE     11 -1.013708e-02 -0.0493302026
## 12 0.011634996 0.5443369  FALSE     12  1.569134e-01 -0.1341952676
## 13 0.233817238 0.6893816   TRUE     13  6.716737e-02 -0.1049477634
## 14 0.328125000 0.4066938   TRUE     14  1.066855e+00 -1.4346091810
## 15 0.172876527 0.6476384  FALSE     15  6.252074e-03 -0.0154955750
## 16 0.035726594 0.4673375  FALSE     16  1.762797e-02 -0.0139701366
## 17 0.103786277 0.3782798  FALSE     17 -6.201614e-02  0.0134073041
## 18 0.099020221 0.5126959  FALSE     18 -1.504853e-02  0.0047022353
## 19 0.049969531 0.7952724  FALSE     19  2.628805e-01 -0.1946399221
## 20 0.042309547 0.1388883  FALSE     20 -2.824373e-01  0.2175404057
## 21 0.088896768 0.4562123  FALSE     21 -3.688814e-02  0.0170573814
## 22 0.139906417 0.4845145  FALSE     22 -5.398524e-03 -0.0363953006
## 23 0.123244734 0.8677526  FALSE     23  4.959300e-02  0.0307869467
## 24 0.094410361 0.7424031  FALSE     24  8.290115e-02 -0.0321849925
## 25 0.146223992 0.9937601  FALSE     25 -4.731501e-03  0.1223169975
## 26 0.130221456 0.6202970  FALSE     26  1.980878e-03  0.0029097920
## 27 0.167790471 0.6249471  FALSE     27  7.888841e-03 -0.0218795723
## 28 0.071965722 0.5943425  FALSE     28  5.246438e-02 -0.0325222647
## 29 0.081323225 1.0228571   TRUE     29  2.751362e-01 -0.1494692096
## 30 0.014508463 0.5838156  FALSE     30  1.927041e-01 -0.1635975266
## 31 0.053018885 0.2377501  FALSE     31 -1.860506e-01  0.1352880335
## 32 0.029046175 0.3638863  FALSE     32 -7.255458e-02  0.0589512098
## 33 0.162242620 0.6650766  FALSE     33 -1.262214e-03  0.0041530798
## 34 0.132556850 0.6677714  FALSE     34  6.059505e-03  0.0119410591
## 35 0.178551177 0.3820554  FALSE     35  8.754680e-02 -0.1973248645
## 36 0.022569166 0.2753347  FALSE     36 -1.592052e-01  0.1321115424
## 37 0.032599401 0.4545602  FALSE     37  1.088912e-02 -0.0087348373
## 38 0.210573343 0.7363150  FALSE     38 -1.694792e-04  0.0002901371
## 39 0.180553493 1.0758991  FALSE     39 -1.212922e-01  0.2658326447
## 40 0.141832900 0.7118773  FALSE     40  2.120955e-03  0.0255435422
## 41 0.136821685 0.6464375  FALSE     41  2.078364e-03  0.0075957977
## 42 0.176362089 0.9771359  FALSE     42 -8.069814e-02  0.1881088906
## 43 0.184500988 0.7929872  FALSE     43 -3.376742e-02  0.0704816226
## 44 0.170411771 0.7709939  FALSE     44 -2.189127e-02  0.0570675536
## 45 0.174044176 1.0144062  FALSE     45 -8.540649e-02  0.2071254848
## 46 0.126977055 0.7793586  FALSE     46  2.671643e-02  0.0264879109
## 47 0.111666667 0.4791162  FALSE     47 -2.578410e-02 -0.0001020081
## 48 0.159382663 0.7180573  FALSE     48 -8.147369e-03  0.0302208170
## 49 0.143044390 0.8033677  FALSE     49  2.465404e-03  0.0535240854
## 50 0.084767048 0.6136005  FALSE     50  4.216360e-02 -0.0214623711
## 51 0.180704259 0.9242946  FALSE     51 -7.303332e-02  0.1597430837
## 52 0.151579091 0.5341649  FALSE     52  6.024415e-03 -0.0399031721
## 53 0.045671539 0.6158990  FALSE     53  1.269288e-01 -0.0961750143
## 54 0.119063943 0.9398757  FALSE     54  7.689820e-02  0.0256667656
## 55 0.072942122 0.3778326  FALSE     55 -7.899820e-02  0.0484124418
## 56 0.116313528 0.6137300  FALSE     56  7.715034e-03  0.0014740874
## 57 0.064184163 0.3072819  FALSE     57 -1.272453e-01  0.0852614548
## 58 0.154441880 0.3697327  FALSE     58  2.257486e-02 -0.1131544595
## 59 0.074629970 0.6177490  FALSE     59  6.165131e-02 -0.0370004469
## 60 0.133499688 0.3389276  FALSE     60 -2.542551e-02 -0.0567618190
## 61 0.014264900 0.2971122  FALSE     61 -1.309178e-01  0.1112144718
## 62 0.048114247 0.1108384  FALSE     62 -2.980198e-01  0.2229387728
##            dffit     cov.r       cook.d        hat
## 1   0.0138407492 1.1039657 9.740193e-05 0.06329864
## 2  -0.1488462240 1.0241514 1.108406e-02 0.02244277
## 3  -0.2472079760 0.9881469 2.995982e-02 0.02710245
## 4   0.3069581297 0.8963022 4.418553e-02 0.01859317
## 5  -0.1777876534 0.9976934 1.563639e-02 0.01886360
## 6  -0.1270274995 1.0185156 8.074288e-03 0.01664614
## 7  -0.2227448242 0.9733180 2.422870e-02 0.01997981
## 8   0.0112061123 1.0547863 6.384594e-05 0.01973798
## 9  -0.0606860849 1.0451865 1.866050e-03 0.01744582
## 10 -0.0049007320 1.0560679 1.221189e-05 0.02075816
## 11 -0.1220895068 1.0279169 7.483058e-03 0.01927595
## 12  0.1570915358 1.0854688 1.246586e-02 0.05967994
## 13 -0.1172099469 1.1198576 6.967189e-03 0.08134073
## 14 -1.4900804172 1.0341155 9.965871e-01 0.22073952
## 15 -0.0218215458 1.0684509 2.420673e-04 0.03253449
## 16  0.0179072604 1.0783637 1.630319e-04 0.04121020
## 17 -0.1058330440 1.0278694 5.631091e-03 0.01639210
## 18 -0.0232954692 1.0507413 2.757901e-04 0.01681411
## 19  0.2735321827 0.9931414 3.666732e-02 0.03267272
## 20 -0.2895136231 0.9988353 4.110154e-02 0.03704428
## 21 -0.0488558922 1.0489961 1.211059e-03 0.01836803
## 22 -0.0861178727 1.0417830 3.747476e-03 0.01963625
## 23  0.1630337396 0.9982021 1.316654e-02 0.01672546
## 24  0.1186186797 1.0247414 7.059425e-03 0.01741083
## 25  0.2469349977 0.9634322 2.960428e-02 0.02137322
## 26  0.0099136994 1.0525712 4.996898e-05 0.01764953
## 27 -0.0322233299 1.0648742 5.276700e-04 0.02992611
## 28  0.0596089690 1.0530906 1.802109e-03 0.02296511
## 29  0.3357647624 0.8833189 5.244288e-02 0.02011522
## 30  0.1930593464 1.0744555 1.875656e-02 0.05721104
## 31 -0.1950759165 1.0256834 1.896833e-02 0.03107490
## 32 -0.0732313390 1.0797612 2.721711e-03 0.04582432
## 33  0.0064858657 1.0631991 2.138918e-05 0.02733825
## 34  0.0365717493 1.0506074 6.792436e-04 0.01805372
## 35 -0.2664735400 1.0060340 3.496939e-02 0.03571111
## 36 -0.1600108954 1.0718837 1.291373e-02 0.05066970
## 37  0.0110250811 1.0809187 6.180350e-05 0.04332165
## 38  0.0003404728 1.0989146 5.894326e-08 0.05890277
## 39  0.3543214769 0.9637801 6.047694e-02 0.03689905
## 40  0.0573011268 1.0497350 1.665024e-03 0.02012902
## 41  0.0197930171 1.0533982 1.991332e-04 0.01891465
## 42  0.2579386653 1.0064661 3.279362e-02 0.03445239
## 43  0.0917557103 1.0690776 4.266040e-03 0.03934348
## 44  0.0820506226 1.0600234 3.411136e-03 0.03124226
## 45  0.2889941512 0.9865578 4.078360e-02 0.03316517
## 46  0.1077225241 1.0289467 5.834711e-03 0.01716698
## 47 -0.0540607040 1.0448159 1.481571e-03 0.01612909
## 48  0.0488804893 1.0587102 1.213067e-03 0.02610909
## 49  0.1163826103 1.0329458 6.812351e-03 0.02045549
## 50  0.0532404784 1.0493528 1.437778e-03 0.01925870
## 51  0.2127146525 1.0322771 2.255685e-02 0.03698991
## 52 -0.0725722561 1.0507183 2.667936e-03 0.02311825
## 53  0.1308785109 1.0548409 8.640717e-03 0.03506255
## 54  0.2093383076 0.9635092 2.133094e-02 0.01637520
## 55 -0.0903153749 1.0456105 4.122944e-03 0.02263204
## 56  0.0188817732 1.0504754 1.812157e-04 0.01622794
## 57 -0.1387398528 1.0361353 9.668937e-03 0.02591683
## 58 -0.1963129648 1.0059864 1.909214e-02 0.02415376
## 59  0.0712889062 1.0495070 2.574300e-03 0.02207590
## 60 -0.1671878071 1.0016271 1.385916e-02 0.01823039
## 61 -0.1311502256 1.0867561 8.704299e-03 0.05741751
## 62 -0.3088119963 0.9774850 4.634179e-02 0.03368455
mp %>% 
  filter(is_inf)
##              x        wx is_inf labels     dfb.1_       dfb.x       dffit
## 1  0.007566706 0.4135600   TRUE      1 0.01383449 -0.01194795  0.01384075
## 4  0.087785001 1.0116196   TRUE      4 0.23494780 -0.11174677  0.30695813
## 13 0.233817238 0.6893816   TRUE     13 0.06716737 -0.10494776 -0.11720995
## 14 0.328125000 0.4066938   TRUE     14 1.06685537 -1.43460918 -1.49008042
## 29 0.081323225 1.0228571   TRUE     29 0.27513624 -0.14946921  0.33576476
##        cov.r       cook.d        hat
## 1  1.1039657 9.740193e-05 0.06329864
## 4  0.8963022 4.418553e-02 0.01859317
## 13 1.1198576 6.967189e-03 0.08134073
## 14 1.0341155 9.965871e-01 0.22073952
## 29 0.8833189 5.244288e-02 0.02011522

Los barrios que parecen influir en esta investigación son 1 (Ciudad Vieja), 4 (Cordón), 13 (Punta Gorda), 14 (Carrasco) y 29 (Aires Puros).

7 - C de Geary Global
geary.test(arbol_barrio$arbol_hab, lwb)
## 
##  Geary C test under randomisation
## 
## data:  arbol_barrio$arbol_hab 
## weights: lwb 
## 
## Geary C statistic standard deviate = 5.439, p-value = 2.679e-08
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##       0.490257895       1.000000000       0.008783379

El C de Geary es de 0.49. Si bien inferior a 1, no próximo a 0. Esto confirma que no podría descartarse la hipótesis de no existe autocorrelación espacial en la cantidad de árboles por habitante.

IcorrC <-sp.correlogram(neighbours=w,
                       var=arbol_barrio$arbol_hab,
                       order=5,
                       method="C")
plot(IcorrC, 
     main = "Correlograma de Geary")

A nivel global, los valores del arbol por habitante tienen vecinos con valores similares.

8 - Moran local (\(I_i\))

Vamos a calcular el índice de Moran local:

LOC_MORAN <- localmoran(arbol_barrio$arbol_hab, 
                        listw = lwb)
Local Indicators of Spatial Association (LISA)
arbol_barrio <- cbind(arbol_barrio, # Base original
             mp[c("x", "wx")], # Variable y lag (estandarizado)
             LOC_MORAN,  # Valores de Moran Local
             attributes(LOC_MORAN)$quadr) %>% # Cuadrantes LISA
  # Renombramos la columna "Pr(z != E(Ii))"
  rename(p = Pr.z....E.Ii..) %>%  
  # Los valores no significativos se diferencian en otra categoría
  mutate(quad = ifelse(p > 0.05, 5, mean),
         quad = factor(quad, 
                       levels = 1:5,
                       labels = c("Low-Low", 
                                  "High-Low", 
                                  "Low-High", 
                                  "High-High", 
                                  "No Signif")))

table(arbol_barrio$quad)
## 
##   Low-Low  High-Low  Low-High High-High No Signif 
##         4         1         0         4        53

Definimos los colores para LISA.

LISA_col <- c("blue2","skyblue1", "lightpink", "red2", "white")
names(LISA_col) <- levels(arbol_barrio$quad)
LISA_col
##     Low-Low    High-Low    Low-High   High-High   No Signif 
##     "blue2"  "skyblue1" "lightpink"      "red2"     "white"

Mostramos estos resultados del hotspots en el scatTerplot.

arbol_barrio %>%  
  st_drop_geometry() %>%  
  ggplot(aes(x = x, y = wx)) + 
  geom_hline(linetype = 2, yintercept = mean(arbol_barrio$wx)) +
  geom_vline(linetype = 2, xintercept = mean(arbol_barrio$x)) +
  geom_point(aes(fill = quad), shape = 21) + 
  geom_smooth(method = lm, se = F, linetype = 2, color = "darkorchid" ) + 
  labs(x = "Variable", 
       y = "Lag",
       title = "Scatterplot de LISA") +
    theme(plot.margin = margin(0.25, 1, 0.25, 0.1, "cm"), #ajustar los margenes del mapa
        panel.background = element_rect(fill = "gray100", colour = "gray100", linewidth = 2, linetype = "solid"), 
        panel.grid.major = element_line(linewidth = 0.5, linetype = "dashed", colour = "gray80"), 
        panel.grid.minor = element_line(linewidth = 0.25, linetype = "dashed", colour = "gray90"), 
        title=element_text(size=12, face = "bold"),
        plot.caption=element_text(face = "italic", colour = "gray35",size=6), 
        axis.text.x = element_text(hjust = 1)) +
  scale_fill_manual(values = LISA_col, drop = F)
## `geom_smooth()` using formula = 'y ~ x'

Por último visualizamos estos resultados en el mapa:

tmap_mode("view")
tm_shape(arbol_barrio) + 
  tm_fill(col = "quad", #"mean"
          alpha = 0.6,
          palette = LISA_col, 
          style = "fixed",
          title="Agrupamientos LISA", 
          popup.vars = c("BARRIO.x",
                         "cant_arboles",
                         "P_TOT_Sum",
                         "arbol_hab")) +
  tm_legend(outside = TRUE) +
  tm_borders(col = "Grey")

Del mapa podemos ver que el área pintada en celeste, correspondiente al barrio Carrasco, es un área con alta proporción de árboles por habitantes, rodeada por áreas que presentan una baja proporción del índice analizado.

Vale comentar que Carrasco es el barrio de mayor riqueza acumulada en Montevideo.

alto_alto <-arbol_barrio %>% 
  filter(quad=="High-High") %>% 
   as.data.frame(.) %>% 
   select(BARRIO.x)

knitr::kable(
  alto_alto, caption = 'Barrios con clusters de arbolado Alto-Alto'
)
Barrios con clusters de arbolado Alto-Alto
BARRIO.x
PUNTA GORDA
PRADO NUEVA SAVONA
REDUCTO
FIGURITA 
bajo_bajo <-arbol_barrio %>% 
  filter(quad=="Low-Low") %>% 
   as.data.frame(.) %>% 
   select(BARRIO.x)

knitr::kable(
  bajo_bajo, caption = 'Barrios con clusters de arbolado Bajo-Bajo'
)
Barrios con clusters de arbolado Bajo-Bajo
BARRIO.x
JARDINES DEL HIPODROMO
PIEDRAS BLANCAS
VILLA GARCIA MANGA RURAL
MANGA

En tanto los barrios pintados en rojo (Punta Gorda, Prado, Reducto y Figurita) son aquellos con mayores niveles de proporción de árboles por cantidad de habitantes, rodeados por vecinos con también valores altos. Mientras que Jardines del Hipódromo, Piedras Blancas, Villa García y Manga son los barrios con bajo indicador rodeados de otros también con bajos valores. Parecería que sería un área indicada para una política de recuperación del arbolado en espacios públicos.

Análisis de procesos puntuales

Si bien se agruparon por barrios para tener una mirada por área, los árboles representan patrones de puntos, por tanto procederemos a analizar si su frecuencia de ocurrencia es aleatoria o responde a algún patrón.

1 - Definición de ventana

Definimos la ventana (clase owin) en la cual vamos a trabajar.

MVD <- summarise(barrios) # Junta todo en un solo objeto del tamaño de la ciudad
MVD <- nngeo::st_remove_holes(MVD) # Remueve agujeros

Ventana <- as.owin(st_geometry(MVD)) # Me quedo con la geometría de la ciudad y lo convierte a ventana con as.owin()
unitname(Ventana) <- "Meter" # Indicamos unidad de medida utilizada

plot(Ventana)

2 - Generación de patrón de puntos

Generamos el patrón de puntos (clase ppp).

arbol_ppp <- as.ppp(st_geometry(arbol), W = Ventana)
arbol_ppp
## Planar point pattern: 135639 points
## window: polygonal boundary
## enclosing rectangle: [551982.7, 589227.2] x [6133499, 6159810] Meter

Para evitar los puntos duplicados movemos de forma aleatoria los puntos (con rjitter()).

set.seed(200)
arbol_ppp <- rjitter(arbol_ppp, 
                     retry=TRUE, 
                     nsim = 1, 
                     radius = 2, 
                     drop = TRUE)

cat("\nExisten puntos duplicados:", any(duplicated.ppp(arbol_ppp)))
## 
## Existen puntos duplicados: FALSE
plot(arbol_ppp, 
     cols = "forestgreen", 
     main = "Dispersión espacial de arbolado montevideano")

Miremos al interior del objeto ppp.

summary(arbol_ppp)
## Planar point pattern:  135639 points
## Average intensity 0.0002568665 points per square Meter
## 
## Coordinates are given to 1 decimal place
## i.e. rounded to the nearest multiple of 0.1 Meter
## 
## Window: polygonal boundary
## single connected closed polygon with 18193 vertices
## enclosing rectangle: [551982.7, 589227.2] x [6133499, 6159810] Meter
##                      (37240 x 26310 Meter)
## Window area = 528053000 square Meter
## Unit of length: 1 Meter
## Fraction of frame area: 0.539
3 - Estadísticas de resumen de primer orden
3.1 - Intensidad del patrón de puntos

Lo primero que podemos saber es cuál es la intensidad del patrón de árboles (cuantos árboles se han plantado por unidad de área).

cat("Intensidad de árboles por metro cuadrado:", 
    intensity(arbol_ppp)  %>%  
      format(scientific = F), "\n")
## Intensidad de árboles por metro cuadrado: 0.0002568665
# Lo pasamos a km para visualizar mejor
cat("Intensidad de árboles por kilómetro cuadrado:", 
    format( (intensity(arbol_ppp) * (1000^2) ), 
            scientific = F, digits = 4, 
            decimal.mark = ",") )
## Intensidad de árboles por kilómetro cuadrado: 256,9

Esto asume que la intensidad del patrón de puntos es homogénea, es decir, que en todos lados hay la misma cantidad de árboles, pero claramente esto no es así.

3.2 - Intensidad por cuadrantes

Una forma de saber si existe un patrón aleatorio es subdividir la ventana en porciones que cubran toda el área de interés (teselación) y calcular la intensidad en cada una de estas áreas.

Optamos por la teselación por hexágonos.

# Generamos hexágonos con lado 1500m
H <- hextess(Ventana, 1500)

QC_arbol_HEX <- quadratcount(arbol_ppp, tess = H) # H tiene que ser del tipo tess, que se genera con la función hextess()

plot(arbol_ppp, 
     main = "Conteo por cuadrante (hexagonal)",
     cols = "forestgreen")
plot(QC_arbol_HEX, add = TRUE, cex = 1)

intensity(QC_arbol_HEX, 
          image = T)  %>% # image es para poder generar la imagen 
  plot(main = "Intensidad por cuadrante (hexagonal)")

3.3 - Test de hipótesis de Poisson

Vamos a evaluar si se trata de un patrón homogéneo de Poisson (CSR) mediante un test de hipótesis:

  • Hipótesis nula (\(H_0\)): la intensidad es homogénea y la distribución de casos en los cuadrantes responde a una distribución de Poisson (CSR)

  • Hipótesis alternativa (\(H_1\)) : la intensidad no es homogénea (en una forma no especificada).

Si consideramos una significación de \(0.05\), ¿el patrón de distribución espacial es aleatorio? La forma por defecto que calcula la significancia está basada en \(\chi^2\), lo cual supone algunas restricciones: (a) el numero de conteos por cuadrado tiene que ser > 5; y (b) Cada uno de los cuadrantes es relativamente similar en superficie.

Dado que no todos los cuadrantes son de igual superficie y además algunos conteos son menores a 5 usamos el método Montecarlo.

quadrat.test(QC_arbol_HEX, 
             method = "MonteCarlo", 
             nsim = 1000)
## 
##  Conditional Monte Carlo test of CSR using quadrat counts
##  Test statistic: Pearson X2 statistic
## 
## data:  
## X2 = 393129, p-value = 0.001998
## alternative hypothesis: two.sided
## 
## Quadrats: 122 tiles (irregular windows)

Existe una probilidad de 99.8% de que rechazar la hipótesis nula esté bien, y los datos no sigan una distribución al azar.

4 - Estadísticas de resumen de segundo orden

Ademas de preguntar si el patrón de puntos es significativamente distinto de la aleatoriedad espacial completa (considerando distintos rangos de distancia o escalas espaciales), estos indicadores permiten determinar si los puntos tienden a formar agregaciones o presentar dispersión (regularidad) en un rango de distancias.

Como el tamaño de la muestra dificulta correr simulaciones, se evaluará la concentración de puntos a partir de la densidad.

Intensidad de Kernel

La función density() genera una imagen “raster” con los valores estimados de la densidad para cada punto de la ventana.

DD <- density(arbol_ppp, 
              edge = F) # No tiene en cuenta el efecto del borde

plot(DD, 
     main='Densidad de arboles en Montevideo')

Conclusiones

Después de realizar los análisis de procesos puntuales y de asociación espacial, parecería que el proceso de arbolado de Montevideo no obedece un orden aleatorio y que habría áreas de prioridad para llevar a cabo políticas de forestación urbana. Tal como se comentó anteriormente, hay ciertas áreas que parecerían prioritarias y estás están explicitadas en el análisis LISA, en aquel cluster identificado como “bajo-bajo”.

Por otro lado, surge del análisis que Ciudad Vieja, si bien no conforma un cluster podría identificarse como un barrio para llevar a cabo una intervención de forestación ya que presenta valores bajos de los dos índices estudiados: arbol_hab y arbol_km2.

Referencias bibliográficas

Bozzo, A. et. al. (2021). Consultoría para apoyo al desarrollo e implementación de planes de arbolado urbano y áreas verdes departamentales considerando la capacidad de producción de viveros, para la adaptación al cambio climático y la variabilidad. Montevideo, Uruguay.

Gabinete Ambiental (2017). Política Nacional de Cambio Climático.

Ochoa de la Torre, J.M. (2010). Ciudad, vegetación e impacto climático. El confort en los espacios urbanos. Palapa: Universidad de Colima, México

Virginia Recagno - Trabajo final de Geoestadística para el posgrado Big Data e Inteligencia Territorial de FLACSO Argentina