metaRange is a framework to build a variety of different process based species distribution models that can include a (basically) arbitrary number of environmental factors, processes, species and species interactions. The common denominator for all models built with metaRange is that they are raster and population based.

Each metaRange simulation object contains one environment that holds and manages all the environmental factors that may influence the simulation (as raster data) and one or more species that are simulated in the environment.

Each species object has two main characteristics: traits and processes. Species traits can be (somewhat arbitrary) pieces of data that describe and store information about the species while processes are functions that describe how the species interacts with itself, time, climate and other species. During each time step of the simulation, the processes of the species are executed in a user defined order and can access and modify the species traits.

While the models built with metaRange can be quite variable in their structure, they are all based on population dynamics. This means that in most cases a trait will not be a single number but a matrix that has the same size as the raster data in the environment, where each value represents the trait value for one population in the corresponding grid cell of the environment. Based on this, most processes will describe population (and meta population) dynamics and not individual based mechanisms. Figure 1 shows an example of some of the different environmental factors, species traits and processes that could be included in a simulation.

Figure 1: Example of some of the different environmental factors, species traits and processes that could be included in a simulation.

A more technical overview of the different components of a simulation and how they interact with each other is shown in Figure 2. While the environment has a temporal dimension (i.e. the different time steps, see Fig.1), the traits (or the global variables of the simulation) have no temporal property. They represent the state of the species (or of the simulation respectively) in one (the current) time step. When a process is executed within a time step, it can access this current state of the simulation and modify it, which results in changing traits over the course of the simulation. Each process can be assigned a priority that is used by the process priority queue to determine in which order the processes are executed within one time step.

Figure 2: Overview of the different components of a simulation and how they interact with each other. Note that the number of species as well as the number of traits and processes per species is not limited, but only a selection is shown for simplicity.

Setting up a simulation

Following is a simple example of how to set up a simulation with metaRange, in which we only use a single species and one environmental factor (habitat quality). At the end of this introduction we will see how the abundance of the species changes in relation to the quality of the habitat each population occupies.

To start, we need to load the packages.

library(metaRange) # does the simulation
#> metaRange version: 1.1.4
library(terra) # handles the raster data processing
#> terra 1.7.83

Loading the landscape

The first step when setting up a simulation is the loading of the environment in which the simulation will take place. This can either be real world data or “theoretical” / generated data and may include for example different climate variables, land cover or elevation.

The simulation expects this data as an SpatRasterDataset (SDS) which is a collection of different raster files that all share the same extent and resolution. Each sub-dataset in this SDS represents one environmental variable and each layer represents one time step of the simulation. In other words, metaRange does not simulate the environmental conditions itself, but expects the user to provide the environmental data for each time step.

To create such a dataset one can use the function terra::sds(). One important note: Since each layer represents the condition in one time step, all the raster files that go into the SDS need to have the same number of layers (i.e. the desired number of time steps the simulation should have). After the SDS is created, the individual sub-datasets should be named, since this is how the simulation will refer to them.

To simplify this introduction, we use an example landscape consisting only of habitat quality data, with 10 time steps (layers) that are all the same (i.e. no environmental change). Luckily the terra package has a built-in demo that we can use for this purpose.

# find the file
raster_file <- system.file("ex/elev.tif", package = "terra")

# load it
r <- rast(raster_file)

# scale it
r <- scale(r, center = FALSE, scale = TRUE)
plot(r, main = "Habitat quality")

Figure 3: The habitat quality of the example landscape. Note: higher value = better habitat quality

Now we can turn this raster with one layer into an SDS that has multiple layer (one for each time step).

r <- rep(r, 10)
landscape <- sds(r)
names(landscape) <- c("habitat_quality")
landscape
#> class       : SpatRasterDataset 
#> subdatasets : 1 
#> dimensions  : 90, 95 (nrow, ncol)
#> nlyr        : 10 
#> resolution  : 0.008333333, 0.008333333  (x, y)
#> extent      : 5.741667, 6.533333, 49.44167, 50.19167  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source(s)   : memory 
#> names       : habitat_quality

Pre-setup

Before creating the simulation, it may be helpful to enable extensive reporting, which will print out a lot of information each time a metaRange function is called. This can be enabled or disabled at any time (i.e. also while the simulation is running), but in order to highlight what each function call in this tutorial does, we enable it at the beginning of the setup.

# 0 = no reporting
# 1 = a bit of info
# 2 = very verbose
set_verbosity(2)

Creating the simulation

After the landscape is loaded, the simulation can be created using the create_simulation() function. The only required argument is source_environment which is the landscape / environment SDS that was created in the first step. One can optionally specify an ID for the simulation and a seed for the random number generator.

sim <- create_simulation(
    source_environment = landscape,
    ID = "example_simulation",
    seed = 1
)
#> number of time steps: 10
#> time step layer mapping: 1, 2, 3, 4, 5, 6, 7, 8, 9, 10
#> added environment
#> class       : SpatRasterDataset 
#> subdatasets : 1 
#> dimensions  : 90, 95 (nrow, ncol)
#> nlyr        : 10 
#> resolution  : 0.008333333, 0.008333333  (x, y)
#> extent      : 5.741667, 6.533333, 49.44167, 50.19167  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source(s)   : memory 
#> names       : habitat_quality
#> 
#> created simulation: example_simulation

If you want to inspect the simulation object, you can either print it, to lists its fields and methods or use the summary() function to get an overview of the simulation state.

sim
#> metaRangeSimulation object
#> Fields: 
#>   $ID
#>   $globals
#>   $environment
#>   $number_time_steps
#>   $time_step_layer
#>   $current_time_step
#>   $queue
#>   $processes
#>   $seed
#> Species: none
#> Methods: 
#>   $species_names()
#>   $add_globals()
#>   $add_species()
#>   $add_traits()
#>   $add_process()
#>   $begin()
#>   $exit()
#>   $set_current_time_step()
#>   $set_time_layer_mapping()
#>   $print()
#>   $summary()
summary(sim)
#> ID: example_simulation 
#> Environment: 
#> Fields: 
#> $current ==== the environment at the current time step
#> classes     : all -> matrix
#> number      : 1
#> names       : habitat_quality
#> $sourceSDS == the source raster data of the environment
#> class       : SpatRasterDataset 
#> subdatasets : 1 
#> dimensions  : 90, 95 (nrow, ncol)
#> nlyr        : 10 
#> resolution  : 0.008333333, 0.008333333  (x, y)
#> extent      : 5.741667, 6.533333, 49.44167, 50.19167  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source(s)   : memory 
#> names       : habitat_quality 
#> Time step layer mapping:  1 2 3 4 5 6 7 8 9 10 
#> Current time step:  1 
#> Seed:  1 
#> Species: 0
#> 
#> Simulation level processes:
#> NULL
#> Gobal variables:
#> NULL
#> Queue:
#> Remaining queue (this time step):  0 
#>  NULL
#> Future queue (next time step):  0 
#>  NULL

Adding species to the simulation

After the simulation is created, species can be added to it using the add_species() function. At this point one has to switch to the syntax of the R6 package that metaRange uses. This means that add_species() is a method of the simulation object and can be called using the $ operator (i.e. by indexing the simulation object and calling a function that is stored inside of it). The only required argument is name which is the name of the species that will be added.

sim$add_species("species_1")
#> adding species
#> name: species_1

This species can now be accessed by using the $ operator again.

sim$species_1
#> Species:  species_1 
#> processes: 
#> NULL
#> traits: 
#> character(0)

Adding traits to species

After species are added to the simulation, traits can be assigned to them using the add_traits() method. The first argument is species which is a character vector of the species names to which the trait should be assigned to. The second argument is population_level, a TRUE/FALSE value, that decides if the trait should be stored with one value per population (i.e. as a matrix of the same size as the landscape) or not (i.e. only one value per species). All following arguments can be supplied in the form of trait_name = trait_value.

For now we only add three traits: abundance (number of individuals in each population), reproduction_rate (how fast the populations can reproduce) and carrying_capacity (maximum number of individuals per grid cell). Note: Traits always represent the “current” state of a species. This means that the abundance we use as input here represents the initial state of the simulation. Over the course of the simulation (i.e. in each time step) the traits can be updated and changed. In this example, the abundance will change each time step while e.g. the reproduction rate stays the same, but in other cases each trait might change with time.

sim$add_traits(
    species = "species_1",
    population_level = TRUE,
    abundance = 100,
    reproduction_rate = 0.5,
    carrying_capacity = 1000
    # ...
    # Note that here could be more traits, there is no limit
)
#> adding traits:
#> [1] "abundance"         "reproduction_rate" "carrying_capacity"
#> 
#> to species:
#> [1] "species_1"
#>

We can check what traits a species has by printing them:

sim$species_1$traits
#> abundance :  num [1:90, 1:95] 100 100 100 100 100 100 100 100 100 100 ...
#> carrying_capacity :  num [1:90, 1:95] 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 ...
#> reproduction_rate :  num [1:90, 1:95] 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 ...

Or plotting them:

plot(sim$species_1, "abundance")

Figure 4: The initial abundance of the species.

Note that the above plot is not very interesting, since the abundance is the same for each population at the beginning of the simulation.

Adding processes

After the species and its traits are added, the processes that describe how the species interacts with its environment can be added, using the add_process() method. The arguments are: species which is a character vector of the species (names) that should receive the process, process_name which is a human readable name for the process and process_fun which is the function that will be called when the process is executed.

One argument that might be confusing is the execution_priority. This is a number that gives the process a priority “weight” and decides in which order the processes are executed within one time step. The smaller the number, the earlier the process will be executed (e.g. 1 gets executed before 2). In the case two (or more) processes have the same priority, it is assumed that they are independent from each other and that their execution order does not matter.

Reproduction

In this example we will only add a single process (reproduction) to the species, that is going to calculate the abundance (for each population) in the next time step, depending on the habitat quality. To do so, we can use a built-in function ricker_reproduction_model() that implements the “classic” Ricker reproduction model (Ricker, W.E. (1954)) [Ref. 1], which describes the population dynamics of a species with non-overlapping generations in discrete time steps. This model features density dependent growth and possibly overcompensatory dynamics (i.e. the populations can, depending on the reproduction rate, become larger than the carrying capacity, which inevitably leads to a decline in the next time step).

Note the use of the self keyword in the function. In this context, self refers to the species that the process is attached to. This means that the function can access the species traits and modify them and also access the environment (each species holds a reference to the simulation it was created in).

sim$add_process(
    species = "species_1",
    process_name = "reproduction",
    process_fun = function() {
        # use a ricker reproduction model
        # to calculate the new abundance
        # and let the carrying capacity
        # depend on the habitat quality
        ricker_reproduction_model(
            self$traits$abundance,
            self$traits$reproduction_rate,
            self$traits$carrying_capacity * self$sim$environment$current$habitat_quality
        )

        # print out the current mean abundance
        print(
            paste0("mean abundance: ", mean(self$traits$abundance))
        )
    },
    execution_priority = 1
)
#> adding process: reproduction
#> to species:
#> [1] "species_1"
#>

Executing the simulation

After the species, traits and processes are added to the simulation, it can be executed via the begin() method.

sim$begin()
#> Starting simualtion.
#> passed initial sanity checks.
#> start of time step: 1
#> |- species_1 : reproduction
#> [1] "mean abundance: 84.1732579542268"
#> |---- 0.00061 secs
#>  10 % done | 0.062 secs remaining (estimate)
#> start of time step: 2
#> |- species_1 : reproduction
#> [1] "mean abundance: 127.598064388471"
#> |---- 0.00069 secs
#>  20 % done | 0.1 secs remaining (estimate)
#> start of time step: 3
#> |- species_1 : reproduction
#> [1] "mean abundance: 185.433176212344"
#> |---- 5e-04 secs
#>  30 % done | 0.04 secs remaining (estimate)
#> start of time step: 4
#> |- species_1 : reproduction
#> [1] "mean abundance: 254.974925085775"
#> |---- 5e-04 secs
#>  40 % done | 0.033 secs remaining (estimate)
#> start of time step: 5
#> |- species_1 : reproduction
#> [1] "mean abundance: 328.335965177599"
#> |---- 5e-04 secs
#>  50 % done | 0.028 secs remaining (estimate)
#> start of time step: 6
#> |- species_1 : reproduction
#> [1] "mean abundance: 394.79908318955"
#> |---- 6e-04 secs
#>  60 % done | 0.041 secs remaining (estimate)
#> start of time step: 7
#> |- species_1 : reproduction
#> [1] "mean abundance: 446.224976588256"
#> |---- 0.00058 secs
#>  70 % done | 0.018 secs remaining (estimate)
#> start of time step: 8
#> |- species_1 : reproduction
#> [1] "mean abundance: 480.705159122178"
#> |---- 0.00055 secs
#>  80 % done | 0.012 secs remaining (estimate)
#> start of time step: 9
#> |- species_1 : reproduction
#> [1] "mean abundance: 501.356439544315"
#> |---- 0.00051 secs
#>  90 % done | 0.0056 secs remaining (estimate)
#> start of time step: 10
#> |- species_1 : reproduction
#> [1] "mean abundance: 512.803082103058"
#> |---- 5e-04 secs
#> 100 % done | 0 secs remaining (estimate)
#> 
#> Simulation: 'example_simulation' finished
#> Exiting the Simulation
#> Runtime: 0.077 secs

Plotting the results

To investigate the results, you can use the plot() function.

# define a nice color palette
plot_cols <- hcl.colors(100, "BluYl", rev = TRUE)
plot(
    sim,
    obj = "species_1", # name of the species
    name = "abundance", # name of the trait to plot
    main = "Species 1: abundance", # optional title
    col = plot_cols # color palette
)

Figure 5: The resulting abundance distribution of species 1 after 10 simulation time steps.

Saving the simulation

To save results you can use the save_species() function. This will save the (possibly specified) traits of a species, either as a raster (.tif) or as a text (.csv) file, whatever is more appropriate for the data. Note that this function does not save the species processes. One should keep a copy of the script that is used to run the simulation to make it repeatable.

save_species(
    sim$species_1,
    traits = c("name", "of", "one_or_more", "traits"),
    path = "path/to/a/folder/"
)

References

Ricker, W.E. (1954) Stock and recruitment. Journal of the Fisheries Research Board of Canada, 11, 559–623. doi:10.1139/f54-039

01: Introduction to metaRange