Diversidad del arbolado viario de Barcelona
Juan Carlos López
2023-11-12
Para contactar con el autor, jcar2lopez[gmail] o https://es.linkedin.com/in/juan-carlos-l%C3%B3pez-l%C3%B3pez-11542959
Objetivos, resumen y conclusiones
La introducción de la librería sf en R ha cambiado el manejo y análisis de datos espaciales en comparación con librerías anteriores como sp y rgdal. Su ventaja principal es el enfoque en objetos espaciales que combinan datos geográficos y atributos en una única estructura. Esto facilita la manipulación, visualización y análisis espacial, con una sintaxis más intuitiva y clara, operaciones espaciales más sencillas y una mayor compatibilidad con otras librerías de análisis de datos en R.
Para este ejercicio, utilizando datos del arbolado viario de Barcelona (disponibles Open Data), se va a calcular la diversidad biológica (índice de biodiversidad específica de Shannon) en una cuadrícula superpuesta sobre el mapa de la ciudad. Esto va más allá de la mera representación de puntos sobre un territorio, pues se trata de realizar operaciones más o menos complejas sobre una segmentación de dicho territorio, generada para responder a los intereses del investigador.
Una primera fase es la creación de los objetos espaciales necesarios. Se parte de un mapa de Barcelona en formato shp, que se lee de forma muy directa, con la función st_read(). El siguiente paso es descargar los datos del arbolado, disponibles en formato csv y convertir la tabla resultante en un objeto sf con la función st_as_sf(). Finalmente se crea una cuadrícula circunscrita al área del mapa de Barcelona, utilizando la función st_make_grid(). Una vez definidas sus dimensiones, se obtiene una cuadricula de 25 x 25 = 625 celdas, con un área aproximada de 630 m2 cada celda.
Una segunda fase se orienta a las operaciones sobre los puntos contenidos en cada celda de la cuadrícula. Mediante la función st_intersects() , usando como argumentos la cuadrícula y la tabla de los árboles, se determina a que celda pertenece cada árbol. Dicha información se incorpora a la tabla de los árboles y se puede utilizar para segmentar los análisis que se quiera aplicar. En este ejercicio se va a calcular el índice de biodiversidad de Shannon (H’), utilizando la función diversity() de la librería vegan, que considera tanto el número de especies como los ejemplares de cada una de ellas. El valor de este índice varía entre 0,5 y 5, considerándose balos los valores inferiores a 2 y altos los superiores a 3. Los resultados se incorporan a la cuadrícula mediante la función left_join() y se representan gráficamente utilizando ggplot() y geom_sf().
La mayor parte de Barcelona cuenta con una diversidad biológica del arbolado viario que se puede considerar media, con valores comprendidos entre 2 y 3, y algunas pocas zonas con valores entre 3 y 4 que se pueden considerar altos. La menor diversidad, por debajo de 2, se corresponde con el distrito del Eixample y algunos puntos de la periferia urbana.
Más allá del caso concreto objeto de este desarrollo, el procedimiento aquí descrito tiene aplicaciones para el análisis espacial de datos geolocalizados de cualquier naturaleza.
Presentación
Trabajo habitualmente con datos de salud, pero recientemente sentí la necesidad de hacer algo diferente. Por ello he preparado este ejercicio que combina dos temas que me atraen: los mapas y la ecología numérica. En cierta medida es una vuelta a mis orígenes, pues mi primer desarrollo fue una plantilla configurada en una hoja de cálculo de una aplicación llamada Lotus 1 2 31.
Utilizando datos del arbolado viario de Barcelona, disponibles en su magnífico portal de Open Data, se va a calcular la diversidad biológica en una cuadrícula superpuesta sobre el mapa de la ciudad.
Esto va más allá de la mera representación de puntos sobre un territorio, pues se trata de realizar operaciones más o menos complejas sobre una segmentación de dicho territorio, generada para responder a los intereses del investigador.
La librería sf
La idea inicial vino de la mano de una entrada del inspirador blog long time ago…, concretamente de la entrada “Are betting houses randomly distributed in Madrid? Point pattern analysis in R”. Lo estaba replicando con datos del arbolado, cuando decidí replantear el desarrollo para hacerlo exclusivamente utilizando la librería sf.
La introducción de la librería sf en R ha revolucionado el manejo y análisis de datos espaciales en comparación con librerías anteriores como sp y rgdal. La ventaja principal de sf es su enfoque en objetos espaciales simples que combinan datos geográficos y atributos en una estructura única. Esto facilita la manipulación, visualización y análisis espacial, con una sintaxis más clara, operaciones espaciales más sencillas y mayor compatibilidad con otras librerías de R.
Preparación del entorno
El primer paso es definir el directorio de trabajo (wd)
setwd("C:/wd/abol_viario_bcn")Cargar librerías
library(sf)
library(readxl)
library(readr)
library(tidyverse)
library(kableExtra)
library(units)
library(vegan)
library(plotly)Mapa de Barcelona
Se parte de un mapa de Barcelona con sus unidades administrativas
(distritos), que se puede descargar desde el portal Open Data en este enlace. Se lee, de forma muy directa, con la función
st_read(), de la librería sf.
mapa.sf <- st_read("data/0301040100_Districtes_UNITATS_ADM.shp", stringsAsFactors = FALSE)## Reading layer `0301040100_Districtes_UNITATS_ADM' from data source
## `C:\wd\abol_viario_bcn\data\0301040100_Districtes_UNITATS_ADM.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 10 features and 48 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 420812.5 ymin: 4574282 xmax: 435480.4 ymax: 4591066
## Projected CRS: ETRS89 / UTM zone 31NEl mapa así creado tiene las coordenadas en sistema UTM. Para seguir
adelante con el resto del ejercicio, es necesario transformarlas al
sistema LONLAT2. Se utiliza la función
st_transform(), que necesita el argumento
crs, con la información del sistema de referencia de
coordenadas. En este caso, para asegurar la correcta superposición del
mapa con el resto de representaciones espaciales que se van a generar,
se va a utilizar el sistema de coordenadas mundial WGS 843 en todos los objetos
espaciales de este documento.
mapa.sf <- st_transform(mapa.sf, "+proj=longlat +ellps=WGS84 +datum=WGS84")Se puede ya disponer de una vista del objeto creado.
mapa.sf %>% ggplot() + geom_sf()
Datos del arbolado
La información se puede descargar desde este enlace, perteneciente al portal Open Data del ayuntamiento de Barcelona.
arb <- readr::read_csv("data/2022_1T_OD_Arbrat_Viari_BCN.csv", show_col_types = FALSE)Un primer paso es seleccionar las columnas de interés, en este caso la latitud y longitud, el nombre científico y el código del distrito.
arb <- arb %>% dplyr::select(latitud, longitud, nom_cientific, codi_districte)Esta es la cabecera de la tabla resultante.
kbl(head(arb)) %>%
kable_paper(bootstrap_options = "striped", full_width = F, position = "left")| latitud | longitud | nom_cientific | codi_districte |
|---|---|---|---|
| 41.43729 | 2.165353 | Fraxinus angustifolia ‘Raywood’ | 08 |
| 41.43734 | 2.165436 | Fraxinus angustifolia ‘Raywood’ | 08 |
| 41.43777 | 2.162530 | Platanus x hispanica | 07 |
| 41.43779 | 2.162501 | Platanus x hispanica | 07 |
| 41.43781 | 2.162471 | Platanus x hispanica | 07 |
| 41.43787 | 2.165110 | Pinus pinea | 08 |
Se eliminan los posibles registros incompletos y ya se pueden obtener algunos datos resumen acerca del arbolado viario de Barcelona.
arb <- na.omit(arb)Por ejemplo, el número de especies diferentes.
length(unique(arb$nom_cientific))## [1] 267Y el número de ejemplares diferentes.
nrow(arb)## [1] 154226El siguiente paso es convertir el dataframe arb en puntos
espaciales con la librería sf. Los argumentos principales son
coords y crs.
coords necesita un vector con las coordenadas
longitudinales y latitudinales. Por defecto, se espera que las columnas
se llamen “long” y “lat”, respectivamente pero, como en este caso, se
pueden especificar nombres diferentes.
Para el argumento CRS se utiliza el definido para el mapa de Barcelona
arb.sf <- st_as_sf(x = arb,
coords = c("longitud", "latitud"),
crs = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")Se puede dar un vistazo a la tabla así creada. Es muy sencilla, pues solo contiene el nombre científico del ejemplar, el código del distrito donde se asienta y la columna geometry, que contiene las coordenadas de su localización.
kbl(head(arb.sf)) %>%
kable_paper(bootstrap_options = "striped", full_width = F, position = "left")| nom_cientific | codi_districte | geometry |
|---|---|---|
| Fraxinus angustifolia ‘Raywood’ | 08 | POINT (2.165353 41.43729) |
| Fraxinus angustifolia ‘Raywood’ | 08 | POINT (2.165436 41.43734) |
| Platanus x hispanica | 07 | POINT (2.16253 41.43777) |
| Platanus x hispanica | 07 | POINT (2.162501 41.43779) |
| Platanus x hispanica | 07 | POINT (2.162471 41.43781) |
| Pinus pinea | 08 | POINT (2.16511 41.43787) |
Existen algunas funciones que proporcionan información que puede ser
interesante sobre los objetos espaciales. class() informa
de la naturaleza del objeto.
class(arb.sf) ## [1] "sf" "tbl_df" "tbl" "data.frame"st_bbox() devuelve los límites del polígono que lo
encierra.
st_bbox(arb.sf) ## xmin ymin xmax ymax
## 2.088837 41.344797 2.226232 41.466312st_crs() proporciona el CRS (coordinate reference
system) utilizado.
st_crs(arb.sf)[1]## $input
## [1] "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"Ahora se puede representar una muestra de los árboles sobre el mapa de Barcelona.
ggplot() +
geom_sf(data = mapa.sf) +
# geom_sf(data = grid.sf, alpha = 0.2) +
geom_sf(data = arb.sf[sample(1:nrow(arb.sf), 1000), ],
aes(color = codi_districte), size = 0.3) +
theme_minimal() +
labs(title = "Arbolado viario de Barcelona",
subtitle = "Muestra de 1.000 ejemplares")
Cuadrícula sobre el mapa de Barcelona
Para crear una cuadricula embebida en el mapa de Barcelona
directamente como un objeto sf, se utilizará la función
st_make_grid()4, si bien primero hay que definir las
dimensiones de la cuadrícula. Utilizaremos una cuadrícula de 25 celdas
de lado, y con ese datos, se establecen las dimensiones de cada celda
(cellsize).
cellsize = c(diff(st_bbox(mapa.sf)[c(1, 3)]), diff(st_bbox(mapa.sf)[c(2, 4)]))/25La expresión anterior calcula el tamaño de celda dividiendo las
diferencias entre las coordenadas x e y del cuadro delimitador,
obtenidas con la función st_bbox(), del
objeto mapa.sf por 25, el número de celdas que va
a utilizar..
Para crear la cuadricula se emplea la función
st_make_grid(), que genera una retícula de polígonos
que cubra el área definida por el objeto
mapa.sf, y con un tamaño de celda definido
por el argumento cellsize. El argumento
what = "polygons" indica que la cuadrícula
debe estar compuesta por polígonos en lugar de puntos o líneas, y el
argumento square = TRUE establece que los
polígonos generados sean cuadrados.
area.grid <- st_make_grid(mapa.sf,
cellsize = cellsize,
what = "polygons",
square = TRUE)El objeto que se obtiene pertenece a la clase “sfc_POLYGON”, una geometría simple, y conviene transformarlo en un data frame espacial para operar posteriormente con el. Además, se va a añadir un valor id a cada celda.
grid.sf <- st_sf(area.grid) %>%
mutate(id = 1:length(lengths(area.grid)))Se examina el objeto así creado.
head(grid.sf)## Simple feature collection with 6 features and 1 field
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 2.052333 ymin: 41.31704 xmax: 2.094504 ymax: 41.32309
## CRS: +proj=longlat +ellps=WGS84 +datum=WGS84
## area.grid id
## 1 POLYGON ((2.052333 41.31704... 1
## 2 POLYGON ((2.059362 41.31704... 2
## 3 POLYGON ((2.06639 41.31704,... 3
## 4 POLYGON ((2.073419 41.31704... 4
## 5 POLYGON ((2.080447 41.31704... 5
## 6 POLYGON ((2.087476 41.31704... 6Aparte de la descripción del objeto y su proyección, tenemos una tabla con solo dos columnas. En area.grid aparece la relación de polígonos que conforman las celdas. Para cada uno de ellos se definen 5 pares de puntos, el primer par es igual al último, y representan las cuatros esquinas de la celda. La columna id contiene el identificador del polígono, que será más adelante de vital importancia para operar con los datos contenidos en la celda.
Es posible obtener información sobre las características de las celdas. Extraemos la primera celda.
celda1 <- grid.sf[[1]][1]Mediante la función st_area() se puede
conocer su tamaño en metros cuadrados:
st_area(celda1) ## 394886.4 [m^2]Es fácil convertir esta medida Km2 usando la librería
units:
units::set_units(st_area(celda1), km^2)## 0.3948864 [km^2]Y la longitud aproximada del lado.
sqrt(st_area(celda1)) ## 628.3999 [m]Ya están todos los elementos disponibles. Es hora de representar de forma conjunta el mapa, la cuadricula y una muestra de los árboles. El orden en el que se van a disponer las diferentes capas de la figura es importante, para evitar que queden ocultas bajo la capa siguiente.
ggplot() +
geom_sf(data = mapa.sf) +
geom_sf(data = grid.sf, alpha = 0.2) +
geom_sf(data = arb.sf[sample(1:nrow(arb.sf), 1000), ],
aes(color = codi_districte), size = 0.3) +
theme_minimal() +
labs(title = "Arbolado viario de Barcelona",
subtitle = "Muestra de 1.000 ejemplares")
Contar puntos en la cuadrícula
El análisis ecológico puede empezar contando los ejemplares (puntos)
que hay en cada celda5. Para ello, se va a utilizar la función
st_intersect(), que determina cuantos
puntos hay en la intersección entre los 2 objetos espaciales. Genera una
lista con 625 elementos, el número de celdas, y para cada una de ellas
informa del número de árboles existentes 6. En la tabla
resultante, area.grid es la geometría de la celda, id
es el identificador de la celda y n_arb el número de ejemplares
que contiene.
grid.sf$n_arb <- lengths(st_intersects(grid.sf, arb.sf))Esta es la tabla resultante, ordenada de forma descendente a partir del número de árboles en la celda. Se observa que la número 447 es la que contiene un mayor número de ejemplares.
grid.sf %>%
arrange(-n_arb) %>%
head() %>%
kable("html", align=c("l", 'c', 'c','c')) %>%
kable_styling(full_width = FALSE, position="left")%>%
row_spec(0, background = "#F3E2A9")| id | n_arb | area.grid |
|---|---|---|
| 447 | 2447 | POLYGON ((2.199931 41.41989… |
| 423 | 2060 | POLYGON ((2.20696 41.41384,… |
| 421 | 1806 | POLYGON ((2.192903 41.41384… |
| 468 | 1626 | POLYGON ((2.171817 41.42595… |
| 312 | 1600 | POLYGON ((2.129646 41.38964… |
| 422 | 1576 | POLYGON ((2.199931 41.41384… |
Para visualizar los resultados, se representan sobre el mapa de Barcelona.
ggplot() +
geom_sf(data = grid.sf, aes(fill = n_arb)) +
scale_fill_viridis_c(option = "B", direction = -1,
legend_title <- "Nº de ejemplares") +
geom_sf(data = mapa.sf, fill = NA, colour = "black", size = 1) +
theme(panel.background = element_rect(fill = "transparent"),
rect = element_rect(fill = "transparent")) +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank() ) +
ggtitle("Densidad del arbolado viario en BCN")
Recortar la cuadrícula
El mapa anterior es suficiente para fines científicos, pero la estética se puede mejorar. Un recurso habitual es eliminar las celdas sin datos.
grid.sf <- filter(grid.sf, n_arb > 0)El objeto grid.sf ha pasado de 625 a 217 elementos. Se representa de nuevo el mapa, asociando cada celda con el número de árboles que contiene, haciendo transparente el fondo (fill = NA) y eliminando las etiquetas de los ejes.
ggplot() +
geom_sf(data = grid.sf, aes(fill = n_arb)) +
scale_fill_viridis_c(option = "B", direction = -1,
legend_title <- "Nº de ejemplares") +
geom_sf(data = mapa.sf, fill = NA, colour = "black", size = 1) +
theme(panel.background = element_rect(fill = "transparent"),
rect = element_rect(fill = "transparent")) +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank() ) +
ggtitle("Densidad del arbolado viario en BCN")
La densidad del arbolado no es homogénea en Barcelona. Las mayores densidades aparecen al norte del distrito de Sant Martí, al sur de Sant Andreu y en la frontera de Les Corts con el Eixample. Sin embargo, como se verá más adelante, más densidad de ejemplares no implica necesariamente una mayor biodiversidad.
Determinar a que celda pertenece cada árbol.
Se utiliza la función st_intersects(),
que utiliza como argumentos el conjunto de puntos espaciales y la
cuadrícula. Se muestran los primeros 3 registros.
int <- sf::st_intersects(arb.sf , grid.sf)
int[1:3]## [[1]]
## [1] 183
##
## [[2]]
## [1] 183
##
## [[3]]
## [1] 182El objeto int es una lista, que contiene tantos elementos como árboles, y cada uno tiene asociado el valor id de la celda a la que pertenece. Los tres primeros elementos de la lista son los árboles 1, 2, y 3, que se encuentran en las celdas 183, 183, y 182 respectivamente. Ahora hay que añadir esa información (la celda a la que pertenece cada árbol) al conjunto de puntos arb.sf.
arb.sf$celda <- grid.sf$id[unlist(int)]Volviendo a la cabecera de la tabla arb.sf, se comprueba que hay una nueva columna, celda, que contiene el identificador de la celda en la que se encuentra el árbol.
(head(arb.sf)) %>%
kable("html", align=c("l", 'c', 'c','c')) %>%
kable_styling(full_width = FALSE, position="left")%>%
row_spec(0, background = "#F3E2A9")| nom_cientific | codi_districte | geometry | celda |
|---|---|---|---|
| Fraxinus angustifolia ‘Raywood’ | 08 | POINT (2.165353 41.43729) | 492 |
| Fraxinus angustifolia ‘Raywood’ | 08 | POINT (2.165436 41.43734) | 492 |
| Platanus x hispanica | 07 | POINT (2.16253 41.43777) | 491 |
| Platanus x hispanica | 07 | POINT (2.162501 41.43779) | 491 |
| Platanus x hispanica | 07 | POINT (2.162471 41.43781) | 491 |
| Pinus pinea | 08 | POINT (2.16511 41.43787) | 492 |
Análisis de la vegetación
¿Alguna vez se ha preguntado alguien si la diversidad del arbolado viario de Barcelona es similar en toda la ciudad? Es posible, pero una manera de resolver la pregunta es calculando en cada celda el índice de diversidad de Shannon-Weaver.
El índice de Shannon, de Shannon-Weaver o de Shannon-Wiener (Según la Wikipedia7 ) se usa en ecología para medir la biodiversidad específica. El índice contempla la cantidad de especies presentes en el área de estudio (riqueza de especies), y la cantidad relativa de individuos de cada una de esas especies (abundancia).
Este índice se representa normalmente como H’ y se expresa con un número positivo, que en la mayoría de los ecosistemas naturales varía entre 0,5 y 5, aunque su valor normal está entre 2 y 3; valores inferiores a 2 se consideran bajos en diversidad y superiores a 3 son altos en diversidad de especies. Los ecosistemas con mayores valores son los bosques tropicales y arrecifes de coral, y los menores las zonas desérticas.
\[ H'=-\sum_{i=1}^S p_i \text{ } \text{ln } p_i \] donde:
S es el número de especies
pi es la proporción de individuos de la especie i respecto al total de individuos (es decir la abundancia relativa de la especie i, ni/N)
ni es el número de individuos de la especie i
N es el número de todos los individuos de todas las especies
Para calcular el índice de diversidad (H’) se va a utilizar la
función diversity() de la librería
vegan8 9.
Para aplicar la función diversity() es
necesario realizar algunas transformaciones a los datos, pues necesita,
entre otras cosas, que las estaciones (los puntos o zonas de
donde proceden los datos) sean las filas, las especies sean las columnas
y la intersección el número de ejemplares10.
| Celda | Especie 1 | Especie 2 | Especie 3 | Especie 4 |
|---|---|---|---|---|
| Celda 1 | X | X | X | X |
| Celda 2 | X | X | X | X |
| Celda 3 | X | X | X | X |
| Celda 4 | X | X | X | X |
| Celda 5 | X | X | X | X |
Partiendo de arb.sf se crea un nuevo
conjunto de datos, data, seleccionando solo las
columnas de interés (la celda y el nombre científico). Para eliminar la
columna geometry, la manera más cómoda es usar la función
st_drop_geometry() de la librería
sf.
data <- arb.sf %>%
select(celda, nom_cientific) %>%
st_drop_geometry()Posteriormente se cuentan los ejemplares de cada especie presentes en cada celda.
La tabla tiene este aspecto.
head(data) %>%
kable("html", align=c('c', 'l','r')) %>%
kable_styling(full_width = FALSE, position="left")%>%
row_spec(0, background = "#F3E2A9")| celda | nom_cientific | n |
|---|---|---|
| 111 | Casuarina cunninghamiana | 11 |
| 111 | Ulmus pumila | 26 |
| 112 | Casuarina cunninghamiana | 7 |
| 112 | Melia azedarach | 14 |
| 112 | Ulmus pumila | 23 |
| 112 | Washingtonia robusta | 3 |
El siguiente paso es pasar de formato “long” (largo) a formato “wide” (ancho). Las celdas son ahora las filas y las especies las columnas.
data <- data %>%
pivot_wider(names_from = nom_cientific, values_from = n)Este es el aspecto de la tabla. La celda 111 contiene 11 ejemplares de Casuarina cunninghamiana, 26 de Ulmus pumilla …
head(data[, 1:4]) %>%
kable("html", align=c("l", 'c', 'c','c','c')) %>%
kable_styling(full_width = FALSE, position="left")%>%
row_spec(0, background = "#F3E2A9")| celda | Casuarina cunninghamiana | Ulmus pumila | Melia azedarach |
|---|---|---|---|
| 111 | 11 | 26 | NA |
| 112 | 7 | 23 | 14 |
| 137 | 3 | 19 | NA |
| 138 | NA | 52 | NA |
| 139 | NA | 3 | NA |
| 140 | NA | NA | NA |
La función para el cálculo del índice de diversidad necesita operar sobre valores numéricos, por lo que hay que reemplazar los valores NA por ceros.
data <- data %>% replace(is.na(.), 0)El siguiente paso es extraer un vector con los identificadores de las celdas y convertir el objeto data en una matriz, manteniendo los identificadores como nombres de filas
Extraer vector con los identificadores de celdas (celdas). Se muestran los primeros 5 elementos.
celdas <- data$celda
celdas[1:5]## [1] 111 112 137 138 139Hay que onvertir el objeto data en una matriz (data2), pero conteniendo solo las especies, es decir, eliminando la primera columna, la que contiene los nombres de las celdas.
data2 <- as.matrix(data[,2:nrow(data)])Este es el aspecto de la matriz.
head(data2[, 1:3]) ## Casuarina cunninghamiana Ulmus pumila Melia azedarach
## [1,] 11 26 0
## [2,] 7 23 14
## [3,] 3 19 0
## [4,] 0 52 0
## [5,] 0 3 0
## [6,] 0 0 0Ahora toca renombrar las filas, asignándoles el identificador de la celda.
rownames(data2) <- celdasY el aspecto ha cambiado a este.
head(data2[, 1:3])## Casuarina cunninghamiana Ulmus pumila Melia azedarach
## 111 11 26 0
## 112 7 23 14
## 137 3 19 0
## 138 0 52 0
## 139 0 3 0
## 140 0 0 0Ya están los datos preparados en el formato que necesita la función
diversity(). En los documentos (entre
otros) “Trabajando con vegan”11 , “Analysis of
community ecology data in R”12 y en el propio manual de la librería, se
ofrece información sobre esta función y las opciones de sus
argumentos.
Se aplica la función de diversidad y se muestran los primeros 5 valores,
diversidad <- diversity(data2, index = "shannon", MARGIN = 1, base = exp(1))
diversidad[1:5]## 111 112 137 138 139
## 0.6085569 1.1697140 1.7597453 2.6253226 2.9967163“diversidad” es un “named vector”, (un vector que tiene asignado un nombre a cada uno de sus elementos) que contiene los valores del índice en cada celda. Para terminar la preparación de los datos, hay que convertirlo en un data frame y añadirle el identificador de la celda.
diversidad <- as.data.frame(diversidad)
diversidad$celda <- rownames(diversidad)Este es el resultado.
head(diversidad) %>%
select(celda, diversidad) %>%
kable("html", align=c("c", 'c', 'c','c','c')) %>%
kable_styling(full_width = FALSE, position="left")%>%
row_spec(0, background = "#F3E2A9")| celda | diversidad | |
|---|---|---|
| 111 | 111 | 0.6085569 |
| 112 | 112 | 1.1697140 |
| 137 | 137 | 1.7597453 |
| 138 | 138 | 2.6253226 |
| 139 | 139 | 2.9967163 |
| 140 | 140 | 1.5326599 |
Hay que cambiar el formato de “celda” a integer, para hacerlo compatible con el que tiene en grid.sf.
diversidad$celda<- as.integer(diversidad$celda)Para terminar, se añade esta información al objeto grid.sf.
grid.sf <- left_join(grid.sf, diversidad, by = c("id" = "celda"))Este es el resultado.
head(grid.sf) %>%
kable("html", align=c("c", 'c', 'c','c','c')) %>%
kable_styling(full_width = FALSE, position="left")%>%
row_spec(0, background = "#F3E2A9")| id | n_arb | diversidad | area.grid |
|---|---|---|---|
| 111 | 37 | 0.6085569 | POLYGON ((2.122618 41.34124… |
| 112 | 47 | 1.1697140 | POLYGON ((2.129646 41.34124… |
| 137 | 76 | 1.7597453 | POLYGON ((2.129646 41.34729… |
| 138 | 786 | 2.6253226 | POLYGON ((2.136675 41.34729… |
| 139 | 199 | 2.9967163 | POLYGON ((2.143703 41.34729… |
| 140 | 57 | 1.5326599 | POLYGON ((2.150732 41.34729… |
Es posible ahora examinar como es la distribución de la diversidad del arbolado de Barcelona.
summary(diversidad$diversidad)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 1.760 2.250 2.121 2.649 3.323p1<- ggplotly(
diversidad %>% arrange(diversidad) %>% mutate(acum = cumsum(diversidad)) %>%
ggplot() +
geom_histogram((aes(x = diversidad)),
bins = 11, color = "darkblue",
fill = "lightblue",
size = 0.5) +
theme_minimal() +
ggtitle("Distribución de la diversidad")
)
p2 <-ggplotly(
diversidad %>% arrange(diversidad) %>%
mutate(cumsum = cumsum(diversidad)/sum(diversidad)) %>%
ggplot(aes(x= diversidad, y=cumsum)) +
geom_line(color = "blue") +
ggtitle("Distribución acumulada") +
theme_minimal()
)
fig <- subplot(p1, p2, margin = 0.05) %>%
layout(title = 'Distribución de la diversidad')
figEl rango es bastante amplio, de 0 a 3.3, con una media y mediana similares, entorno al 2,2, valor que corresponde a una diversidad media. A grosso modo, se puede estimar que un 20% de la superficie de Barcelona presenta una diversidad que se puede considerar baja, mientras que un 4% muestra una alta diversidad. En conjunto, y sin pretender haber realizado un revisión bibliográfica exahustiva de este tema, se puede señalar que la diversidad del arbolado de Barcelona es algo superior a la que presentan otras ciudades similares donde se han realizado este tipo de análisis13.
Por último, se puede analizar, aunque muy someramente, la posible relación entre la densidad arbórea y la diversidad. un primer paso es representar gráficamente dicha relación.
# left_join(diversidad, grid.sf, by = c("celda" = "id")) %>%
grid.sf %>%
ggplot(aes(x = n_arb, y = diversidad)) +
geom_point(color = "darkorange", size = 2) +
geom_smooth(method="nls", se=FALSE, formula=y~a*log(x)+k,
method.args=list(start=c(a=1, k=1)), linewidth = 2) +
labs(x = "Árboles por celda",
y = "Diveridad biológica") +
ggtitle(("Relación diversidad - densidad arbórea")) +
theme_minimal()
A la figura anterior se le ha añadido la linea de tendencia correspondiente a una formulación logarítmica de la relación entre ambas variables, dado que la distribución de los puntos así lo sugiere: un cambio importante en el inicio de la línea (incremento rápido de la diversidad a medida que aumenta el número de ejemplares), seguido de una fase de estabilización (aunque el número de ejemplares sigue aumentando, la diversidad se mantiene bastante estable).
Otra discusión sería si la relación planteada es significativa. Los resultados de un primer modelo parecen apuntar en esa dirección. Sin embargo, es algo que debe tomarse con cautela; y al fin y al cabo, el análisis en profundidad de estos hallazgos no constituye el objetivo primordial del presente trabajo.
Representación de los resultados
Ya solo queda representar sobre el mapa de Barcelona los resultados obtenidos. Es necesario considerar que el orden de las capas y la asignación de transparencias son críticas para lograr una correcta visualización.
ggplot() +
geom_sf(data = grid.sf, aes(fill = diversidad)) +
scale_fill_viridis_c(option = "G", direction = -1,
alpha = 0.7,
legend_title <- "H") +
geom_sf(data = mapa.sf, fill = NA, colour = "black", size = 1.5) +
theme(panel.background = element_rect(fill = "transparent"),
rect = element_rect(fill = "transparent")) +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank() ) +
ggtitle("Diversidad del arbolado viario en BCN")
Una opción es convertir los valores de la diversidad en categorías.
Una forma sencilla es mediante la función
cut(). Es importante añadir el argumento
include.lowest = TRUE), porque puede haber celdas con
valor cero (todos los ejemplares pertenecen a la misma especie).
grid.sf <- grid.sf %>%
mutate(cat_diversidad = cut(diversidad,
breaks = seq(0, 3.5, by = 0.5),
labels = c("0-0.5", "0.5-1", "1-1.5", "1.5-2",
"2-2.5", "2.5-3", "3-3.5"),
include.lowest = TRUE))Se representan los resultados:
ggplot() +
geom_sf(data = grid.sf, aes(fill = cat_diversidad)) +
scale_fill_brewer(palette = "Greens", "H'") +
geom_sf(data = mapa.sf, fill = NA, colour = "black", size = 1) +
theme(panel.background = element_rect(fill = "transparent"),
rect = element_rect(fill = "transparent")) +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank() ) +
ggtitle("Diversidad del arbolado viario en BCN")
Otra posibilidad, interesante en función del tipo de comunicación que se vaya a efectuar, es etiquetar los resultados, para facilitar su interpretación a un público no especializado.
# Crear las etiquetas
grid.sf <- grid.sf %>%
mutate(cat_diversidad_2 = cut(diversidad,
breaks = seq(0, 4, by = 1),
labels = c("Muy baja", "Baja",
"Media", "Alta"),
include.lowest = TRUE))
# Representar
ggplot() +
geom_sf(data = grid.sf, aes(fill = cat_diversidad_2)) +
scale_fill_manual(values=c('khaki','gold','darkolivegreen2','chartreuse4')) +
geom_sf(data = mapa.sf, fill = NA, colour = "black", size = 1) +
theme(panel.background = element_rect(fill = "transparent"),
rect = element_rect(fill = "transparent")) +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank() ) +
labs(fill='Diversidad') +
ggtitle("Diversidad del arbolado viario en BCN")
Si bien en conjunto, Barcelona se caracteriza por tener una diversidad biológica del arbolado viario que se puede considerar media - alta, existe una cierta heterogeneidad espacial en su distribución. Tal vez lo más llamativo sea el caso del distrito del Eixample, situado en el centro de Barcelona, con una diversidad claramente inferior.
ver https://gis.stackexchange.com/questions/296170/r-shapefile-transform-longitude-latitude-coordinates↩︎
Una explicación sencilla se encuentra en https://dateandtime.info/es/citycoordinates.php?id=3128760↩︎
https://rdrr.io/cran/sf/man/st_make_grid.html (sintaxis función)↩︎
El concepto “Point in polygon” está ampliamente tratado en el ámbito del análisis espacial. La web gis.stackexchange.com registra miles de entradas sobre este tema.↩︎
https://gis.stackexchange.com/questions/323698/counting-points-in-polygons-with-sf-package-of-r↩︎
https://cran.r-project.org/web/packages/vegan/vignettes/FAQ-vegan.html↩︎
https://cran.r-project.org/web/packages/vegan/vignettes/diversity-vegan.pdf↩︎
Otra posibilidad de cálculo puede consultarse en https://stackoverflow.com/questions/74708843/vegan-package-diversity-calculates-results-incorrectly↩︎
“Trabajando con vegan”; Isabela Cano Osorio y otros. https://rpubs.com/isabelacano↩︎
“Analysis of community ecology data in R”; Jinliang Liu. http://jinliang.weebly.com/uploads/2/5/7/8/25781074/analysis_of_community_ecology_data_in_r.pdf↩︎
Ver por ejemplo el trabajo “Mapping the diversity of street tree inventories across eight cities internationally using open data” https://www.sciencedirect.com/science/article/pii/S1618866721001242↩︎