NLMR_virtualspecies_ibis-isdm.Rmd

---
title: "MAMBO sim framework tests"
author: "Simon Rolph"
date: "`r Sys.Date()`"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

Requirements

```{r}

#devtools::install_github("ropensci/NLMR")
#install.packages("virtualspecies")

#remotes::install_github("IIASA/ibis.iSDM") #0.07


#install.packages("INLA",repos=c(getOption("repos"),INLA="https://inla.r-inla-download.org/R/stable"), dep=TRUE)
#install.packages("inlabru")


library(NLMR)
library(virtualspecies)

library(dplyr)
library(tidyr)


library(ibis.iSDM)
library(INLA)
library(inlabru)
library(xgboost)
library(terra)
library(uuid)
library(assertthat)

library(sf)

sessionInfo()
```


```
R version 4.2.2 (2022-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 17763)

Matrix products: default

locale:
[1] LC_COLLATE=English_United Kingdom.1252  LC_CTYPE=English_United Kingdom.1252   
[3] LC_MONETARY=English_United Kingdom.1252 LC_NUMERIC=C                           
[5] LC_TIME=English_United Kingdom.1252    

attached base packages:
[1] parallel  stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] sf_1.0-9             assertthat_0.2.1     uuid_1.1-0           terra_1.7-39        
 [5] xgboost_1.7.5.1      inlabru_2.8.0        INLA_23.04.24        foreach_1.5.2       
 [9] Matrix_1.5-1         ibis.iSDM_0.0.7      tidyr_1.3.0          dplyr_1.1.2         
[13] virtualspecies_1.5.1 raster_3.6-14        sp_1.5-1             NLMR_1.1.1          

loaded via a namespace (and not attached):
 [1] tidyselect_1.2.0   xfun_0.37          purrr_1.0.1        splines_4.2.2     
 [5] lattice_0.20-45    colorspace_2.1-0   vctrs_0.6.2        generics_0.1.3    
 [9] htmltools_0.5.4    yaml_2.3.7         utf8_1.2.3         rlang_1.1.1       
[13] e1071_1.7-13       pillar_1.9.0       glue_1.6.2         withr_2.5.0       
[17] DBI_1.1.3          lifecycle_1.0.3    munsell_0.5.0      gtable_0.3.3      
[21] evaluate_0.20      codetools_0.2-18   knitr_1.42         fastmap_1.1.0     
[25] class_7.3-20       fansi_1.0.4        proto_1.0.0        Rcpp_1.0.10       
[29] KernSmooth_2.23-20 backports_1.4.1    scales_1.2.1       classInt_0.4-8    
[33] checkmate_2.1.0    jsonlite_1.8.4     farver_2.1.1       digest_0.6.31     
[37] ggplot2_3.4.2      grid_4.2.2         cli_3.6.0          tools_4.2.2       
[41] magrittr_2.0.3     proxy_0.4-27       tibble_3.2.1       pkgconfig_2.0.3   
[45] data.table_1.14.8  rmarkdown_2.21     rstudioapi_0.14    iterators_1.0.14  
[49] R6_2.5.1           units_0.8-1        compiler_4.2.2  
```

Generate a random landscape

```{r}
#create 3 variables
env_variables <- terra::rast(
  list(
    env1 = terra::rast(NLMR::nlm_mpd(ncol = 1000,nrow = 1000,roughness = 0.4)),
    env2 = terra::rast(NLMR::nlm_mpd(ncol = 1000,nrow = 1000,roughness = 0.5)),
    env3 = terra::rast(NLMR::nlm_mpd(ncol = 1000,nrow = 1000,roughness = 0.61))
    )
  )

# resample
# target raster
s <- rast(nrows=nrow(env_variables), ncols=ncol(env_variables), xmin=-1.5, xmax=0, ymin=53, ymax=54,nlyrs=3,names = names(env_variables),crs = "EPSG:4326")

#put the values in the new raster
values(s) <- values(env_variables)
env_variables <- s

#visualise
plot(env_variables)

```

## generate some virtual species

```{r}
#generate species using virtual species package and the simualted landscapes
virtual_species1 <- virtualspecies::generateRandomSp(raster::raster(env_variables),species.prevalence = 0.1)
virtual_species2 <- virtualspecies::generateRandomSp(raster::raster(env_variables),species.prevalence = 0.3)
virtual_species3 <- virtualspecies::generateRandomSp(raster::raster(env_variables),species.prevalence = 0.5)

#compile into one raster
species_occ_true <- rast(
  list(
    sp1 = rast(virtual_species1$probability.of.occurrence),
    sp2 = rast(virtual_species2$probability.of.occurrence),
    sp3 = rast(virtual_species3$probability.of.occurrence))
)

#visualise
plot(species_occ_true)

```

## Sample


```{r}
get_samples <- function(env_var,sp_occ,n_deployments,n_samples){
  deployment_locations <- spatSample(env_var,
                                     size = n_deployments,
                                     as.points=T,
                                     xy=T,
                                     replace=T)
  
  #sample the species
  deployments <- terra::extract(sp_occ,deployment_locations,xy=T)
  
  #go long data frame
  deployments <- pivot_longer(deployments,cols = names(sp_occ),names_to="species",values_to="sp_occ")
  
  names(deployments)[1] <- "deployment_id"
  
  #repeat for n_samples
  samples <- do.call("rbind", replicate(n_samples, deployments, simplify = FALSE))
  
  #turn to T/F detection 
  samples$detect_test <- runif(nrow(samples))
  samples$detected <- samples$detect_test<samples$sp_occ
  
  #add a sample ID
  samples$sample_id <- 1:nrow(samples)
  
  samples %>% arrange(deployment_id)
}

#unstructured sampling
cit_sci_samples <- get_samples(env_variables,species_occ_true,200,1)
cit_sci_samples 

#AMI tram
ami_trap_samples <- get_samples(env_variables,species_occ_true,4,50)
ami_trap_samples 

```

Fit a model using ibis.iSDM


```{r}
#create a background raster layer that just has value 1 for all cells
background <- setValues(env_variables,1)[[1]]

# First we define a distribution object using the background layer
mod <- distribution(background)

# turn the citizen science samples into points (sf)
virtual_species <- cit_sci_samples %>% 
  filter(detected==T,
         species == "sp1") %>% 
  st_as_sf(coords =c("x","y"),crs = 4326)

predictors <- env_variables

# Then lets add species data to it. 
# This data needs to be in sf format and key information is that
# the model knows where occurrence data is stored (e.g. how many observations per entry) as
# indicated by the field_occurrence field.
mod <- add_biodiversity_poipo(mod, virtual_species,
                                      name = "Virtual test species",
                                      field_occurrence = "detected")

# Then lets add predictor information
# Here we are interested in basic transformations (scaling), but derivates (like quadratic)
# for now, but check options
mod <- add_predictors(mod, 
                      env = predictors,
                      transform = "scale", derivates = "none")

# Finally define the engine for the model
# This uses the default data currently backed in the model,
# !Note that any other data might require an adaptation of the default mesh parameters used by the engine!
mod <- engine_inlabru(mod)

# Print out the object to see the information that is now stored within
print(mod)

plot(mod$biodiversity)


fit <- train(mod,
             runname =  "Test INLA run",
             verbose = T #
             )

# Plot the mean of the posterior predictions
plot(fit, "mean")

#compare to original
plot(terra::rast(virtual_species1$probability.of.occurrence))

# Print out some summary statistics
summary(fit)
```