Species richness, change, and beta diversity

Author

Florencia Grattarola

Published

September 20, 2024

Calculate species richness, beta diversity and spatial and temporal dissimilarity

library(knitr)
library(terra)
library(vegan)
library(tmap)
tmap_mode('view')
library(sf)
sf::sf_use_s2(FALSE)
library(tidyverse)

Base maps

equalareaCRS <-  '+proj=laea +lon_0=-73.125 +lat_0=0 +datum=WGS84 +units=m +no_defs'
latam <- st_read('data/latam.gpkg', layer = 'latam', quiet = T)
countries <- st_read('data/latam.gpkg', layer = 'countries', quiet = T)
latam_land <- st_read('data/latam.gpkg', layer = 'latam_land', quiet = T)
latam_raster <- rast('data/latam_raster.tif', lyrs='latam')

Data

hyagouaroundi_P.pred <- readRDS('data/Ppred_samples/hyagouaroundi_P.pred_sample.rds')
cbrachyurus_P.pred <- readRDS('data/Ppred_samples/cbrachyurus_P.pred_sample.rds')
ebarbara_P.pred <- readRDS('data/Ppred_samples/ebarbara_P.pred_sample.rds')
lwiedii_P.pred <- readRDS('data/Ppred_samples/lwiedii_P.pred_sample.rds')
pbrasiliensis_P.pred <- readRDS('data/Ppred_samples/pbrasiliensis_P.pred_sample.rds')

Species richness

Code

time1 <- data.frame(hyagouaroundi=robustbase::colMedians(hyagouaroundi_P.pred[,1:2180]), 
                    cbrachyurus=robustbase::colMedians(cbrachyurus_P.pred[,1:2180]),
                    ebarbara=robustbase::colMedians(ebarbara_P.pred[,1:2180]),
                    lwiedii=robustbase::colMedians(lwiedii_P.pred[,1:2180]), 
                    pbrasiliensis=robustbase::colMedians(pbrasiliensis_P.pred[,1:2180]))

time2 <- data.frame(hyagouaroundi=robustbase::colMedians(hyagouaroundi_P.pred[,(2180+1):4360]), 
                    cbrachyurus=robustbase::colMedians(cbrachyurus_P.pred[,(2180+1):4360]),
                    ebarbara=robustbase::colMedians(ebarbara_P.pred[,(2180+1):4360]),
                    lwiedii=robustbase::colMedians(lwiedii_P.pred[,(2180+1):4360]), 
                    pbrasiliensis=robustbase::colMedians(pbrasiliensis_P.pred[,(2180+1):4360]))

# sum of occurrence probabilities per grid cell
richness_time1 <- time1 %>% rowSums()
richness_time2 <- time2 %>% rowSums()

# values are the same for all species (see check below) thus I can use one species as example
cbrachyurus_expert_gridcell <- readRDS('data/species_POPA_data/cbrachyurus_expert_gridcell.rds')
PO_time1 <- readRDS('data/species_POPA_data/data_cbrachyurus_PO_time1.rds') %>%
  cbind(expert=cbrachyurus_expert_gridcell$dist_exprt) %>% 
  filter(!is.na(env.bio_10) &!is.na(env.bio_14) & !is.na(env.elev) & !is.na(env.grass) & 
           !is.na(acce) & !is.na(count) & !is.na(expert)) 

PO_time2  <- readRDS('data/species_POPA_data/data_cbrachyurus_PO_time2.rds')  %>%
  cbind(expert=cbrachyurus_expert_gridcell$dist_exprt) %>% 
  filter(!is.na(env.bio_10) &!is.na(env.bio_14) & !is.na(env.elev) & !is.na(env.grass) & 
           !is.na(acce) & !is.na(count) & !is.na(expert)) 
PO_time1_time2 <- rbind(PO_time1 %>% mutate(time=1), PO_time2 %>% mutate(time=2)) 

preds.richness <- data.frame(PO_time1_time2[,c('X', 'Y', 'pixel', 'area')], 
                             richness_time1, richness_time2)

# rasterise
rast <- latam_raster
rast[] <- NA
rast1 <- rast2 <- rast3 <- terra::rast(rast)

rast1[preds.richness$pixel] <- preds.richness$richness_time1
rast2[preds.richness$pixel] <- preds.richness$richness_time2

rast1 <- rast1 %>% terra::mask(., terra::vect(latam_land))
rast2 <- rast2 %>% terra::mask(., terra::vect(latam_land))

names(rast1) <- 'time1'
names(rast2) <- 'time2'

Plots

Code

time1MAP <- 
  # tm_graticules(alpha = 0.3) +  
  tm_shape(rast1) +
  tm_raster(palette = 'Oranges', midpoint = NA, 
            # alpha = 0.7,
            style= 'cont', breaks = seq(0, 5, 1), 
            title='SR time1') + 
  # tm_shape(countries) +
  # tm_borders(col='grey60', alpha = 0.4) + 
  # tm_layout(legend.outside = T, frame.lwd = 0.3, scale=1.2, legend.outside.size = 0.1) +
  tm_view(set.view = 3) 

time1MAP

Code

# tmap::tmap_save(tm = time1MAP, filename = 'docs/figs/time1MAP_SR.svg', device = svglite::svglite)

# Map of the of the probability of occurrence in the second period (time2: 2014-2021)
time2MAP <- 
  # tm_graticules(alpha = 0.3) +  
  tm_shape(rast2) +
  tm_raster(palette = 'Purples', midpoint = NA, 
            alpha = 0.7,
            style= 'cont', breaks = seq(0, 5, 1),
            title='SR time2')  + 
  # tm_shape(countries) +
  # tm_borders(col='grey60', alpha = 0.4) + 
  # tm_layout(legend.outside = T, frame.lwd = 0.3, scale=1.2, legend.outside.size = 0.1) +
  tm_view(set.view = 3) 

time2MAP

Code

# tmap::tmap_save(tm = time2MAP, filename = 'docs/figs/time2MAP_SR.svg', device = svglite::svglite)

diffMAP <- 
  # tm_graticules(alpha = 0.3) + 
  tm_shape(rast2 - rast1) +
  tm_raster(palette = 'PiYG', midpoint = 0, 
            #alpha = 0.7,
            style = 'cont', breaks = seq(-4, 2, 1),
            title='SR change') +
  # tm_shape(countries) +
  # tm_borders(col='grey60', alpha = 0.4) + 
  # tm_layout(legend.outside = T, frame.lwd = 0.3, scale=1.2, legend.outside.size = 0.1) +
  tm_view(set.view = 3) 

diffMAP

Code

# tmap::tmap_save(tm = diffMAP, filename = 'docs/figs/diffMAP_SR.svg', device = svglite::svglite)

SR_change <- as_tibble(rast2 - rast1) %>% 
  rename(diffSR=time2) %>% 
  mutate(change=ifelse(diffSR<(-0.112), 'neg', 
                       ifelse(diffSR>0.112, 'pos', 'zero')))

SR_change_hist <- ggplot(SR_change) +
  geom_histogram(aes(x=diffSR, fill=change), bins = 20) +
  scale_fill_manual(values=c('#d01c8b', '#4dac26', '#f7f7f7')) +
  labs(x='Change in SR', y='N of grid cells') + 
  ggpubr::theme_pubclean()

SR_change_hist

IUCN species richness map

Code

cthous_IUCN <- sf::st_read('big_data/cthous_IUCN.shp', quiet = T) %>% sf::st_transform(crs=equalareaCRS)
hyagouaroundi_IUCN <- sf::st_read('big_data/hyagouaroundi_IUCN.shp', quiet = T) %>% sf::st_transform(crs=equalareaCRS)
cbrachyurus_IUCN <- sf::st_read('big_data/cbrachyurus_IUCN.shp', quiet = T) %>% sf::st_transform(crs=equalareaCRS)
ebarbara_IUCN <- sf::st_read('big_data/ebarbara_IUCN.shp', quiet = T) %>% sf::st_transform(crs=equalareaCRS)
lwiedii_IUCN <- sf::st_read('big_data/lwiedii_IUCN.shp', quiet = T) %>% sf::st_transform(crs=equalareaCRS)
pbrasiliensis_IUCN <- sf::st_read('big_data/pbrasiliensis_IUCN.shp', quiet = T) %>% sf::st_transform(crs=equalareaCRS)

# Latam_grid <- st_make_grid(Latam, cellsize = 100000) %>% 
#   st_intersection(Latam)  %>% 
#   st_sf(gridID=1:length(.), geometry= .) %>% 
#   st_make_valid() %>% st_cast() 

carnivores <- do.call(rbind, list(hyagouaroundi_IUCN %>% select(binomial),
                                  cbrachyurus_IUCN %>% select(binomial), 
                                  ebarbara_IUCN %>% select(binomial), 
                                  lwiedii_IUCN %>% select(binomial), 
                                  pbrasiliensis_IUCN %>% select(binomial))) %>% 
  st_cast()

# SR_carnivores <- st_join(Latam_grid, carnivores) %>% 
#   group_by(gridID) %>% 
#   summarise(SR=ifelse(n_distinct(binomial, na.rm=T)!=0, n(), 0)) %>% 
#   st_cast() %>% sf::st_transform(crs=4326)  

SR_carnivores <- vect(carnivores) %>% 
  rasterize(., latam_raster, fun='count', 'binomial') 
names(SR_carnivores) <- 'SR'

Plot

Code

tmap_mode('view')
tmap_mode('plot')

IUCN_SR <- 
  #tm_graticules(alpha = 0.3) +  
  tm_shape(SR_carnivores) +
  tm_raster(palette = 'Greys',
            alpha = 0.7,
            style = 'cont', 
            breaks = seq(0, 5, 1),
            title='SR IUCN range maps') +
  # tm_shape(countries) +
  # tm_borders(col='grey60', alpha = 0.4) + 
  # tm_layout(legend.outside = T, 
  #           frame.lwd = 0.3, scale=1.2, legend.outside.size = 0.1) +
  tm_view(set.view = 3) 

IUCN_SR

Code

# tmap::tmap_save(tm = IUCN_SR, filename = 'docs/figs/IUCN_SR.svg', device = svglite::svglite)

Beta diversity

Code

# ratio gamma / alpha -1
gamma <- ncol(time1)
alpha1 <- mean(richness_time1)
beta1 <- round(gamma/alpha1, 3)

# confidence interval
sd1 <- sd(richness_time1)
beta1.high <- round(gamma/(alpha1+sd1), 3)
beta1.low <- round(gamma/(alpha1-sd1), 3) 

tibble(gamma=gamma, alpha1=alpha1, 
       sd1=sd1, beta1=beta1, 
       beta1.high=beta1.high, 
       beta1.low=beta1.low) %>% 
  kable()

gamma	alpha1	sd1	beta1	beta1.high	beta1.low
5	2.61632	1.505596	1.911	1.213	4.502

Code

# ratio gamma / alpha -1
gamma <- ncol(time2)
alpha2 <- mean(richness_time2)
beta2 <- round(gamma/alpha2, 3)

# confidence interval
sd2 <- sd(richness_time2)
beta2.high <- round(gamma/(alpha2+sd2), 3)
beta2.low <- round(gamma/(alpha2-sd2), 3) 

tibble(gamma=gamma, alpha2=alpha2, 
       sd2=sd2, beta2=beta2, 
       beta2.high=beta2.high, 
       beta2.low=beta2.low) %>% 
  kable()

gamma	alpha2	sd2	beta2	beta2.high	beta2.low
5	2.395093	1.467106	2.088	1.295	5.388

Plot

Code

beta_w_time1 <- gamma/richness_time1
beta_w_time2 <- gamma/richness_time2

betas_time <- tibble(beta = c(beta_w_time1, beta_w_time2),
                    time = c(rep('time1', length(beta_w_time1)),
                             rep('time2', length(beta_w_time2))))
## betas
betas_plot <- ggplot(betas_time, aes(x=time, y=beta, fill=time)) +
  geom_boxplot(show.legend = F) + 
  scale_fill_manual(values = c('#e66101', '#5e3c99')) + 
  scale_y_continuous(limits = quantile(richness_time1, c(0.1, 0.9)))  +
  labs(x='', y=expression(beta)) + 
  ggpubr::theme_pubclean()

betas_plot

Code

## alphas
alphas <- data.frame(site=1:nrow(time1), 
                     alpha=c(richness_time1, richness_time2),
                     time=c(rep('time1', length(alpha1)),
                            rep('time2', length(alpha2))))

alphas_plot <- ggplot(alphas, aes(x = time, y=alpha, fill=time)) + 
  geom_boxplot(show.legend = F) +
  scale_fill_manual(values = c('#e66101', '#5e3c99')) + 
  labs(x='', y=expression(alpha)) + 
  ggpubr::theme_pubclean()

alphas_plot

Code

# ggsave(alphas_plot, filename='docs/figs/alphas_plot.svg', 
#        width = 3, height = 5, device = 'svg', dpi=200)

Spatial dissimilarity

Code

b.jac1 <- vegdist(time1, method = 'jaccard')
b.jac2 <- vegdist(time2, method = 'jaccard')

# dist
XY <- preds.richness[1:2180,1:2]
dist_XY <- dist(XY, method = 'euclidian')
dist_XY <- as.matrix(dist_XY)

XY_sf <- st_as_sf(XY, coords=c('X','Y'), crs=equalareaCRS)

gridsToKeep <- XY_sf %>%
  pull(geometry) %>%
  st_distance() %>%
  as.data.frame.table() %>%
  filter(Freq < units::set_units(500000, 'm')) %>% 
  mutate(index_keep = paste0(Var1,'-',Var2)) %>% 
  select(index_keep)

####
# select grid closer than 200km 
dist_XY_200km <- dist_XY %>%
  as_tibble() %>% 
  mutate(row = row_number()) %>% 
  pivot_longer(-row, names_to = 'column',values_to = 'values', names_prefix = 'V') %>% 
  mutate(index = paste0(row,'-',column)) %>% 
  filter(index %in% gridsToKeep$index_keep) %>% 
  select(dist=values)

# distance decay
b.jac1_200km <- b.jac1 %>% as.matrix() %>% 
  as_tibble() %>% 
  mutate(row = row_number()) %>% 
  pivot_longer(-row, names_to = 'column',values_to = 'values', names_prefix = 'V') %>% 
  mutate(index = paste0(row,'-',column)) %>% 
  filter(index %in% gridsToKeep$index_keep) %>% 
  select(jac_time1=values)

b.jac2_200km <- b.jac2 %>% as.matrix() %>% 
  as_tibble() %>% 
  mutate(row = row_number()) %>% 
  pivot_longer(-row, names_to = 'column',values_to = 'values', names_prefix = 'V') %>% 
  mutate(index = paste0(row,'-',column)) %>% 
  filter(index %in% gridsToKeep$index_keep) %>% 
  select(jac_time2=values) 

spatial_beta <- bind_cols(dist_XY_200km, b.jac1_200km, b.jac2_200km)

Plot

Code

beta_boxplot <- spatial_beta %>% 
  rename(time1=jac_time1, time2=jac_time2) %>% 
  mutate(dist=as.factor(round(dist/1000,0))) %>% 
  filter(dist!=0) %>% 
  pivot_longer(-dist, names_to = 'time',values_to = 'dissimilarity') %>% 
  ggplot() +
  geom_boxplot(aes(y = dissimilarity, x=dist, fill=time), alpha=0.8, outlier.shape = NA) +
  scale_fill_manual(values = c('#e66101', '#5e3c99')) + 
  stat_summary(fun = median, geom = 'path', linewidth=1,
               mapping = aes(y = dissimilarity, x=dist, group=time, col=time), 
               show.legend = F) +
  scale_color_manual(values = c('#e66101', '#5e3c99')) + 
  labs(x='Distance (km)', y='Beta diversity (Ružička index)', col='', fill='') +
  theme_bw() 

beta_boxplot

Code

# ggsave(beta_boxplot, filename='docs/figs/beta_boxplot.svg', 
#        width = 7, height = 5, device = 'svg', dpi=300)

Temporal dissimilarity

Code

calculateBetaTemporal <- function(comm_time1, comm_time2, index='jaccard'){
  beta.temp <- data.frame(row=numeric(), beta.temp=numeric())
  for(i in 1:nrow(comm_time1)){
    comm.matrix <- bind_rows(comm_time1[i,], comm_time2[i,])
    beta.temp_i <- data.frame(row=i, beta.temp=as.numeric(vegdist(comm.matrix, index)))
    beta.temp <- rbind(beta.temp, beta.temp_i)
  }
  return(beta.temp)
}

beta.temp <- calculateBetaTemporal(time1, time2)

preds.richness <- bind_cols(preds.richness[1:2180,], beta.temp=beta.temp[,2]) %>% 
  select(pixel, beta.temp)

rast <- latam_raster
rast[] <- NA
rast3 <- terra::rast(rast)

rast3[preds.richness$pixel] <- preds.richness$beta.temp

rast3 <- rast3 %>% terra::mask(., terra::vect(latam_land))

names(rast3) <- 'beta.temp'

Plot

Code

beta_map <- 
  # tm_graticules(alpha = 0.3) +  
  tm_shape(rast3) +
  tm_raster(palette = '-RdPu', midpoint = NA, 
            alpha=0.7,
            style= 'cont')  + 
  # tm_shape(countries) +
  # tm_borders(col='grey60', alpha = 0.4) + 
  # tm_layout(legend.outside = T, frame.lwd = 0.3, scale=1.2)  +
  tm_view(set.view = 3) 

beta_map

Code

# tmap::tmap_save(tm = beta_map, filename = 'docs/figs/beta_map.svg', device = svglite::svglite)