## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(echo = TRUE)

## ----eval=TRUE, message=FALSE, echo=TRUE, warning=FALSE, results = FALSE------
library(USE)
library(terra)
library(sf)
library(ggplot2)

## ----eval=FALSE, message=FALSE------------------------------------------------
# Worldclim <- geodata::worldclim_global(var='bio', res=10, path=getwd())
# envData <- terra::crop(Worldclim, terra::ext(-12, 25, 36, 60))

## ----eval=TRUE, echo=FALSE, message=FALSE, results=FALSE----------------------
envData <- USE::Worldclim_tmp
envData <- terra::rast(envData,  type="xyz")

## ----eval=TRUE, message=FALSE, warning=FALSE, fig.height = 8, fig.width = 8, fig.align='center'----
myCoef <- c(-8, 1.7, -0.1, -0.001) #0.001
names(myCoef) <- c("intercept", "bio1", "quad_bio1", "bio12" )
myVS.po <- USE::SpatialProba(coefs = myCoef,
                             env.rast = envData, 
                             quadr_term = "bio1",
                             marginalPlots = TRUE)
myVS.po$plot

## ----eval=TRUE, message=FALSE, warning=FALSE, fig.height = 8, fig.width = 8, fig.align='center'----
# Convert the probability of occurrence raster to a raster of presence-absenceusing a random binomial draw
new.pres <- terra::app(myVS.po$rast, fun = function(x) {replicate(1, rbinom(n = length(x), 1, x))})

#Sample true occurrences and true absences using a stratified sampling, keeping a sample prevalence fixed to 1
set.seed(123)
myVS.pa <- terra::spatSample(new.pres, 150, "stratified", xy=TRUE)
myVS.pa <- as.data.frame(myVS.pa)
names(myVS.pa) <- c("x", "y", "Observed")

## ----eval=TRUE, message=FALSE, warning=FALSE, fig.height = 8, fig.width = 8, fig.align='center'----
ggplot()+
  tidyterra::geom_spatraster(data = new.pres)+
  tidyterra::scale_fill_whitebox_c("deep", direction=1, na.value = "transparent",
                                   breaks=c(0,1)) +
  geom_point(data=myVS.pa, 
             aes(x=x, y=y, color=as.factor(Observed)),
             alpha=1, size=2, shape= 19)+
     scale_colour_manual(name=NULL,
                        values=c("1"='steelblue',"0"='#A41616'))+
  labs(x="Longitude", 
       y="Latitude", 
       fill="Virtual species",
       colour = "Presence/Absence")+
  theme_light()+
  theme(legend.position = "bottom",  
        legend.background=element_blank(),
        legend.box="vertical",
        panel.grid = element_blank(),
        text = element_text(size=14),  
        legend.text=element_text(size=14), 
        aspect.ratio = 1, 
        panel.spacing.y = unit(2, "lines"))

## ----eval=TRUE----------------------------------------------------------------
myPres <- subset(myVS.pa, myVS.pa$Observed==1)
myPres <- st_as_sf(myPres , coords=c("x", "y"), crs=4326)

## ----eval=TRUE----------------------------------------------------------------
rpc <- rastPCA(envData, stand = TRUE)
dt <- na.omit(as.data.frame(rpc$PCs[[c("PC1", "PC2")]], xy = TRUE))
dt <- sf::st_as_sf(dt, coords = c("PC1", "PC2"))

## ----eval=FALSE, echo=TRUE----------------------------------------------------
# myRes <- USE::optimRes(sdf=dt,
#                     grid.res=c(1:10),
#                     perc.thr = 20,
#                     showOpt = FALSE,
#                     cr=5)

## ----eval=TRUE, echo=FALSE, message=FALSE, results="hide"---------------------
myRes <- list()
myRes$Opt_res <- 5

## ----eval=TRUE----------------------------------------------------------------
myRes$Opt_res

## ----eval=TRUE, message=FALSE-------------------------------------------------
myObs <- USE::uniformSampling(sdf=dt, 
                              grid.res=myRes$Opt_res,
                              n.tr = 5,
                              sub.ts = TRUE,
                              n.ts = 2,
                              plot_proc = FALSE)

## ----eval=TRUE----------------------------------------------------------------
head(myObs$obs.tr)

## ----eval=TRUE, message=FALSE, warning=FALSE, fig.height = 6, fig.width = 6, fig.align='center'----
env_pca <- c(rpc$PCs$PC1, rpc$PCs$PC2)
env_pca <- na.omit(as.data.frame(env_pca))

ggplot(env_pca, aes(x=PC1))+
  geom_density(aes(color="Environment"), linewidth=1 )+
  geom_density(data=data.frame(st_coordinates(myObs$obs.tr)), 
               aes(x=X,  color="Uniform"), linewidth=1)+
  scale_color_manual(name=NULL, 
                     values=c('Environment'='#1E88E5', 'Uniform'='#D81B60'))+     
  labs(y="Density of PC-scores")+
  ylim(0,1)+
  theme_classic()+
  theme(legend.position = "bottom",  
        text = element_text(size=14),  
        legend.text=element_text(size=12))

ggplot(env_pca, aes(x=PC2))+
  geom_density(aes(color="Environment"), linewidth=1 )+
  geom_density(data=data.frame(st_coordinates(myObs$obs.tr)), 
               aes(x=Y,  color="Uniform"), linewidth=1)+
  scale_color_manual(name=NULL, 
                     values=c('Environment'='#1E88E5', 'Uniform'='#D81B60'))+     
  labs(y="Density of PC-scores")+
  ylim(0,1)+
  theme_classic()+
  theme(legend.position = "bottom",  
        text = element_text(size=14),  
        legend.text=element_text(size=12))

## ----eval=TRUE, message=FALSE-------------------------------------------------
myGrid.psAbs <- USE::paSampling(env.rast=envData,
                                pres=myPres,
                                thres=0.75,
                                H=NULL,
                                grid.res=as.numeric(myRes$Opt_res),
                                n.tr = 5,
                                prev=1,
                                sub.ts=TRUE,
                                n.ts=5,
                                plot_proc=FALSE,
                                verbose=FALSE)

## ----eval=TRUE, message=FALSE, warning=FALSE, fig.height = 6, fig.width = 6, fig.align='center'----
ggplot(env_pca, aes(x=PC1))+
  geom_density(aes(color="Environment"), linewidth=1 )+
  geom_density(data=data.frame(st_coordinates(myGrid.psAbs$obs.tr)), 
               aes(x=X,  color="Uniform"), linewidth=1)+
  geom_density(data=terra::extract(c(rpc$PCs$PC1, rpc$PCs$PC2), myPres, df=TRUE), 
               aes(x=PC1, color="Presence"), linewidth=1 )+
  scale_color_manual(name=NULL, 
                     values=c('Environment'='#1E88E5', 'Uniform'='#D81B60', "Presence"="black"))+     
  labs(y="Density of PC-scores")+
  ylim(0,1)+
  theme_classic()+
  theme(legend.position = "bottom",  
        text = element_text(size=14),  
        legend.text=element_text(size=12))

ggplot(env_pca, aes(x=PC2))+
  geom_density(aes(color="Environment"), linewidth=1 )+
  geom_density(data=data.frame(st_coordinates(myGrid.psAbs$obs.tr)), 
               aes(x=Y,  color="Uniform"), linewidth=1)+
  geom_density(data=terra::extract(c(rpc$PCs$PC1, rpc$PCs$PC2), myPres, df=TRUE), 
               aes(x=PC2, color="Presence"), linewidth=1 )+
  scale_color_manual(name=NULL, 
                     values=c('Environment'='#1E88E5', 'Uniform'='#D81B60', "Presence"="black"))+     
  labs(y="Density of PC-scores")+
  ylim(0,1)+
  theme_classic()+
  theme(legend.position = "bottom",  
        text = element_text(size=14),  
        legend.text=element_text(size=12))

## ----eval=TRUE, message=FALSE, warning=FALSE, fig.height = 8, fig.width = 8, fig.align='center'----
ggplot()+
  tidyterra::geom_spatraster(data = new.pres)+
  tidyterra::scale_fill_whitebox_c("deep", direction=1, na.value = "transparent",
                                   breaks=c(0,1)) +
  geom_sf(data=myPres, 
          aes(color= "Presences"), 
          alpha=1, size=2, shape= 19)+
  geom_sf(data=st_as_sf(st_drop_geometry(myGrid.psAbs$obs.tr), 
                        coords = c("x","y"), crs=4326), 
          aes(color="Pseudo-absences"), 
          alpha=0.8, size=2, shape = 19 )+
    scale_colour_manual(name=NULL,
                        values=c('Presences'='steelblue','Pseudo-absences'='#A41616'))+
  labs(x="Longitude", 
       y="Latitude", 
       fill="Virtual species")+
  theme_light()+
  theme(legend.position = "bottom",  
        legend.background=element_blank(),
        legend.box="vertical",
        panel.grid = element_blank(),
        text = element_text(size=14),  
        legend.text=element_text(size=14), 
        aspect.ratio = 1, 
        panel.spacing.y = unit(2, "lines"))

## ----eval=TRUE, message=FALSE, warning=FALSE, fig.height = 6, fig.width = 9, fig.align='center'----
USE::thresh.inspect(env.rast=envData,
                    pres=myPres,
                    thres=c(0.1, 0.25, 0.5, 0.75, 0.9),
                    H=NULL
                    )