Skip to content

Commit

Permalink
Function for realising state as pres/abs
Browse files Browse the repository at this point in the history
  • Loading branch information
simonrolph committed Nov 8, 2023
1 parent 53bebaf commit 006f3f2
Show file tree
Hide file tree
Showing 14 changed files with 203 additions and 71 deletions.
4 changes: 3 additions & 1 deletion R/SimulationObjectClass.R
Original file line number Diff line number Diff line change
Expand Up @@ -4,18 +4,20 @@ setClass("SimulationObject",
background = "ANY",
state_env = "ANY",
state_target = "ANY",
state_target_realised = "ANY",
effort = "ANY",
detect = "ANY",
report = "ANY"
)
)

# Create a constructor for the SimulationObject class
SimulationObject <- function(background, state_env = NULL, state_target = NULL, effort = NULL, detect = NULL, report = NULL) {
SimulationObject <- function(background, state_env = NULL, state_target = NULL, state_target_realised= NULL, effort = NULL, detect = NULL, report = NULL) {
new("SimulationObject",
background = background,
state_env = state_env,
state_target = state_target,
state_target_realised = state_target_realised,
effort = effort,
detect = detect,
report = report
Expand Down
20 changes: 13 additions & 7 deletions R/sim_detect_equal.R
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
#' Simulates detection where all detections occur at an equal probability and identified correctly
#' Simulates detection / identification where all detections occur at an equal probability and identified correctly
#'
#' @param simulation_object an R object of class 'SimulationObject' containing all the necessary information for the simulation
#' @param prob a numeric probability of each target being detected
Expand All @@ -7,21 +7,27 @@
#' \dontrun{
#' sim_detect_equal()
#' }
sim_detect_equal <- function(simulation_object, prob = 1) {
sim_detect_equal <- function(simulation_object, prob = 0.5) {
background <- simulation_object@background
state_env <- simulation_object@state_env
state_target <- simulation_object@state_target
state_target <- simulation_object@state_target_realised
effort <- simulation_object@effort

detections_all <- data.frame()

for (i in 1:dim(state_target)[3]) {
#loop through each of the targets
for (i in 1:(dim(state_target)[3])) {
detections <- effort
detections$state <- terra::extract(state_target[[i]], effort)$abundance
detections$target <- i

detections$detected <- runif(nrow(detections)) < prob
detections$identified <- TRUE
detections$state_realised <- unname(terra::extract(state_target[[i]], effort,ID=F,raw=T))

#detect based probability value provided as argument
detections$detected <- detections$state_realised * (runif(nrow(detections)) < prob)

#in this basic example all are identified correctly
detections$identified_as <- detections$target
detections$identified_correct <- detections$identified_as==detections$target

detections_all <- rbind(detections_all, detections)
}
Expand Down
20 changes: 13 additions & 7 deletions R/sim_effort_uniform.R
Original file line number Diff line number Diff line change
@@ -1,24 +1,30 @@
#' Simulates a uniform effort
#'
#' @param simulation_object an R object of class 'SimulationObject' containing all the necessary information for the simulation
#' @param n_visits a number indicating the average number of visits
#' @param n_samplers a number indicating the number of samplers
#' @param n_visits a number indicating the average number of visits per sampler
#' @param n_sample_events a number of sampling events made per visit
#' @param replace a logical indicating whether to sample with replacement
#' @return An updated simulation object with the newly calculated effort in the correct slot
#' @examples
#' \dontrun{
#' sim_effort_uniform(simulation_object, 100, FALSE)
#' sim_effort_uniform(simulation_object, 100, 20,1 FALSE)
#' }
sim_effort_uniform <- function(simulation_object, n_visits = 100, replace = FALSE) {
state_target <- simulation_object@state_target
sim_effort_uniform <- function(simulation_object, n_samplers = 1, n_visits = 1, n_sample_events=1, replace = FALSE) {

visited_cells <- sample(terra::cells(state_target), size = n_visits, replace = replace)
#which cells are visits
state_target <- simulation_object@state_target
visited_cells <- rep(sample(terra::cells(state_target), size = n_samplers*n_visits, replace = replace),each = n_sample_events)

# capture data
sim_effort_points <- as.data.frame(terra::xyFromCell(state_target, visited_cells))
sim_effort_points$sampler_id <- 1
sim_effort_points$sampler_id <- rep(1:n_samplers,each = n_visits*n_sample_events)
sim_effort_points$visit_id <- rep(1:n_visits,n_samplers,each = n_sample_events)
sim_effort_points$sample_id <- rep(1:n_sample_events,n_samplers*n_visits)

sim_effort_points$cell_id <- visited_cells

effort_sf <- sf::st_as_sf(sim_effort_points, coords = c("x", "y"), crs = terra::crs(state_target))

simulation_object@effort <- effort_sf

# Return the updated simulation_object
Expand Down
27 changes: 27 additions & 0 deletions R/sim_state_target_realise_binomial.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,27 @@
#' Realizes the state_target into binary from the continuous probability values
#'
#' @param sim_obj a SimulationObject
#' @return A SimulationObject with a binary state_target
#' @examples
#' \dontrun{
#' sim_state_target_binary(sim_obj)
#' }
sim_state_target_realise_binomial <- function(sim_obj) {
state_target <- sim_obj@state_target
binary_state_target <- state_target

for (i in 1:dim(sim_obj@state_target)[3]){
# Get the probability values from the state target
prob_values <- terra::values(state_target[[i]])

# Simulate binary values from the binomial distribution based on the probability values
binary_values <- rbinom(length(prob_values), 1, prob_values)

terra::values(binary_state_target[[i]]) <- binary_values
}

# Update the SimulationObject with the binary state_target
sim_obj@state_target_realised <- binary_state_target

sim_obj
}
10 changes: 5 additions & 5 deletions R/sim_state_target_uniform.R
Original file line number Diff line number Diff line change
@@ -1,18 +1,18 @@
#' Simulates a uniform state of the target
#'
#' @param simulation_object an R object of class 'SimulationObject' containing all the necessary information for the simulation
#' @param abundance a number indicating the abundance of the target in each cell
#' @param value a number indicating the value of the target in each cell
#' @return An updated simulation object with the newly calculated state of the target in the correct slot
#' @examples
#' \dontrun{
#' sim_state_target_uniform(simulation_object, 42)
#' sim_state_target_uniform(simulation_object, 0.5)
#' }
sim_state_target_uniform <- function(simulation_object, abundance = 10) {
sim_state_target_uniform <- function(simulation_object, value = 0.5) {
background <- simulation_object@background

sim_state <- background[[1]]
terra::values(sim_state) <- abundance
names(sim_state) <- "abundance"
terra::values(sim_state) <- value
names(sim_state) <- "target_1"

simulation_object@state_target <- sim_state

Expand Down
21 changes: 15 additions & 6 deletions README.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,8 @@ The `simulation_object` includes the following components:

- `@background`: Background extent and resolution of the simulated reality
- `@state_env`: Simulated state of the environment
- `@state_target`: Simulated state of the target
- `@state_target`: Simulated state of the target (as a continuous probability)
- `@state_target_realised`: Simulated state of the target (as a realised absolute/binary value)
- `@effort`: Simulated sampling effort
- `@detect`: Simulated detection information
- `@report`: Simulated reporting information
Expand All @@ -56,9 +57,10 @@ The functions used at each stage are as follows:

* `sim_state_env_...(simulation_object, ...)`
* `sim_state_target_...(simulation_object, ...)`
* `sim_state_target_realise_...(simulation_object, ...)`
* `sim_effort_...(simulation_object, ...)`
* `sim_detect_...(simulation_object, ...)`
* `sim_detect_...(simulation_object, ....)`
* `sim_report_...(simulation_object, ....)`

You could use the `targets` R package to create reproducible workflows for simulating your data.

Expand Down Expand Up @@ -122,7 +124,7 @@ library(terra)
library(sf)
# Create the background
background <- terra::rast(matrix(0,1000,600))
background <- terra::rast(matrix(0,30,30))
# Create the simulation object
sim_obj <- SimulationObject(background = background)
Expand All @@ -131,7 +133,10 @@ sim_obj <- SimulationObject(background = background)
sim_obj <- sim_state_env_gradient(sim_obj)
# Simulate the target state
sim_obj <- sim_state_target_uniform(sim_obj, abundance = 42)
sim_obj <- sim_state_target_uniform(sim_obj, value = 0.5)
#realise the state
sim_obj <- sim_state_target_realise_binomial(sim_obj)
# Simulate the sampling effort
sim_obj <- sim_effort_uniform(sim_obj, n_visits = 100, replace = FALSE)
Expand All @@ -142,8 +147,12 @@ sim_obj <- sim_detect_equal(sim_obj, prob = 0.5)
# Simulate the reporting
sim_obj <- sim_report_equal(sim_obj, prob = 0.8, platform = "iRecord")
plot(sim_obj@state_target) # State of the target
plot(sim_obj@effort$geometry, add = TRUE) # Effort
plot(sim_obj@state_target_realised) # State of the target
plot(sim_obj@effort$geometry, add = TRUE,pch=16) # Effort
plot(sim_obj@detect$geometry[sim_obj@detect$detected == FALSE], col = "red", pch = 4, add = TRUE) # Highlight the non-detections
plot(sim_obj@report$geometry[sim_obj@report$reported], col = "yellow", add = TRUE) # Highlight reported records as yellow
#add a legend
legend(1, 5, legend=c("Sampled", "Not detected","Reported"),
col=c("black","red", "yellow"), pch=c(16,4,1), cex=0.8,bg='grey')
```
93 changes: 66 additions & 27 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -33,25 +33,39 @@ Install from GitHub

## How to use the R package

For each of the 5 processes there are choices of functions to use
depending on your need. For each the processes there is the most basic
For each of the 5 processes, there are choices of functions to use
depending on your needs. For each process, there is the most basic
version for demonstration purposes.

The functions all follow this basic schema whereby all the objects from
the previous stage are arguments in the subsequent functions, whether or
not they are actually used in the calculations within the function:

- `state_env <- sim_state_env_______(background)`
- `state_target <- sim_state_target____(background, state_env, ...)`
- `effort <- sim_effort__________(background, state_env, state_target, ...)`
- `detect <- sim_detect__________(background, state_env, state_target, effort, ...)`
- `report <- sim_detect__________(background, state_env, state_target, effort, detect, ...)`

There are no species STRIDER R objects, this is intentional as to allow
flexibility ad interoperability. The outputs at each step are `terra`
SpatRasters or `sf` feature collections (POINT) see figure above, so if
you can use custom R scripts to generate the outputs of any of the
steps.
the previous stage, along with the `background` object, are combined
into a single `simulation_object`. This object is then used as an
argument in the subsequent functions, whether or not they are actually
used in the calculations within the function.

The `simulation_object` includes the following components:

- `@background`: Background extent and resolution of the simulated
reality
- `@state_env`: Simulated state of the environment
- `@state_target`: Simulated state of the target
- `@effort`: Simulated sampling effort
- `@detect`: Simulated detection information
- `@report`: Simulated reporting information

You can access and manipulate the `simulation_object` at each step to
generate the outputs of the corresponding processes. The outputs at each
step are `terra` SpatRasters or `sf` feature collections (POINT), as
shown in the figure above. You can use custom R scripts to generate the
outputs of any of the steps, ensuring flexibility and interoperability.

The functions used at each stage are as follows:

- `sim_state_env_...(simulation_object, ...)`
- `sim_state_target_...(simulation_object, ...)`
- `sim_effort_...(simulation_object, ...)`
- `sim_detect_...(simulation_object, ...)`
- `sim_detect_...(simulation_object, ....)`

You could use the `targets` R package to create reproducible workflows
for simulating your data.
Expand Down Expand Up @@ -158,17 +172,42 @@ library(sf)
## Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE

``` r
background <- terra::rast(matrix(0,1000,600)) # create background
state_env <- sim_state_env_gradient(background) #environment
state_target <- sim_state_target_uniform(background,state_env,42) #target
effort <- sim_effort_uniform(background,state_env,state_target,n_visits=100,replace=F) #effort
detections <-sim_detect_equal(background,state_env,state_target,effort,prob=0.5) #detection
reports <- sim_report_equal(background,state_env,state_target,effort,detections,prob=0.8,platform="iRecord") #reports

plot(state_target) #state of target
plot(effort$geometry,add=T) #effort
plot(detections$geometry[detections$detected==F],col="red",pch=4,add=T) #highlight the non-detections
plot(reports$geometry[reports$reported],col="yellow",add=T) # highlight reported records as yellow
# Create the background
background <- terra::rast(matrix(0,30,30))

# Create the simulation object
sim_obj <- SimulationObject(background = background)

# Simulate the environment state
sim_obj <- sim_state_env_gradient(sim_obj)

# Simulate the target state
sim_obj <- sim_state_target_uniform(sim_obj, value = 0.5)

#realise the state
sim_obj <- sim_state_target_realise_binomial(sim_obj)

# Simulate the sampling effort
sim_obj <- sim_effort_uniform(sim_obj, n_visits = 100, replace = FALSE)

# Simulate the detection
sim_obj <- sim_detect_equal(sim_obj, prob = 0.5)
```

## [1] 1

``` r
# Simulate the reporting
sim_obj <- sim_report_equal(sim_obj, prob = 0.8, platform = "iRecord")

plot(sim_obj@state_target_realised) # State of the target
plot(sim_obj@effort$geometry, add = TRUE,pch=16) # Effort
plot(sim_obj@detect$geometry[sim_obj@detect$detected == FALSE], col = "red", pch = 4, add = TRUE) # Highlight the non-detections
plot(sim_obj@report$geometry[sim_obj@report$reported], col = "yellow", add = TRUE) # Highlight reported records as yellow

#add a legend
legend(1, 5, legend=c("Sampled", "Not detected","Reported"),
col=c("black","red", "yellow"), pch=c(16,4,1), cex=0.8,bg='grey')
```

![](README_files/figure-gfm/example-1.png)<!-- -->
Binary file modified README_files/figure-gfm/example-1.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
6 changes: 3 additions & 3 deletions man/sim_detect_equal.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

16 changes: 13 additions & 3 deletions man/sim_effort_uniform.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

22 changes: 22 additions & 0 deletions man/sim_state_target_realise_binomial.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

Loading

0 comments on commit 006f3f2

Please sign in to comment.