3. eucop_data_preparation on Ursus spelaeus
walkthrough
eucop_data_preparation allows users to format raw
occurrences according to the preferred strategy. In the following lines,
we show how eucop_data_preparation automatically imports
the EutherianCop occurrences, defines the backgorund area, and creates
the absence points to model the spatial distribution of the target
species.
a. Load the occurrence of the target species
First, cave bear occurrences are converted into spatial features
object. The spatial extent is reduced to mainland Europe (up to 50°
longitude). Four bioclimatic variables (“bio1,”bio4”,“bio11”,“bio19”)
are selected as predictors. This is accomplished by running
eucop_data_preparation by setting
species_name = Ursus spelaeus,
variables = "bio", and
which.vars = c("bio1,"bio4","bio11","bio19"). Given the
importance of using an equal-area projection for accurate area-based
calculations, eucop_data_preparation internally reprojects
raw data to the Mollweide World projection. The use of the Mollweide
rather than other projection systems is due to its suitability for
global-scale applications.
library(RRgeo)
library(terra)
library(sf)
library(rnaturalearth)
library(ggplot2)
library(ggtext)
library(viridis)
load(system.file("exdata/Ursus_occs.Rda", package="RRgeo"))ne_countries()->globalmap
subset(globalmap,continent=="Europe")->euromap
sf_use_s2(FALSE)
st_crop(euromap, xmin=-20, ymin=2 , xmax=50, ymax=75)->euromap
st_crop(occs, xmin=-20, ymin=2 , xmax=50, ymax=75)->occs
st_transform(euromap, st_crs("ESRI:54009"))->euromap
st_transform(occs, st_crs("ESRI:54009"))->occs
p1<-ggplot()+
geom_sf(data=euromap,col="grey40",fill="white")+
geom_sf(data=occs,fill="darkorchid2",size=2.5,color="black",pch=21)+
labs(title="*Ursus spelaeus* occurrences")+
theme(panel.background = element_rect(fill="aliceblue",colour = "black"),
panel.grid = element_blank(),
axis.text=element_text(size=10),
plot.title = element_markdown(size=12,hjust=0.5),
plot.margin = unit(c(0.1,0.1,0.1,0.1),"cm"))
plot(p1)
b. Minimum Convex Polygon and buffer area creation
Next, eucop_data_preparation calculates the Minimum
Convex Polygon (MCP) which encloses all the fossil localities where
U. spelaeus was present. Then, to account for species’
dispersal abilities and overcome potential bias in the record
preservation, it creates a buffer area around the polygon with a radius
equal to 10% of the maximum distance between the occurrences of U.
spelaeus. The 10% value is set by default and can be changed by
setting the argument in the list.
buff<-0.1 # this is the element "buff" within the argument `bk_points`
sf::st_convex_hull(st_union(occs))->pol
max(st_distance(occs))*buff->buf
st_buffer(pol,dist=as.numeric(buf))->pol_with_buffer
p2<-p1+
geom_sf(data=pol,col="forestgreen",fill="transparent",lwd=1)+
labs(title="*Ursus spelaeus* occurrences and MCP")+
theme(panel.background = element_rect(fill="aliceblue",colour = "black"),
panel.grid = element_blank(),
axis.text= element_text(size=10),
plot.title = element_markdown(size=12,hjust=0.5),
plot.margin = unit(c(0.1,0.1,0.1,0.1),"cm"))
p3<-p1+
geom_sf(data=pol,col="forestgreen",fill="transparent",lwd=1,linetype="dashed")+
geom_sf(data=pol_with_buffer,col="forestgreen",fill="transparent",lwd=1)+
labs(title="*Ursus spelaeus* occurrences and <br> MCP with buffer area")+
theme(panel.background = element_rect(fill="aliceblue",colour = "black"),
panel.grid = element_blank(),
axis.text= element_text(size=10),
plot.title = element_markdown(size=12,hjust=0.5),
plot.margin = unit(c(0.1,0.1,0.1,0.1),"cm"))
plot(p2)
plot(p3)
c. Absence points
Once the background area is defined, all the paleoclimatic/vegetational variables selected by the user are cropped accordingly. Then, the selection of absence points depends on the chosen strategy.
Under the “background” approach (ENFA - presence-only), the values
associated to each grid cell included in the background area represent
the multidimensional space of variables. This is accomplished by setting
bk_points = list(buff = 0.1, bk_strategy = "background", bk_n = 10000)
where all values are defaults. Here, the bk_n argument
defines the maximum number of absence points in each time bin.
Notice: The inclusion of background points is repeated
for each time bin where the target species occur as it is a necessary
requirement for applying ENFA. Consequently, the output files might take
up much disk space. However, this might be a problem only if a very long
list of species_name is provided. In any case, we invite
the users to ensure having sufficient available disk space on their
device before running the analyses.
Unlike ENFA, the “pseudoabsence” approach requires the selection of
specific locations within the accessible area of the species, where
absence points are to be sampled. To end it,
eucop_data_preparation develops on the approach described
in Mondanaro et al. (2021). It starts by generating a density map on all
the occurrences of the target species (irrespective of the time bin), so
that there are more pseudoabsences where the presences are denser. Then,
for each time bin, it samples a number of pseudoabsences which is
proportional to the number of presences in the bin itself, up to a
maximum total number (i.e. summed over all time bins) defined by
bk_n. This is accomplished by setting
bk_points = list(buff = 0.1, bk_strategy = "pseudoabsence", bk_n = 10000).
In the chunk below, we used the map of a random “bio” variable at 35 kyr
as a geographical mask to plot the density map, after converting it by
using the extent and the coordinate raster system (crs) of the
background area.
latesturl<-RRgeo:::get_latest_version("12734585")
curl::curl_download(paste0(latesturl,"/files/X35kya.tif?download=1"),
destfile = "X35kya.tif", quiet = FALSE)
rast("X35kya.tif")->map35
map35[[1]]->mapproject(map,st_crs(pol_with_buffer)$proj4string,res = 50000)->map
mask(crop(map,vect(pol_with_buffer)),vect(pol_with_buffer))->map
RRgeo:::density_background(pres.locs=occs, MASK=map, rm.pres=TRUE)->dens.ras
as.data.frame(dens.ras,xy=TRUE)->dens.probs
colnames(dens.probs)[3]<-"lyr.1"
p_dens<-ggplot()+
geom_sf(data=euromap)+
geom_tile(data=dens.probs,aes(x=x,y=y,fill=lyr.1))+
scale_fill_viridis()+
geom_sf(data=pol_with_buffer,col="forestgreen",fill="transparent",lwd=1)+
labs(title="*Ursus spelaeus* density map",
fill = "Probability")+
theme(panel.background = element_rect(fill="aliceblue",colour = "black"),
panel.grid = element_blank(),
axis.text= element_text(size=10),
axis.title = element_blank(),
plot.title = element_markdown(size=12,hjust=0.5),
plot.margin = unit(c(0.1,0.1,0.1,0.1),"cm"))
plot(p_dens)
Example: we picked the 35 kyr time interval and a random “bio”
variable to produce a toy example of what is described above. This time
interval includes 7 occurrences of U. spelaeus. To calculate
the density map, eucop_data_preparation uses an
RRgeo internal function named
density_background. Additionally, the function excludes the
occupied cells from the density map (the blank cells in the figure) to
prevent misleading conclusions about habitat suitability or the species’
environmental preferences. After partitioning 10000 pseudoabsences over
time, our sample of pseudoabsences for the 35 kyr time interval includes
386 datapoints.
ggplot()+
geom_sf(data=euromap,col="grey40",fill="white")+
geom_sf(data=psabs.points,aes(fill="pseudoabsence"),size=2.5,color="black",pch=21)+
geom_sf(data=pres.points,aes(fill="presence"),size=2.5,color="black",pch=21)+
scale_fill_manual(values=c("darkorchid2","gold2"))+
geom_sf(data=pol_with_buffer,col="forestgreen",fill="transparent",lwd=1)+
labs(title="*Ursus spelaeus* at 35 kya")+
theme(panel.background = element_rect(fill="aliceblue",colour = "black"),
panel.grid = element_blank(),
axis.text=element_text(size=10),
legend.title = element_blank(),
legend.key = element_rect(fill = "white", colour = "white"),
plot.title = element_markdown(size=12,hjust=0.5),
plot.margin = unit(c(0.1,0.1,0.1,0.1),"cm"))