Effort slot output is in long format with rows for each target. Updat…

…ed readme and vignettes.

R/sim_effort.R

@@ -30,12 +30,29 @@ sim_effort <- function(simulation_object, fun, sf=NULL, ...) {
   }
 
   #get values from env, suitability, realised
-  extracted_values <- terra::extract(simulation_object@state_env,effort_sf)
+  extracted_values <- terra::extract(simulation_object@state_env,effort_sf,ID=F)
   effort_sf[,names(extracted_values)] <- extracted_values
-  extracted_values <- terra::extract(simulation_object@state_target_suitability,effort_sf)
-  effort_sf[,paste0("suit_",names(extracted_values))] <- extracted_values
-  extracted_values <- terra::extract(simulation_object@state_target_realised,effort_sf)
-  effort_sf[,paste0("real_",names(extracted_values))] <- extracted_values
+
+  #loop through each target
+  targets_sf <- list()
+  for (target in names(sim_obj@state_target_suitability)){
+    target_sf <- effort_sf
+    target_sf$target <- target
+
+    #extract suitability
+    extracted_values <- terra::extract(simulation_object@state_target_suitability[target],effort_sf,ID=F)
+
+    target_sf[,"target_suitability"] <- extracted_values
+
+    #extract realised
+    extracted_values <- terra::extract(simulation_object@state_target_realised[target],effort_sf,ID=F)
+    target_sf[,"target_realised"] <- extracted_values
+
+    targets_sf[[target]] <- target_sf
+  }
+
+  effort_sf <- do.call(rbind,targets_sf)
+  rownames(effort_sf) <- NULL
 
   # validity checks
   fun_args <- as.list(match.call())

README.Rmd

@@ -104,7 +104,7 @@ background <- terra::rast(matrix(0,30,30))
 
 # Create the simulation object
 sim_obj <- SimulationObject(background = background)
-sim_obj <- sim_state_env(sim_obj,fun = "uniform",value = 20)
+sim_obj <- sim_state_env(sim_obj,fun = state_env_gradient,from = 0,to = 1)
 sim_obj
 ```
 
@@ -129,7 +129,15 @@ The BYOD (Bring Your Own Data) function is `sim_state_env(spatraster = [your_ras
 Here we define the state of the target or targets. We define two versions of this: a continuous variable representing a probability of occurrence (slot `@state_target`), and a realised absolute value (slot `@state_target_realised`) which could contain a binary (0 or 1) representing species occupancy or a positive integer representing abundance. Here's an example: 
 
 ```{r realisation_diagram, echo=F}
-sim_obj <- sim_state_target_suitability(sim_obj,fun = "uniform")
+# optimal environment is 0.5
+state_target_suitability_example <- function(sim_obj){
+  out <-  sim_obj@state_env*2
+  out[out>1] <- 2-out[out>1]
+  out
+}
+
+
+sim_obj <- sim_state_target_suitability(sim_obj,state_target_suitability_example)
 sim_obj <- sim_state_target_realise(sim_obj,fun = "binomial")
 plot(sim_obj)
 ```
@@ -215,7 +223,7 @@ suit_fun <- function(sim_obj){
   names(target_suitability) <- "frog" #give my layer a name
 
   # set suitability under certain critera eg. 0.7 when rainfall>500 and altitude < 50
-  target_suitability[sim_obj@state_env$rainfall>500 & sim_obj@state_env$altitude<50] <- 0.4
+  target_suitability[[sim_obj@state_env$rainfall>500 & sim_obj@state_env$altitude<50]] <- 0.4
   target_suitability[sim_obj@state_env$rainfall>800 & sim_obj@state_env$altitude<40] <- 0.9
 
   target_suitability #return just the suitability layer
@@ -231,21 +239,27 @@ This function must take the SimulationObject as its first argument. This means y
 Here is an example which runs through a very simple example and plots the output. 
 
 ```{r example}
+rm(sim_obj)
 library(STRIDER)
 library(terra)
 library(sf)
 
 # Create the background
-background <- terra::rast(matrix(0,30,30))
+background <- terra::rast(matrix(0,50,50))
 
 # Create the simulation object
 sim_obj <- SimulationObject(background = background)
 
 # Simulate the environment state
-sim_obj <- sim_state_env(sim_obj,fun="gradient")
+sim_obj <- sim_state_env(sim_obj,fun = state_env_gradient,from = 0,to = 1)
 
 # Simulate the target state
-sim_obj <- sim_state_target_suitability(sim_obj,fun = "uniform", value = 0.5)
+state_target_suitability_example <- function(sim_obj){ 
+  out <-  sim_obj@state_env*2
+  out[out>1] <- 2-out[out>1]
+  out # optimal environment is 0.5
+}
+sim_obj <- sim_state_target_suitability(sim_obj,state_target_suitability_example)
 
 #realise the state
 sim_obj <- sim_state_target_realise(sim_obj,fun = "binomial")

README.md

@@ -144,7 +144,7 @@ background <- terra::rast(matrix(0,30,30))
 
 # Create the simulation object
 sim_obj <- SimulationObject(background = background)
-sim_obj <- sim_state_env(sim_obj,fun = "uniform",value = 20)
+sim_obj <- sim_state_env(sim_obj,fun = state_env_gradient,from = 0,to = 1)
 sim_obj
 ```
 
@@ -168,8 +168,8 @@ sim_obj
     ## coord. ref. :  
     ## source(s)   : memory
     ## name        : env 
-    ## min value   :  20 
-    ## max value   :  20 
+    ## min value   :   0 
+    ## max value   :   1 
     ## 
     ## Slot "state_target_suitability":
     ## NULL
@@ -189,15 +189,18 @@ sim_obj
     ## Slot "metadata":
     ## $state_env
     ## $state_env$fun
-    ## [1] "uniform"
+    ## state_env_gradient
     ## 
-    ## $state_env$value
-    ## [1] 20
+    ## $state_env$from
+    ## [1] 0
+    ## 
+    ## $state_env$to
+    ## [1] 1
     ## 
     ## 
     ## 
     ## Slot "hash":
-    ## [1] "253f8b7c8478a288ccd2b954fa3f1196"
+    ## [1] "841156773684e37801aefddb92881877"
 
 ### Simulating the environmental state
 
@@ -311,21 +314,21 @@ sim_obj <- sim_effort(sim_obj,fun = "basic", n_samplers=2, n_visits = 1, n_sampl
 sim_obj@effort
 ```
 
-    ## Simple feature collection with 4 features and 10 fields
+    ## Simple feature collection with 4 features and 8 fields
     ## Geometry type: POINT
     ## Dimension:     XY
-    ## Bounding box:  xmin: 11.5 ymin: 5.5 xmax: 27.5 ymax: 13.5
+    ## Bounding box:  xmin: 4.5 ymin: 8.5 xmax: 11.5 ymax: 13.5
     ## CRS:           NA
-    ##   sampler visit unit cell_id          geometry ID env suit_ID suit_target_1
-    ## 1       1     1    1     492 POINT (11.5 13.5)  1  20       1           0.5
-    ## 2       1     1    2     492 POINT (11.5 13.5)  2  20       2           0.5
-    ## 3       2     1    1     748  POINT (27.5 5.5)  3  20       3           0.5
-    ## 4       2     1    2     748  POINT (27.5 5.5)  4  20       4           0.5
-    ##   real_ID real_target_1
-    ## 1       1             1
-    ## 2       2             1
-    ## 3       3             0
-    ## 4       4             0
+    ##   sampler visit unit cell_id       env target target_suitability
+    ## 1       1     1    1     492 0.3793103    env          0.7586207
+    ## 2       1     1    2     492 0.3793103    env          0.7586207
+    ## 3       2     1    1     635 0.1379310    env          0.2758621
+    ## 4       2     1    2     635 0.1379310    env          0.2758621
+    ##   target_realised          geometry
+    ## 1               1 POINT (11.5 13.5)
+    ## 2               1 POINT (11.5 13.5)
+    ## 3               0   POINT (4.5 8.5)
+    ## 4               0   POINT (4.5 8.5)
 
 The function for simulating effort start is `sim_effort`
 
@@ -388,7 +391,7 @@ suit_fun <- function(sim_obj){
   names(target_suitability) <- "frog" #give my layer a name
 
   # set suitability under certain critera eg. 0.7 when rainfall>500 and altitude < 50
-  target_suitability[sim_obj@state_env$rainfall>500 & sim_obj@state_env$altitude<50] <- 0.4
+  target_suitability[[sim_obj@state_env$rainfall>500 & sim_obj@state_env$altitude<50]] <- 0.4
   target_suitability[sim_obj@state_env$rainfall>800 & sim_obj@state_env$altitude<40] <- 0.9
 
   target_suitability #return just the suitability layer
@@ -408,21 +411,27 @@ Here is an example which runs through a very simple example and plots
 the output.
 
 ``` r
+rm(sim_obj)
 library(STRIDER)
 library(terra)
 library(sf)
 
 # Create the background
-background <- terra::rast(matrix(0,30,30))
+background <- terra::rast(matrix(0,50,50))
 
 # Create the simulation object
 sim_obj <- SimulationObject(background = background)
 
 # Simulate the environment state
-sim_obj <- sim_state_env(sim_obj,fun="gradient")
+sim_obj <- sim_state_env(sim_obj,fun = state_env_gradient,from = 0,to = 1)
 
 # Simulate the target state
-sim_obj <- sim_state_target_suitability(sim_obj,fun = "uniform", value = 0.5)
+state_target_suitability_example <- function(sim_obj){ 
+  out <-  sim_obj@state_env*2
+  out[out>1] <- 2-out[out>1]
+  out # optimal environment is 0.5
+}
+sim_obj <- sim_state_target_suitability(sim_obj,state_target_suitability_example)
 
 #realise the state
 sim_obj <- sim_state_target_realise(sim_obj,fun = "binomial")

tests/testthat/test-simulation_basic.R

@@ -27,8 +27,6 @@ sim_obj <- sim_detect(sim_obj,fun="equal", prob = 0.5)
 # 5 Simulate reporting within the simulation object
 sim_obj <- sim_report(sim_obj, fun="equal", prob = 0.8, platform = "iRecord")
 
-plot(sim_obj)
-
 #1
 test_that("Test creating a uniform environment", {
   expect_true(class(sim_obj) == "SimulationObject")

vignettes/custom_functions.Rmd

@@ -0,0 +1,159 @@
+---
+title: "More complex simulations"
+output: 
+  rmarkdown::html_vignette:
+    fig_width: 7
+    fig_height: 6
+vignette: >
+  %\VignetteIndexEntry{More complex simulations}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>"
+)
+```
+
+## Introduction
+
+This vignette demonstrates a simulation workflow to investigate the benefits and issues of co-locating recording sites for different taxonomic groups, such as butterflies, plants, and birds. We will generate species counts/abundance data across multiple sites, simulate dependencies between species, and incorporate temporal trends. This approach helps assess whether using the same recording sites across different taxa groups is beneficial or problematic and explore optimal strategies for site co-location.
+
+## Setup
+
+First, we load the required packages and set a seed for reproducibility.
+
+```{r setup}
+library(STRIDER)
+library(terra)
+library(sf)
+library(dplyr)
+library(ggplot2)
+set.seed(42)
+```
+
+## Custom Environmental State
+
+We define custom environmental variables with different spatial patterns to simulate realistic scenarios.
+
+```{r custom environmental state}
+dim_x <- 500
+dim_y <- 500
+
+# Define environmental variables with different spatial patterns
+env1 <- terra::rast(matrix(rep(seq(from = 1, to = dim_x), times = dim_y), dim_x, dim_y))
+env1 <- env1 / max(values(env1))
+
+env2 <- terra::rast(matrix(rep(seq(from = dim_x, to = 1), each = dim_y), dim_x, dim_y))
+env2 <- env2 / max(values(env2))
+
+env3 <- terra::rast(matrix(runif(dim_x * dim_y), dim_x, dim_y))
+env3 <- env3 / max(values(env3))
+
+env4 <- terra::rast(matrix(sin(seq(0, pi, length.out = dim_x)), dim_x, dim_y))
+env4 <- env4 / max(values(env4))
+
+env5 <- terra::rast(matrix(cos(seq(0, pi, length.out = dim_y)), dim_x, dim_y))
+env5 <- env5 / max(values(env5))
+
+custom_env <- c(env1, env2, env3, env4, env5)
+names(custom_env) <- c("rainfall", "temperature", "urban_density", "elevation", "aspect")
+```
+
+Next, we integrate these custom environmental variables into the simulation object.
+
+```{r integrate custom environment}
+background <- terra::rast(matrix(0, dim_x, dim_y))
+sim_obj <- SimulationObject(background = background)
+sim_obj <- sim_state_env(sim_obj, spatraster = custom_env)
+```
+
+## Custom Suitability Functions for Multiple Species
+
+We define custom suitability functions for multiple species, each influenced differently by the environmental variables.
+
+```{r custom suitability functions}
+# Custom suitability function for multiple species
+custom_suitability_function <- function(simulation_object) {
+  # Extract environmental variables from the simulation object
+  env <- simulation_object@state_env
+  
+  # Define the suitability functions for multiple species
+  species_suitability_functions <- list(
+    butterfly = function(env) (env$rainfall * env$temperature + env$urban_density^2) / 3,
+    plant = function(env) (env$rainfall + env$temperature^2 - env$urban_density) / 3,
+    bird = function(env) (env$rainfall + env$elevation + env$aspect) / 3,
+    mammal = function(env) (env$temperature * env$elevation - env$urban_density) / 3,
+    amphibian = function(env) (env$rainfall + env$elevation^2 - env$aspect) / 3
+  )
+  
+  # Initialize an empty list to store the suitability layers
+  suitability_layers <- list()
+  
+  # Loop through each species and calculate its suitability layer
+  for (species in names(species_suitability_functions)) {
+    suitability <- species_suitability_functions[[species]](env)
+    suitability[suitability<0] <- 0
+    suitability[suitability>1] <- 1
+    names(suitability) <- species
+    suitability_layers[[species]] <- suitability
+  }
+  
+  # Combine all suitability layers into a single SpatRaster object
+  suitability_raster <- rast(suitability_layers)
+  
+  return(suitability_raster)
+}
+
+sim_obj <- sim_state_target_suitability(sim_obj, fun = custom_suitability_function)
+```
+
+Realised suitability
+
+```{r}
+sim_obj <- sim_state_target_realise(sim_obj,fun="binomial")
+```
+
+## Custom Effort Function
+
+We create a custom effort function that simulates sampling effort using a combination of urban density and random sampling.
+
+```{r custom effort function}
+custom_effort_function <- function(sim_obj, n_sites, n_taxa_per_site = 1) {
+  sites <- sample(cells(sim_obj@background), n_sites, replace = TRUE)
+  coords <- xyFromCell(sim_obj@background, sites)
+  effort_df <- data.frame(
+    x = coords[, 1],
+    y = coords[, 2]
+  )
+  effort_sf <- st_as_sf(effort_df, coords = c("x", "y"))
+  
+  site_taxa_combination <- expand.grid(sites = 1:n_sites,target = names(sim_obj@state_target_suitability))
+  
+  effort_sf$taxa <- sample(names(sim_obj@state_target_suitability), n_taxa_per_site)
+  
+  return(effort_sf)
+}
+
+sim_obj <- sim_effort(sim_obj, fun = custom_effort_function, n_sites = 100)
+```
+
+## Custom Detection and Reporting Functions
+
+We could define custom detection and reporting functions to introduce variability in the detection and reporting probabilities, but just use the default for now.
+
+```{r custom detection and reporting functions}
+# Simulate the detection
+sim_obj <- sim_detect(sim_obj,fun = "equal", prob = 0.5)
+
+# Simulate the reporting
+sim_obj <- sim_report(sim_obj,fun = "equal", prob = 0.8, platform = "iRecord")
+```
+
+## Visualisation
+
+```{r visualisation}
+plot(sim_obj)
+```