In RS-eco/bdc: Bavarian Data Cube

knitr::opts_chunk$set(
  collapse = T, comment = NA, warning=F, message=F, eval=F,
  echo=T, error=F, comment = "#>", fig.path="../figures/", fig.width=7, fig.height=7
)

Data Analysis

Here, I perform a first species richness calculation of the ASK observations for Bavaria.

Load data

First, we load and subset the data:

rm(list=ls()); invisible(gc())

#+ data_setup, results="hide"
# Load packages
library(data.table); library(dtplyr); library(dplyr, warn.conflicts = FALSE)
library(magrittr); library(tidyr)

########################
# Set file directory
filedir <- sub("/vignettes", "", getwd())

# Load ASK database
my_db <- DBI::dbConnect(RSQLite::SQLite(), dbname = paste0(filedir, "/rawdata/ASK.db"))

# Pull part of ASK database including data on species 
ask_data <- dplyr::tbl(my_db, paste("ask_art")) %>% collect()

# Load background grid
tk25 <- dplyr::tbl(my_db, "geo_tk25_quadranten") %>% dplyr::collect() %>%
  tidyr::unite(col="quadrant", c("KARTE", "QUADRANT"), sep="/", remove = TRUE) %>% 
  dplyr::select(-c("KARTE_QUAD")) 

# Load taxonomy
load(paste0(filedir, "/data/taxonomyStd.rda"))

# Load ask_fuo
ask_fuo <- dplyr::tbl(my_db, "ask_fuo") %>% collect()

# Load Red List data
fuf_art <- dplyr::tbl(my_db, "fuf_art") %>% collect()
DBI::dbDisconnect(my_db)

# combine gridded map from dat_ask with full locations of recorded species and taxonomy
dat_ask <- ask_data %>% dplyr::select(-c(gkk_rw, gkk_hw, objnr, rw, hw)) %>% 
  dplyr::left_join(ask_fuo) %>%  dplyr::left_join(fuf_art %>% dplyr::select(-art), by="art_id") %>% 
  dplyr::left_join(tk25) %>% dplyr::left_join(taxonomyStd, by=c("art" = "scientificName"))
dat_ask <- dat_ask %>% tidyr::drop_na(class, order)
dat_ask$KARTE_QUAD <- as.numeric(sub("/", "", dat_ask$quadrant))

#nrow(dat_ask)
ask_tk25 <- dat_ask %>% dplyr::select(KARTE_QUAD) %>% drop_na()
#nrow(ask_tk25)

ask_gkk <- ask_data %>% dplyr::select(gkk_rw, gkk_hw) %>% drop_na()
#nrow(ask_gkk)

# Create custom taxon vector (Aves, Lepidoptera, Odonata, Orthoptera)
dat_ask$class_order <- "Aves"
dat_ask$class_order[dat_ask$class == "Insecta"] <- dat_ask$order[dat_ask$class == "Insecta"]
unique(dat_ask$class_order)
unique(dat_ask$family)

dat_lepi <- dat_ask %>% filter(class_order == "Lepidoptera") %>% 
  filter(family == "Papilionidae" | family == "Hesperiidae" | family == "Pieridae" | 
           family == "Nymphalidae" | family == "Lycaenidae" | family == "Riodinidae")
dat_non_lepi <-  dat_ask %>% filter(class_order != "Lepidoptera")
dat_ask <- bind_rows(dat_lepi, dat_non_lepi)

dat_ask %<>% dplyr::select(XQMITTE, YQMITTE, KARTE_QUAD, jahr, mon, rld,
                           class, class_order, scientificNameStd, art)
rm(ask_fuo, ask_data, my_db); invisible(gc())

# ## Subset data by year
dat_ask <- dat_ask %>% filter(class == "Insecta") %>%
  filter(jahr >= 1985, jahr <= 2019) %>% dplyr::select(-class); invisible(gc())

# Summarise number of species per taxon
dat_all <- dat_ask %>% group_by(class_order) %>%
  summarise(n_total=n_distinct(scientificNameStd))

# Get species overview for Red List Species
dat_rld <- dat_ask %>% filter(rld %in% c(0,1,2,3)) %>%
  group_by(class_order) %>%
  summarise(n_rld=n_distinct(scientificNameStd))

# Summary of number of species
left_join(dat_rld, dat_all) %>%
  mutate(perc_rld = n_rld/n_total*100)

Calculate richness

Next we split the data into the different taxonomic groups, calculate the species richness and turn them into a raster stack:

# Calculate species richness for each taxonomic group
sum_dat <- dat_ask %>% group_by(XQMITTE, YQMITTE, class_order) %>% 
  summarise(sr_std=n_distinct(scientificNameStd)) %>% drop_na()

list_dat <- sum_dat %>% group_by(class_order) %>% group_split(.keep=F)
r_names <- sum_dat %>% group_by(class_order) %>% group_keys() %>% unlist()
rast_dat <- lapply(list_dat, function(x){
  dat <- raster::rasterFromXYZ(x, crs=sp::CRS("+init=epsg:31468"), res=c(6100, 5550), digits=0)
  raster::extend(dat, y=raster::extent(c(4279033, 4638933, 5236691, 5608687)))
})
rast_dat <- raster::stack(rast_dat)
names(rast_dat) <- r_names
rast_dat

# Calculate total species richness
sum_dat <- dat_ask %>% group_by(XQMITTE, YQMITTE) %>% 
  summarise(sr_std=n_distinct(scientificNameStd)) %>% drop_na()
rast_dat2 <- raster::rasterFromXYZ(sum_dat, crs=sp::CRS("+init=epsg:31468"), 
                                   res=c(6100, 5550), digits=0)
rast_dat2 <- raster::extend(rast_dat2, 
                            y=raster::extent(c(4279033, 4638933, 5236691, 5608687)))

rast_dat <- raster::stack(list(rast_dat, rast_dat2))
names(rast_dat)[raster::nlayers(rast_dat)] <- "Total"

Plot species richness by group:

raster::plot(rast_dat)

Save output as .tif files:

raster::writeRaster(rast_dat, 
            filename=paste0(filedir, "/rawdata/sr_ask_1985_2019_epsg31468_tk4tel.tif"),
            format="GTiff", options=c("COMPRESS=NONE", "TFW=YES"), datatype='INT2U', 
            overwrite=T, bylayer=T, suffix="names", prj=T, setStatistics=T)

Calculate Red List Germany richness

Next we split the data into the different taxonomic groups, calculate the species richness and turn them into a raster stack:

# Calculate species richness for each taxonomic group
sum_dat <- dat_ask %>% filter(rld %in% c(0,1,2,3)) %>% 
  group_by(XQMITTE, YQMITTE, class_order) %>% 
  summarise(sr_std=n_distinct(scientificNameStd)) %>% drop_na()

xy_all <- dat_ask %>% group_by(XQMITTE, YQMITTE, class_order) %>% group_keys() %>%
  mutate(xy_all=NA)
sum_dat <- full_join(sum_dat, xy_all) %>% dplyr::select(-c(xy_all))

list_dat <- sum_dat %>% group_by(class_order) %>% group_split(.keep=F)
r_names <- sum_dat %>% group_by(class_order) %>% group_keys() %>% unlist()
rast_dat <- lapply(list_dat, function(x){
  dat <- raster::rasterFromXYZ(x, crs=sp::CRS("+init=epsg:31468"), res=c(6100, 5550), digits=0)
  dat <- raster::extend(dat, y=raster::extent(c(4279033, 4638933, 5236691, 5608687)))
  return(dat)
})
rast_dat <- raster::stack(rast_dat)
names(rast_dat) <- r_names
rast_dat

# Calculate total species richness
sum_dat <- dat_ask %>% filter(rld %in% c(0,1,2,3)) %>% group_by(XQMITTE, YQMITTE) %>% 
  summarise(sr_std=n_distinct(scientificNameStd)) %>% drop_na()
xy_all <- dat_ask %>% group_by(XQMITTE, YQMITTE) %>% group_keys() %>%
  mutate(xy_all=NA)
sum_dat <- full_join(sum_dat, xy_all) %>% dplyr::select(-c(xy_all))
rast_dat2 <- raster::rasterFromXYZ(sum_dat, crs=sp::CRS("+init=epsg:31468"), 
                                   res=c(6100, 5550), digits=0)
rast_dat2 <- raster::extend(rast_dat2, y=raster::extent(c(4279033, 4638933, 5236691, 5608687)))

rast_dat <- raster::stack(list(rast_dat, rast_dat2))
names(rast_dat)[raster::nlayers(rast_dat)] <- "Total"

Plot species richness by group:

raster::plot(rast_dat)

Save output as .tif files:

raster::writeRaster(rast_dat, 
                    filename=paste0(filedir, 
                                    "/rawdata/sr_ask_rld_1985_2019_epsg31468_tk4tel.tif"),
                    format="GTiff", options=c("COMPRESS=NONE", "TFW=YES"), datatype='INT2U', 
                    overwrite=T, bylayer=T, suffix="names", prj=T, setStatistics=T)

Apply beals() smoothing

Please consider this study:

De Cáceres, M. & Legendre, P. 2008. Beals smoothing revisited. Oecologia 156: 657--669.

before using beals() smoothing.

# Bring data into correct format
sum_dat <- dat_ask %>% group_by(XQMITTE, YQMITTE, scientificNameStd) %>% 
  summarise(no_obs=n()) %>% drop_na() %>% ungroup() %>%
  pivot_wider(names_from=scientificNameStd, values_from=no_obs, values_fill=0)
rm(dat_ask); gc()

# Run beals smoothing for first 10 species, need to remove XQMITTE & YQMITTE 
# (thus selection starts with column 3)

## Remove target species: Yields lower values.
library(vegan)
beals_dat1 <- beals(sum_dat[,3:12], 
                    include = FALSE)

## Keep target species (default)
beals_dat2 <- beals(sum_dat[,3:12], 
                    include=TRUE)

## Remove the bias of tarbet species: Yields lower values.
beals_dat3 <- beals(sum_dat[,3:12], 
                    type=3, include=FALSE)
## Uses abundance information?

Plot smoothed values against presence or absence of species

pa <- decostand(sum_dat[,3:12], "pa")
par(mfrow=c(1,3))
boxplot(as.vector(beals_dat1) ~ unlist(pa), xlab="Presence", ylab="Beals (1)", notch=T)
boxplot(as.vector(beals_dat2) ~ unlist(pa), xlab="Presence", ylab="Beals (2)", notch=T)
boxplot(as.vector(beals_dat3) ~ unlist(pa), xlab="Presence", ylab="Beals (3)", notch=T)
rm(list=ls()); invisible(gc())