In RS-eco/edc: European Data Cube

knitr::opts_chunk$set(
  collapse = TRUE, comment = "#>", results="asis",
  warning=F, message=F, fig.width=8, fig.height=8
)

# Load packages
library(dplyr); library(magrittr); library(tidyr)
library(sp); library(sf)
library(ggplot2); library(patchwork)

# Load edc package
library(edc)

# Load shapefile of Europe
data("europe", package="edc")

Load EBCC data

# EBCC data for selected bird species
data("ebcc_Sylvia_curruca")
data("ebcc_Turdus_torquatus")
data("ebcc_Phylloscopus_bonelli")
data("ebcc_Vanellus_vanellus")

# For other species, load the entire EBCC data, with:
# data("ebcc_eur")

# Quick look at data
#head(ebcc_Vanellus_vanellus)

# Turn species into one data.frame
ebcc_sub <- bind_rows(list(ebcc_Sylvia_curruca, ebcc_Turdus_torquatus, 
                           ebcc_Phylloscopus_bonelli, ebcc_Vanellus_vanellus)) %>%
  as.data.frame()

# Convert data into a SpatialPointsDataFrame
coordinates(ebcc_sub) <- ~decimalLongitude+decimalLatitude # set coordinates
raster::projection(ebcc_sub) <- CRS("+proj=longlat +datum=WGS84") # set projection

Load Environmental data for Europe

# Corine data for whole of Europe
data("corine_lc_eur_p44deg")

# Select land-cover data for 2018
corine_lc_eur_p44deg <- corine_lc_eur_p44deg %>% select(c(x,y,var, `2018`)) %>% spread(var, `2018`)
#coordinates(clc_eur_grid) <- ~x+y # set coordinates
#projection(clc_eur_grid) <- CRS("+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs") # set 
#gridded(clc_eur_grid) <- TRUE  
#clc_eur_grid <- projectRaster(clc_eur_grid, crs=CRS("+proj=longlat +datum=WGS84"))
#clc_eur_grid
#plot(clc_eur_grid)

# Elevation data 
data("srtm_csi_eur_p44deg")

# Elevation - Rastergrid 
#coordinates(srtm_eur_grid) <- ~x+y # set coordinates
#gridded(srtm_eur_grid) <- TRUE                # grid the data
#projection(srtm_eur_grid) <- CRS("+proj=longlat +datum=WGS84") # set projection
# Turn SpatialPixelsDataFrame into Raster Stack
#srtm_eur_grid <- stack(srtm_eur_grid)
#plot(srtm_eur_grid[[1]])
#srtm_eur_grid

#Climate data
data("cordex_bioclim_eur_p44deg")

# Split data by scenario and timeframe to be able to use rasterFromXYZ
groupkeys <- cordex_bioclim_eur_p44deg %>% 
  tidyr::unite(col="gcm_rcm_rcp", c("gcm", "ensemble", "rcm", "rs", "rcp"), sep="_") %>% 
  group_by(gcm_rcm_rcp, time_frame) %>% group_keys()
cordex_bioclim_eur_p44deg <- cordex_bioclim_eur_p44deg %>% 
  tidyr::unite(col="gcm_rcm_rcp", c("gcm", "ensemble", "rcm", "rs", "rcp"), sep="_") %>% 
  group_split(gcm_rcm_rcp, time_frame, .keep=F)

Merge Environmental data

# Turn data into raster objects
srtm <- srtm_csi_eur_p44deg %>% raster::rasterFromXYZ(crs="+init=EPSG:4326")
corine_lc <- corine_lc_eur_p44deg %>% raster::rasterFromXYZ(crs="+init=EPSG:4326")
bioclim <- cordex_bioclim_eur_p44deg[[1]] %>% raster::rasterFromXYZ(crs="+init=EPSG:4326")

# Merge data
srtm <- raster::extend(srtm, corine_lc)
r_all_dat <- raster::stack(srtm, corine_lc, bioclim)

Extract environmental data for EBCC points

ebcc_env <- raster::extract(r_all_dat, ebcc_sub, sp=T) %>% as.data.frame()

p1 <- ebcc_env %>% ggplot() + geom_violin(aes(x=occurrenceStatus, y=bio1)) + theme_bw() + 
  facet_grid(.~species)
p2 <- ebcc_env %>% ggplot() + geom_violin(aes(x=occurrenceStatus, y=bio12)) + theme_bw() + 
  facet_grid(.~species)
p3 <- ebcc_env %>% ggplot() + geom_violin(aes(x=occurrenceStatus, y=alt_mean)) + 
  theme_bw() + facet_grid(.~species)
p1 / p2 / p3
rm(list=ls()); gc()