In ajpelu/respine: Evaluate RESilience of Mediterranean PINE plantations

library("rgbif")
library("rgdal")
library("sp")
library("raster")
library("dplyr")
library("sf")
library("knitr")

Para establecer los valores iniciales de riqueza en cada uno de los tipos de cobertura vegetal vamos a analizar los datos de campo procedentes de los inventarios forestales SINFONEVADA [@PerezLuque2014]. Para ello consultaremos este conjunto de datos en GBIF donde está documentado.

El objetivo es obtener una tabla de riqueza por ecosistema. Los pasos que realizaremos son:

Estudiar el conjunto de datos. Consultar el identificador único (uuid) del dataset: db6cd9d7-7be5-4cd0-8b3c-fb6dd7446472
Descargar el dataset
Obtener los valores de riqueza por ecosistema. Para ello utilizaremos el mapa de ecosistemas generado por el Observatorio de Cambio Global de Sierra Nevada [@Bonet2010]

Obtener el conjunto de datos

### uuid of Sinfonevada
sinfo_uuid <- 'db6cd9d7-7be5-4cd0-8b3c-fb6dd7446472'

### See metadata
sinfo_meta <- datasets(uuid = sinfo_uuid)

# Get table of ocurrences
sf <- occ_data(datasetKey=sinfo_meta$data$key, limit = 8000)

Numero de plots y especies por plot

Obtendremos el numero de plots (agrupando por localización)
Para cada plot calcularemos la riqueza de especies

# Get only the fields of interest  
df <- sf$data %>% dplyr::select(decimalLatitude, decimalLongitude, scientificName)

# How many species by plot
richness_loc <- df %>%
  group_by(decimalLatitude, decimalLongitude) %>% 
  count() %>% 
  as.data.frame() %>% 
  tibble::rownames_to_column(var='id_plot') %>%
  rename(rich = n)

Riqueza por ecosistemas

Analizaremos la riqueza media de los plots por ecosistema
Obtener los ecosistemas de OBSNEV, que están en este enlace. También se han incluido en la carpeta /data del paquete.

# Prepare Ecosystems data 
wd <- getwd()

temporalwd <- setwd(tempdir())
unzip('../inst/extdata/ecosistemas_sn.zip', exdir = temporalwd)

eco <- readOGR(dsn=temporalwd, layer = 'ecosistemas', encoding="UTF-8", verbose = TRUE)

# Transform projection
eco_t <- spTransform(eco, CRS("+init=epsg:4326"))

Ahora creamos una capa vectorial de puntos para los plots del invetnario forestal, y le asignamos su correspondiente proyección.

# Create an spatial point dataframe for the plots 
richness_sp <- SpatialPointsDataFrame(richness_loc[,c("decimalLongitude", "decimalLatitude")],
                                      richness_loc)
projection(richness_sp) <-  CRS("+init=epsg:4326")

Seguidamente realizamos una intersección entre capas. Estamos interesados en obtener la tipología de ecosistema para cada plot del inventario forestal.

# See this example
# https://gis.stackexchange.com/questions/226035/join-spatial-point-data-with-multiple-polygon-data-using-r

# Convert to sf-objects
richness_sp.sf <- st_as_sf(richness_sp)
eco_t.sf <- st_as_sf(eco_t)

# Keep all "meuse.sf", sort by row.names(meuse.sf). Default overlay is "intersects".
aux <- st_join(richness_sp.sf, eco_t.sf[,c('COD_ECOSIS', 'ECOSISTE_1')])

# Convert back to Spatial*
richness_sp_eco <- as(aux, "Spatial")

Vamos a reagrupar los ecosistemas

aux <- aux %>% mutate(newECO = recode_factor(COD_ECOSIS, 
                                                    `8`="Pine plantations",
                                                    `2`="High-mountain meadows",
                                                    `3`="High-mountain shrubland",
                                                    `5`="Mid-mountain shrubland",
                                                    `1`="Pastures",
                                                    `6`="Aquatic systems", 
                                                    `NA`="NA", 
                                                    .default = 'Natural Forests'))

Finalmente computamos los valores de riqueza por ecosistemas y por agregados

richSinfo <- aux %>%
  group_by(ECOSISTE_1, COD_ECOSIS) %>%
  summarise(mean = mean(rich),
            sd = sd(rich),
            se = sd/sqrt(length(rich)),
            n = length(rich),
            min = min(rich),
            max = max(rich),
            median = median(rich)) %>%
    as.data.frame() %>%
  dplyr::select(-geometry)

richSinfo_agg <- aux %>%
  group_by(newECO) %>%
  summarise(mean = mean(rich),
            sd = sd(rich),
            se = sd/sqrt(length(rich)),
            n = length(rich),
            min = min(rich),
            max = max(rich), 
            median = median(rich)) %>%
  as.data.frame() %>%
  dplyr::select(-geometry)

En [@GomezAparicio2009] se analizan el mismo dataset, y se obtienen, tras aplicar los modelos, los siguientes resultados:

richPot <- data.frame(cbind(
  eco = c('Plantations', 'Quercus ilex forests', 'Natural deciduous forests'),
  n = c(442, 45, 26),
  potRich = c(13.09, 14.92, 17.55),
  lowerInterval = c(12.82, 13.72, 15.62), 
  upperInterval = c(13.34, 16.11, 19.66)))