ecoforge_Rpackage example

Authors
David Ferrer Ferrando, Pelayo Acevedo & Pedro Tarroso

Published
2024/03/19 - 2026/05/27

---

Case Study: Simulating virtual species and sampling designs with the ecoforge package

In this vignette we present a complete case study illustrating how to use the ecoforge package to build spatially explicit simulations of virtual species. This type of virtual ecology framework is particularly useful for evaluating spatial ecology models (such as species distribution models or occupancy models) under known and controlled conditions. For example evaluating the effect of different sampling strategies on model performance, or generating reproducible synthetic datasets for methodological research in ecology.

To facilitate the comprehension of the full simulation workflow and the possibilities offered by the package functions, the example is structured in the following steps:

1. Install and load the required packages (ecoforge and dependencies).

2. Generate the environmental spatial layers, and modify or transform them to define habitat characteristics of the study area.

3. Define the target species’ habitat suitability, combining environmental layers to build a spatial suitability surface.

4. Generate the species’ distribution, placing individuals across the study area according to habitat suitability, and optionally introducing possible distribution bias.

5. Sample the species’ distribution, designing the survey strategy (point counts or line transects, random or uniform sampling units), and obtaining count or presence/absence (or detection/non-detection) data.

6. Export the simulated data in formats ready for use in ecological models.

7. Additional tip assessing correlation between environmental variables before model fitting.

To see more about ecoforge package (and other vignettes) see ecoforge website

1. Getting started

First of all, install and load ecoforge package and all the dependencies required for the correct function of these simulations.

# install.packages("dismo")
# install.packages("fields")
# install.packages("devtools")
# devtools::install_github("NatureProtectionDFF/Ecoforge-R-package")
library(ecoforge)
library(dismo)
library(fields)

set.seed(1234) # To create reproducible results

2. Generate the environmental spatial layers

These functions let us control and determine some parameter values that define spatial variables. Inside of these functions, we employ a previous developed function gstat() of the package gstat.

2.1 Climate variables

We start creating clime conditions. To do this, we have a first function generateClim() that allows us to create climate variables (in this case, temperature and precipitation). When creating these variables, several parameters can be controlled, such as the mean value Tmean & Pmean, the range of variation Tsd & Psd, and the range of spatial autocorrelation Trange & Prange. For temperature, daily and annual ranges are also considered T_amp & T_ar, both with maximum and minimum range values. Precipitation also follows an annual range Pav. Both annual ranges follow a latitudinal gradient: in this case, the higher the latitude, the greater the range. To learn more about the different components of this function (and the following ones), read the description and details of the function available on the website, ecoforge Functions reference.

clim <- generateClim(nrow=200, ncol=200, Tmean=15, Tsd=4, Trange=45, 
                     Tamp_max=10, Tamp_min=1, Tar_max=10, Tar_min=5,
                     Pmean=500, Psd=300, Prange=45, Pav_max=50, 
                     Pav_min=10)

[using unconditional Gaussian simulation]
[using unconditional Gaussian simulation]
[using unconditional Gaussian simulation]

clim_rst <- rast(clim, type="xyz")
cols <- rev(terrain.colors(100))
plot(clim_rst, col= cols)

clim_rst <- c(clim_rst$temp, clim_rst$prec) #save a raster with only mean Temperature and mean precipitation

But what we have obtained is an average value in a fixed image (a single value per cell). In this fixed image, we have a single value per cell. To obtain a value that varies throughout the year for each cell, we create monthly values for temperature and precipitation with the generateMonthlyTemp() and generateMonthlyPrec() functions, respectively. In the case of temperature, we will obtain a maximum, average, and minimum temperature for each month. In the case of precipitation, we obtain the cumulative value for each month. The daily and annual high and low values are determined by the values assigned within the ranges we set before.

m_temp <- generateMonthlyTemp(clim$temp, clim$daily_amp, clim$annual_range, err_sdev=1)

Generating data for January 
Generating data for February 
Generating data for March 
Generating data for April 
Generating data for May 
Generating data for June 
Generating data for July 
Generating data for August 
Generating data for September 
Generating data for October 
Generating data for November 
Generating data for December 

m_prec <- generateMonthlyPrec(clim$prec, clim$prec_variation, err_sdev=0.1)

Here we can visualize this temperature monthly variation (minimum, maximum and average) represented by the mean of all the cells simulated.

#Extract each type of temperature data (min, max and average) and obtain the mean value for all the simulated area
tmin <- colMeans(m_temp$tmin, na.rm=TRUE)
tavg <- colMeans(m_temp$tavg, na.rm=TRUE)
tmax <- colMeans(m_temp$tmax, na.rm=TRUE)

#Visualize the average temperature
plot(1:12, tavg, type="l", col="black", lwd=2,
     ylim=range(c(tmin, tmax)),
     xlab="Month", ylab="Temperature")

#Add min and max temperature values
lines(1:12, tmin, col="blue", lty=2)
lines(1:12, tmax, col="red", lty=2)

#Add the legend
legend("topright",
       legend=c("Tmin","Tavg","Tmax"),
       col=c("blue","black","red"),
       lty=c(2,1,2))

And here we show the monthly variation in accumulated precipitation. It can be observed that the months with higher temperatures coincide with lower precipitation levels. This pattern reflects the characteristics of a Mediterranean ecosystem.

#Obtain the accumulated value
s_prec <- colSums(m_prec, na.rm=TRUE)

#Visualize the average temperature
plot(1:12, s_prec, type="l", col="black", lwd=2,
     ylim=range(c(min(s_prec), max(s_prec))),
     xlab="Month", ylab="Precipitation")

#Add the legend
legend("bottomright",
       legend="Precipitation",
       col="black",
       lty=1)

All of this information, in addition to being useful in itself, can also be used to create other types of covariates. To generate the standard bioclimatic variables available in open-source databases, we need the monthly values for cumulative precipitation, along with the maximum and minimum temperatures for each month. To derive these variables from the data, we use the biovars() function from the dismo package.

bv <- biovars(prec=m_prec, tmin=m_temp$tmin, tmax=m_temp$tmax)
bv_rst <- rast(cbind(clim[,1:2], bv), type="xyz")
plot(bv_rst, col= cols)

2.2 Non-climate variables

But if we want to create other types of environmental (non-climate) variables, we can use another function, generateVarb(). Just as we did with the climate variables, here we can also assign the mean value Vmean, its range of variation Vsd, and the range of spatial autocorrelation Vrange. It can be used like any other type of variable (e.g., human footprint, vegetation, percentage of pollutants in the soil, or distribution of another species). For this example we consider it as vegetation.

vegt <- generateVar(nrow=200, ncol=200, Vmean=8, Vsd=3, Vrange=10)

[using unconditional Gaussian simulation]

vegt_rst <- rast(vegt, type="xyz")
plot(vegt_rst, col= cols)

2.3 Edit environmental variables

Once the layers have been created, we can modify certain aspects. For example, we can create a reference layer with the desired attributes and use it to modify those attributes in the layers we’re interested in. In this case, we’re going to change the extent.

#Create the reference layer
r_base.t <- rast(ncol=200, nrow=200, xmin=-0, xmax=200, ymin=0, ymax=200)

#Change the extension of the predictors created above
bv_rst@pntr[["extent"]] <- r_base.t@pntr[["extent"]]
clim_rst@pntr[["extent"]] <- r_base.t@pntr[["extent"]]
vegt_rst@pntr[["extent"]] <- r_base.t@pntr[["extent"]]

Another thing we can do is change the resolution of the layers. In this case, we can reduce the resolution using the aggregate() function from the Terra package.

#climate variables
bv_40 <- aggregate(bv_rst, fact = 5, fun = mean)  #from raster 200x200 units to raster 40x40 units
bv_20 <- aggregate(bv_rst, fact = 10, fun = mean) #from raster 200x200 units to raster 20x20 units
plot(bv_rst, col= cols)

plot(bv_40, col= cols)

plot(bv_20, col= cols)

#vegetation variable
vegt_40 <- aggregate(vegt_rst, fact = 5, fun = mean)  #from raster 200x200 units to raster 40x40 units
vegt_20 <- aggregate(vegt_rst, fact = 10, fun = mean) #from raster 200x200 units to raster 20x20 units
plot(vegt_rst, col= cols)

plot(vegt_40, col= cols)

plot(vegt_20, col= cols)

3. Define the target species’ habitat suitability (or expected abundance)

In this section, we are going to simulate a layer that indicates the suitability of each cell for a target species, based on previously defined predictor variables (such as those we just created in the previous section).

To do this, we start by standardizing all the predictor variables to ensure that each one contributes equally to the model, regardless of the variable’s range and scale.

#vegetation variables
vegt_std <- scale(vegt_rst)
vegt40_std <- scale(vegt_40)  
vegt20_std <- scale(vegt_20)  

#climate variables
clim_std <- scale(clim_rst)

#bioclim variables (what we usually have available when collect data)
bv_std <- scale(bv_rst)
bv40_std <- scale(bv_40)
bv20_std <- scale(bv_20)

To develop this simulation, we want to account for a certain margin of error. That is why we created an error layer with a normal distribution.

r_error <- bv_std$bio1
r_error[] <- rnorm(ncell(r_error), mean = 0, sd = 0.1)
plot(r_error, col= cols)

Once we have all the necessary components, we proceed to simulate the suitability layer for our target species (LAMBDA, which will show the spatial distribution pattern of our species). To do this, we use a Generalized Linear Model (GLM). For this example, we have considered the temperature, precipitation, and vegetation predictor layers as independent variables, but we could also use the bioclimatic layers we have created, other layers we might want to create, or even existing layers that have been adjusted with the same attributes. It is also important to note that when simulating this LAMBDA, different weights can be applied to each variable, as demonstrated in another vignette available on GitHub, Ferrer-Ferrando et. al., 2025, Ecography example.

LAM <- 0.5*clim_std$temp + 0.5*clim_std$prec + 0.5*vegt_std + r_error

plot(LAM, col= cols)

hist(LAM, main="LAMBDA")

However, we see that this result shows both negative and positive values for LAMBDA, and follows a normal distribution (many intermediate values and few values at the extremes). For this reason, we have developed the logistic() function, which performs a logistic transformation on this distribution. With this transformation, we are able to address both issues: 1) We can define the range of values that LAMBDA takes, starting at zero and setting the maximum value we want (L; in this example, we have set it to 1); 2) And we can better highlight the optimal areas from the less optimal ones, obtaining more values at the extremes and fewer in the center (in this example, considering a scenario with many zeros and low values, and the rest more uniform). This second part can be adjusted using two components: the point of inflection x0, and the softness of the curve k.

LAM_lgt <- logistic(LAM, x0 = 0.5, L = 1, k = 2) 

par(mfrow = c(1,1))
plot(LAM[], LAM_lgt[]) #visualize the relationship between original and transformed data

hist(LAM_lgt, main="LAMBDA")

The point of inflection x0 refers to the value within the range where the curve of the logistic function changes. With this transformation, we can also obtain species with higher or lower range of suitable conditions. For example, small range may represent specialist species (with narrow suitable conditions), whereas wide range may represent generalist species (with broader suitable conditions).

#High tolerance example
LAM_lgt_high <- logistic(LAM, -1.5, 1, 2) 

par(mfrow = c(1,1))
plot(LAM[], LAM_lgt_high[], main="Generalist species") 

#Low tolerance example
LAM_lgt_low <- logistic(LAM, 2.5, 1, 2) 

par(mfrow = c(1,1))
plot(LAM[], LAM_lgt_low[], main="Specialist species") 

#Map visualization of "normal", "specialist", and "generalist" species
par(mfrow = c(2,2))
plot(LAM_lgt, main="Normal Species", col= cols) 
plot(LAM_lgt_high, main="Generalist Species", col= cols) 
plot(LAM_lgt_low, main="Specialist Species", col= cols)

rm(LAM_lgt_high)
rm(LAM_lgt_low)

The softness of the curve k indicates the degree of variation along the point of inflection, with the curve being either more abrupt or smoother. This can also refer to the species’ tolerance, where a more gradual change indicates a species that is more tolerant of unfavorable conditions, while an abrupt change indicates a species that is less tolerant of unfavorable conditions.

#High tolerance example
LAM_lgt_high <- logistic(LAM, 0.5, 1, 1) 

par(mfrow = c(1,1))
plot(LAM[], LAM_lgt_high[], main="High tolerance Species") 

#Low tolerance example
LAM_lgt_low <- logistic(LAM, 0.5, 1, 20) 

par(mfrow = c(1,1))
plot(LAM[], LAM_lgt_low[], main="Low tolerance Species") 

#Map visualization of "normal", "specialist", and "generalist" species
par(mfrow = c(2,2))
plot(LAM_lgt, main="Intermediate Species", col= cols) 
plot(LAM_lgt_high, main="High tolerance Species", col= cols) 
plot(LAM_lgt_low, main="Low tolerance Species", col= cols)

rm(LAM_lgt_high)
rm(LAM_lgt_low)

However, even after the logistic transformation, we still do not have exact zero values (they are very low, but not zero). To address this, we need to round the values to obtain them.

LAM_lgt_R <- round(LAM_lgt, 2)

plot(LAM[], LAM_lgt_R[], pch = 16, cex = 0.5) #visualize the relationship between original and transformed data rounded

To conclude this section, here is a comparison of the maps showing the original LAMBDA value and both transformations (logistic and rounding).

par(mfrow = c(2,2))
plot(LAM, main="Normal LAMBDA", col= cols) 
plot(LAM_lgt, main="Logistic LAMBDA", col= cols) 
plot(LAM_lgt_R, main="Logistic LAMBDA Rounded", col= cols)

4. Generate the species’ distribution

In this section, our goal is to generate the species’ distribution in point format, representing individuals of the species in our study area.

4.1 Generate individual presences

Using the generatePoints() function, we can distribute the individuals based on the suitability layer we obtained earlier, following a Poisson Point Process. With this function, you can specify the number of animals you want to generate npoints (as long as the layer contains values between 0 and 1).

#Create the dots (individuals in space)
occurs <- generatePoints(LAM_lgt_R, npoints = 40000)
occurs

Planar point pattern: 40011 points
window: rectangle = [0, 200] x [0, 200] units

#Visualize the individuals on the maps
par(mfrow = c(1,1))
plot(LAM_lgt_R, col= cols)
points(occurs, cex = 0.3, pch=19, col='red')

4.2 Bias the distribution of individuals

If we want to simulate a bias in the distribution of the individuals on the spatial area (because lack of equilibrium with the environment), or a bias in the sampling process (like in real cases occur sometimes), we can do it removing some points (randomly, based on abundance values or by other causes or covariates).

First, we start choosing the number of individuals that we are going to keep based on percentage

xy <- data.frame(cbind(occurs$x, occurs$y)) #extract the coordinates of each dot (individual)
occurssp <- SpatialPoints(xy) #transform to spatial object
occurssp@bbox #see the extention of the dots

            min      max
X1 0.0006515658 199.9976
X2 0.0002552627 199.9752

n <- occurs$n*0.3 #determine the number of dots that we want to keep. 
n <- round(n, 0) #round to a complete number (not decimal)

Now, using the percentage indicated above, we choose this percentage from the total number of presences based on probability (could be based on spatial covariates or coordinates axes, between others)

4.2.1 Bias based on coordinate axis

In this case we determine that the x axis determine the probability of presence of the individuals, where lower values on the x axis determine lower probability of presence of individuals.

#Select from the total of individuals distributed in the space, the percentage indicated (size) based on the x axis probability (prob)
sampled <- sample(occurssp, size = n, prob = occurs$x) 
xy <- data.frame(sampled@coords); names(xy) <- c("x", "y")
samp_ppp <- ppp(xy$x, xy$y, window = occurs$window)

#Visualize the difference between all the individuals and the percentage selected
plot(LAM_lgt_R, main="Total individuals distributed", col= cols)
points(occurs, cex = 0.3, pch=19, col='red')

plot(LAM_lgt_R, main="Percentage selected based on x axis", col= cols)
points(samp_ppp, cex = 0.3, pch=19, col='blue')

4.2.2 Bias based on spatial covariate

In this other case, we determine that a new spatial covariate, like proximity to cities, could determine detect more individuals, so more proximity to cities determine higher probability of detect presence of individuals.

#Determine some cities on the simulated area
set.seed(1234)
city_dots <- randomSamp(n= 8, xrange = c(0,200), yrange = c(0,200), 5) 
v.city_dots <- vect(city_dots, geom = c("x", "y"))

#Visualize the cities on the simulated area
plot(LAM_lgt_R, main="Cities locations", col= cols)
points(v.city_dots, col='red', pch=16)

#Calculate the distance of each cell to the closest city
r.city_dots <- rasterize(v.city_dots, LAM_lgt_R, background = 0) 
d.city_dots <- distance(r.city_dots, target = 0, exclude = NULL, method="haversine") 

Warning: [distance] unknown CRS. Results can be wrong

plot(d.city_dots, main="Distance to the closest city", col= cols)

#Calculate the inverse of the distance to obtain the proximity
interval_dist <- minmax(d.city_dots)
prox.city <- (d.city_dots*(-1)) + interval_dist[2,] 
plot(prox.city, main="Proximity to the closest city", col= cols)

#Calculate the probability of each individual to be selected based on the proximity to the city
v.occurssp <- vect(occurssp)
locs.prox_city <- extract(prox.city, v.occurssp, method = "simple")

#Select from the total of individuals distributed in the space, the percentage indicated (size) based on this proximity to the city (prob)
sampled <- sample(occurssp, size = n, prob = locs.prox_city$last) 
xy <- data.frame(sampled@coords); names(xy) <- c("x", "y")
samp_ppp <- ppp(xy$x, xy$y, window = occurs$window)

#Visualize the difference between all the individuals and the percentage selected
plot(LAM_lgt_R, main="Total individuals distributed", col= cols)
points(occurs, cex = 0.3, pch=19, col='red')

plot(LAM_lgt_R, main="Percentage selected based on city prox", col= cols)
points(samp_ppp, cex = 0.3, pch=19, col='blue')

5. Sample the species’ distribution on the study area

In this section, we will determine the sampling points and prepare the databases we are going to use for subsequent analyses. To do this, we start by assigning the points in our study area where samples are to be taken.

5.1 Determine the grid for sampling

To do this, we have the function randomSamp() that allows us to define the area to be sampled xrange & yrange, the number of samples n and the resolution cellsize (minimum distance between sample points). From this grid, the function will randomly select the location of samples (based on the number of samples that we have specified) from all available points.

#In this case we select a total of 400 locations, with a grid of 200 x 200 units and with a distance between points of 5 units (cellsize)
sam_units <- randomSamp(n= 400, xrange = c(0,200), yrange = c(0,200), cellsize = 5)

#Visualize the locations selected
plot(LAM_lgt_R, main= "Locations selected", col= cols)
points(sam_units, col='red', pch=16, cex=0.5)

#This is optional, but we can obtain the grid structure used to select the sampled points
  #Create our grid employed on the loc sampling selection
sam_grid <- uniformSamp(xrange = c(0,200), yrange = c(0,200), cellsize = 5)

#Visualize the grid and the sampling units selected
plot(LAM_lgt_R, main= "Random locations selected", col= cols)
points(sam_grid, col='blue', pch=16, cex=0.3)
points(sam_units, col='red', pch=16, cex=0.5)

#Add the legend
legend("topright", inset = c(-0.16, 0.35), 
       legend = c("Available points", "Sampled points"), 
       col = c("blue", "red"), pch = 16, pt.cex = c(0.7, 1),
       bty = "n", xpd = TRUE)

Although sampling can sometimes be uniform, such as with a camera trap grid. That is why we have the uniformSamp() function, which allows us to set this uniform sampling by specifying the area to be sampled xrange & yrange, and the resolution cellsize.

#Create our grid employed on the loc sampling selection
uniform_grid <- uniformSamp(xrange = c(0,200), yrange = c(0,200), cellsize = 40)

#Visualize the grid and the sampling units selected
plot(LAM_lgt_R, main= "Uniform locations selected", col= cols)
points(uniform_grid, col='red', pch=16, cex=0.5)

5.2 Data collection (Abundance values)

Once we have determined the sampling points, we move on to collecting the data itself. In this case, data collection can be done in two different ways: as fixed points or as transects.

5.2.1 Point Sampling Method

For this type of sampling, we have the pointSample() function, which allows us to collect information from each unit as a fixed point (for example, as if it were a camera trap or fixed observation points). With this function, you must specify the points you want to sample sites, the location of the individuals to be sampled occ, and the detection radius radius. The function saves a matrix with the individuals as rows and the sampling units as columns, assigning the Euclidean distance from each individual to each sampling unit if it is less than the detection range.

#Pass the occurrences to dataframe
occus <- data.frame(occurs$x, occurs$y)

#Obtain the individual locations (with euclidean distance) for all the sampled points used (sam_units), with a detection range of 2.5 units
point_samp <- pointSample(sites = sam_units, occ = occus, radius = 2.5)

Once we have these matrix with the individuals detected per sampled point and their distances, we save all this information in a data frame with the occurrence ID (the individual detected on the sampling point), the survey ID (the sampling point unit), and the distance of the detection. This type of data can be used for distance sampling models or hierarchical models, and the data could be weighted based on this distance value.

# Get row (occurrence) and column (site) indices with non-NA distances
dist_mat <- point_samp$distance
idx <- which(!is.na(dist_mat), arr.ind = TRUE)

clear_df <- data.frame(
  occurrences = idx[, 1],        # row index = occurrence ID
  sample_id   = idx[, 2],        # col index = site ID
  distance    = dist_mat[idx]    # distance value
)

head(clear_df, 10)

However, if you want to directly access a data frame containing the number of animals detected per sampling unit (abundance values), you can use the countOcc() function, which, by feeding the output directly from the pointSample() function, allows you to obtain this database.

pABN_samp <- countOcc(dmat = point_samp$distance)
head(pABN_samp)

To better understand what kind of data we are collecting, we can visualize the sampling using a map of the entire study area. Here we can see the sampled points with the detection area (dashed line circle) and the individuals detected with a color pattern based on the euclidean distance from the point.

plot(occurs, cex=0.2, asp=1)
points(sam_units, cex=2.5, pch=21, col='red', bg = NULL)

detrange <- 2.5 #determine the detection range of the transect

for (i in 1:400) {
  d <- point_samp$distance[,i]
  msk <- !is.na(d)
  col <- as.integer(cut(d[msk], seq(0, detrange, length.out=5)))
  points(occurs[msk,], pch=16, cex=0.6, col=rainbow(10)[col]) 
}

To visualize better, we make a new figure with a zoom on a specific region of the study area.

ampl <- terra::ext(50, 100, 75, 125)
zoom(LAM_lgt_R, ampl, col= cols)
points(occurs, cex=0.4, asp=1, col = "blue")

points(sam_units, cex=5, pch=21, col='red', bg = NULL)
for (i in 1:400) {
  d <- point_samp$distance[,i]
  msk <- !is.na(d)
  col <- as.integer(cut(d[msk], seq(0, detrange, length.out=5)))
  points(occurs[msk,], pch=16, cex=0.6, col=rainbow(10)[col]) 
}

5.2.2 Transect Sampling Method

For transect sampling methods, we have the transectSample() function, which allows us to collect information from each unit as a continuous transect (for example, as if conducting nighttime censuses using spotlights or tracking clues). With this function, you must specify the points you want to sample sites, the location of the individuals to be sampled obs, the legth of the transect tlen, and the detection range drng. In addition, we can specify the angle of the transect; in this example we used the default option angle = 'random' (random angle per transect), but we can set the angle for each transect using a vector with same number of entries as there are in sites. The function saves a matrix with the individuals as rows and the sampling units as columns, assigning the Euclidean distance from each individual to each sampling unit (the orthogonal distance from the transect) if it is less than the detection range.

#pass the occurrences to dataframe
occus <- data.frame(occurs$x, occurs$y)

ns <- 200 #determine the number of samples
sampling <- randomSamp(ns, c(0,200), c(0,200), 10) #choose these samples inside a grid

tlength <- 5 #determine the legth of the transect
detrange <- 2.5 #determine the detection range of the transect

#Obtain the individual locations (with euclidean distance) for all the sampled transects used (sam_units)
transect_samp <- transectSample(sites = sampling, occ = occus, tlen = tlength, drng = detrange)

Once we have these matrix with the individuals detected per sampled point and their distances, we can do the same thing as we did with fixed points, save all this information in a data frame with the occurrence ID, the survey ID, and the distance of the detection. As well as being used for distance sampling models or hierarchical models.

# Get row (occurrence) and column (site) indices with non-NA distances
dist_mat <- transect_samp$distance
idx <- which(!is.na(dist_mat), arr.ind = TRUE)

clear_df <- data.frame(
  occurrences = idx[, 1],        # row index = occurrence ID
  sample_id   = idx[, 2],        # col index = site ID
  distance    = dist_mat[idx]    # distance value
)

head(clear_df, 10)

We can also obtain the data frame containing the number of animals detected per sampling unit for this type of sampling procedure. By considering the countOcc() function, to which, in this case, the result of the transectSample() function is passed directly.

tABN_samp <- countOcc(dmat = transect_samp$distance)
head(tABN_samp)

Here we present the figure to better understand what kind of data we are collecting. First, the entire study area, with the sampled transects (black lines), surrounded by the detection area (dashed line squares) and the individuals detected with a color pattern based on the euclidean distance from the orthogonal distance from the transect.

plot(occurs, cex=0.2, asp=1)
points(sampling, cex=0.5, pch=16, col='red')

for (i in 1:ns) {
  pnts <- transect_samp$transect
  lines(pnts$x[i, c(2,5)], pnts$y[i, c(2,5)])
  polygon(pnts$x[i,], pnts$y[i,], lty=2)
  
  d <- transect_samp$distance[,i]
  msk <- !is.na(d)
  col <- as.integer(cut(d[msk], seq(0, detrange, length.out=5)))
  points(occurs[msk,], pch=16, cex=0.6, col=rainbow(10)[col]) 
}

Secondly, the figure with zoom, to visualize better.

ampl <- terra::ext(50, 100, 50, 100)
zoom(LAM_lgt_R, ampl, col= cols)
points(occurs, cex=0.4, asp=1, col = "blue")
for (i in 1:ns) {
  pnts <- transect_samp$transect
  lines(pnts$x[i, c(2,5)], pnts$y[i, c(2,5)])
  polygon(pnts$x[i,], pnts$y[i,], lty=2)
  
  d <- transect_samp$distance[,i]
  msk <- !is.na(d)
  col <- as.integer(cut(d[msk], seq(0, detrange, length.out=5)))
  points(occurs[msk,], pch=16, cex=0.6, col=rainbow(10)[col]) 
}

5.3 Transform abundance into presence/absence data

Not always we are interested in obtaining counts, we could prefer or need presence/absence data, or detections considering imperfect detections. In this section, we present the function abundance2pres(). With this function it is possible to transform the abundance values obtained in the previous sections to detections (detected/not-detected binomial value per sampled unit). The functions takes into account the imperfect detection per individual, so higher the number of individuals per sampled unit, higher the probability that detect at least one individual on the sampled unit. For this example, we will use the abundance values obtained on the 5.2.1 section.

abun <- pABN_samp #Abundance values dataset used as initial dataframe
ipd <- 0.5 #Probability of individual detection

sample_prove <- abundance2pres(abun, ipd)

#see the number of sampling units without and with detections 
#(0 = no-detections, 1 = detections)
table(sample_prove$occurrences)

  0   1 
 11 389 

Once we have these transformation done per sampled unit, we save this information in a data frame with the coordinates, the survey ID (the sampled unit), the abundance value (count of individual presents on the sampled unit), and the occurrence values (binomial value of detection or not-detection of the species on the sampled unit).

sampling_units <- data.frame(cbind(x = sam_units$x,
                                   y = sam_units$y,
                                   sampleID = sample_prove$sampleID,
                                   abundance = sample_prove$abundance, 
                                   occurrences = sample_prove$occurrences))
head(sampling_units, 10)

Now, We can visualize each type of data collected showing the total area of study. To compare easily the differences between them.

First, we visualize the abundance values.

# Configuration of the plotting window
par(mfrow = c(1,1))  # Adjust margins
plot(LAM_lgt_R, 
     main = "Spatial Distribution of Abundance", 
     xlab = "Longitude", 
     ylab = "Latitude", 
     legend = FALSE,
     col = cols)  # Background raster

# Red points for real locations
points(occurs, col = "red", pch = 16, cex = 0.05)

# Create color palette (white to dark blue)
pal <- colorRampPalette(c("white", "darkblue"))

# Adjust point sizes based on abundance
point_sizes <- sqrt(sampling_units$abundance + 1) * 0.2 # Scale size proportionally

# Assign colors based on abundance values
points(sampling_units, 
       col = "black", 
       bg = pal(100)[cut(sampling_units$abundance, breaks = 100, labels = FALSE)], 
       cex = point_sizes, 
       pch = 21)

# Add color bar for abundance
colorTable <- designer.colors(100, c("white", "darkblue"))
image.plot(add = TRUE, legend.only = TRUE, 
           zlim = range(sampling_units$abundance, na.rm = TRUE),
           col = colorTable, 
           legend.args = list(text = "Abundance values", 
                              col = "black", cex = 1.5, side = 4, line = 2))

Next, we visualize the detection, non-detection values.

par(mar = c(5, 5, 4, 6))
plot(LAM_lgt_R,
     xlab = "Longitude", 
     ylab = "Latitude",
     legend = FALSE,
     col = cols)

points(occurs, col = "red", pch = 16, cex = 0.05)

# Define specific colors for categorical values (detected not-detected)
category_colors <- c("blue", "green")  # blue for 0, green for 1
names(category_colors) <- c(0, 1)  # Associate colors with occurrences values

# Assign colors based on occurrence values
points(sampling_units, col = "black", 
       bg = category_colors[as.character(sampling_units$occurrences)], 
       pch = 21)

# Add the legend with categoric values
legend("topright", inset = c(-0.05, 0.35), 
       legend = c("0", "1"), 
       col = "black", pt.bg = category_colors, 
       cex = 1.5, bty = "n", pch = 21, xpd = TRUE)

mtext('Detection Values', cex = 1.5, side = 4, line = 0.5)

6. Export the simulated data

In this section, we will show you how to save and export the various layers and data sets that have been created in this vignette. But before running the following code, make sure where is your working directory using here(). If you want to change, modify with setwd().

6.1 Save Sampled data

We start saving the data frame with all the information obtained in the samples (coordinates, occurrence values (detections), abundance values, sample ID). For this data frame as well as the other elements we have been creating, we will create folders to organize them. All of this is made easier by the file.path(), dir.exist(), and dir.create() functions from the base package.

#Create a folder for this layer
savepath <- file.path(here("Data_created"))
if (!dir.exists(savepath)) {
  dir.create(savepath, recursive=TRUE)
  print(paste("Creating directory", savepath))
} else {
  print("Directory already exists. Overwriting.")
}

[1] "Directory already exists. Overwriting."

write.table(sampling_units, here('Data_created/sampling_units.csv'), row.names = FALSE, sep = ",")

But could be interesting save a database with only detections, because this type of format is used in some methodologies such as Maxent.

final_data <- sampling_units[sampling_units$occurrences == 1,]
final_data <- data.frame(final_data[,c(1,2)])
h <- data.frame(ID = cbind(c(rep.int("Virtual_Species", nrow(final_data)))))
final_data <- list(ID=h, x=final_data$x, y=final_data$y)
write.table(final_data, here("Data_created/occ_virtualspecies.csv"), row.names = FALSE, sep = ",")

6.2 Save Suitability layer

We continue saving the suitability layer, which has been used to determine the distribution of the individuals.

terra::writeRaster(LAM_lgt_R, here("Data_created/Expected_abundance.asc"), NAflag = -9999, overwrite=TRUE)

6.3 Save BIOS layers for the analysis

Previously we have standardized the BIOS layers, now we save all of them changed to a range with positive values (keeping all the values >0), avoiding negative and zero values, that could produce problems in softwares like Maxent.

bv_post <-  bv_std + abs(min(bv_std)) + 1

#Create a folder for this layers
savepath <- file.path(here("Data_created/BIOS"))
if (!dir.exists(savepath)) {
  dir.create(savepath, recursive=TRUE)
  print(paste("Creating directory", savepath))
} else {
  print(paste("Directory already exists. Overwriting."))
}

[1] "Directory already exists. Overwriting."

#Save all the BIO rasters with format ASCII
for (i in 1:19) {
  route_rst<-paste("Data_created/BIOS/BIO",i,"_post.asc", sep="")
  writeRaster(bv_post[[i]], here(route_rst), overwrite= TRUE)
}

And we save the BIOS layers at lower resolution, with the same process.

We save the layers at 20x20 units resolution.

bv_post <-  bv20_std + abs(min(bv20_std)) + 1

#Create a folder for this layers
savepath <- file.path(here("Data_created/20BIOS"))
if (!dir.exists(savepath)) {
  dir.create(savepath, recursive=TRUE)
  print(paste("Creating directory", savepath))
} else {
  print(paste("Directory already exists. Overwriting."))
}

[1] "Directory already exists. Overwriting."

#Save all the BIO rasters with format ASCII
for (i in 1:19) {
  route_rst<-paste("Data_created/20BIOS/20BIO",i,"_post.asc", sep="")
  writeRaster(bv_post[[i]], here(route_rst), overwrite= TRUE)
}

And continue with the layers at 40x40 units resolution.

bv_post <-  bv40_std + abs(min(bv40_std)) + 1

#Create a folder for this layers
savepath <- file.path(here("Data_created/40BIOS"))
if (!dir.exists(savepath)) {
  dir.create(savepath, recursive=TRUE)
  print(paste("Creating directory", savepath))
} else {
  print(paste("Directory already exists. Overwriting."))
}

[1] "Directory already exists. Overwriting."

#Save all the BIO rasters with format ASCII
for (i in 1:19) {
  route_rst<-paste("Data_created/40BIOS/40BIO",i,"_post.asc", sep="")
  writeRaster(bv_post[[i]], here(route_rst), overwrite= TRUE)
}

6.4 Save Vegetation layers

We do the same as in the previous section with the vegetation raster layers.

Original high resolution layer.

#Create a folder for this layer
savepath <- file.path(here("Data_created/VEG"))
if (!dir.exists(savepath)) {
  dir.create(savepath, recursive=TRUE)
  print(paste("Creating directory", savepath))
} else {
  print(paste("Directory already exists. Overwriting."))
}

[1] "Directory already exists. Overwriting."

#Save the VEG raster with format ASCII
writeRaster(vegt_std, here("Data_created/VEG/VEGETATION.asc"), overwrite=TRUE)

Layer at 20x20 units resolution.

#Create a folder for this layer
savepath <- file.path(here("Data_created/20VEG"))
if (!dir.exists(savepath)) {
  dir.create(savepath, recursive=TRUE)
  print(paste("Creating directory", savepath))
} else {
  print(paste("Directory already exists. Overwriting."))
}

[1] "Directory already exists. Overwriting."

#Save the VEG raster with format ASCII
writeRaster(vegt20_std, here("Data_created/20VEG/20VEGETATION.asc"), overwrite=TRUE)

Layer at 40x40 units resolution.

#Create a folder for this layer
savepath <- file.path(here("Data_created/40VEG"))
if (!dir.exists(savepath)) {
  dir.create(savepath, recursive=TRUE)
  print(paste("Creating directory", savepath))
} else {
  print(paste("Directory already exists. Overwriting."))
}

[1] "Directory already exists. Overwriting."

#Save the VEG raster with format ASCII
writeRaster(vegt40_std, here("Data_created/40VEG/40VEGETATION.asc"), overwrite=TRUE)

6.5 Save original climate variable (Temperature and Precipitation)

And finally, we do the same with the layers related with climate (temperature and precipitation).

#Create a folder for this layer
savepath <- file.path(here("Data_created/CLIM"))
if (!dir.exists(savepath)) {
  dir.create(savepath, recursive=TRUE)
  print(paste("Creating directory", savepath))
} else {
  print(paste("Directory already exists. Overwriting."))
}

[1] "Directory already exists. Overwriting."

#Save the VEG raster with format ASCII
writeRaster(clim_std$temp, here("Data_created/CLIM/TEMEPRATURE.asc"), overwrite=TRUE)
writeRaster(clim_std$prec, here("Data_created/CLIM/PRECIPITATION.asc"), overwrite=TRUE)

7. Creating a dendrogram. Correlation between covariates

And to conclude this vignette, here is an example of how to select predictor variables. Because it is common that we are interested in different covariates to consider in our models, but these covariates could be correlated, and at high levels this could be a problem. Here we present one way to visualize this correlation and detect clusters of correlation (to choose only one per cluster; visualized as blue boxes). In this case, we use the BIOS that we have created (because we know that they are correlated), to visualize the correlation and choose those that are of interest and with low level of correlation between them. In this case, we consider the limit a correlation value of 0.8 (or a dissimilarity higher than 0.2, employed in the dendogram). All this thanks to as.dendrogram() function of the stats package.

bv_post <-  bv_std + abs(min(bv_std)) + 1
bv_post_df <- as.data.frame(bv_post) #save a data.frame of the variables in columns
bv_post_df <- na.omit(bv_post_df) #eliminate NAs
library(corrplot)

corrplot 0.95 loaded

BIOS_cor <- cor(bv_post_df, use = "pairwise.complete.obs") #create a correlation matrix
BIOS_disim <- 1 - abs(BIOS_cor) #change correlation values to dissimilarity values
BIOS_dist <- as.dist(BIOS_disim) #consider the value as distance
dendrogram <- as.dendrogram(hclust(BIOS_dist, method ="complete")) #create the dendogram
plot(dendrogram) #visualize

#Add the limit of dissimilarity
abline(h = 0.2, col = "red", lwd = 2, lty = 2)
text(x = 1, y = 0.31, labels = "d = 0.2", col = "red", pos = 4)

#Identify the clusters with lower dissimilarity
rect.hclust(hclust(BIOS_dist, method ="complete"), h = 0.2, border = "blue")

Session info

sessionInfo()

R version 4.5.3 (2026-03-11 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 26200)

Matrix products: default
  LAPACK version 3.12.1

locale:
[1] LC_COLLATE=Spanish_Spain.utf8  LC_CTYPE=Spanish_Spain.utf8   
[3] LC_MONETARY=Spanish_Spain.utf8 LC_NUMERIC=C                  
[5] LC_TIME=Spanish_Spain.utf8    

time zone: Europe/Madrid
tzcode source: internal

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

other attached packages:
 [1] corrplot_0.95         here_1.0.2            fields_17.1          
 [4] RColorBrewer_1.1-3    viridisLite_0.4.3     spam_2.11-3          
 [7] dismo_1.3-16          raster_3.6-32         sp_2.2-1             
[10] ecoforge_0.2          knitr_1.51            spatstat.random_3.4-5
[13] spatstat.geom_3.7-3   spatstat.univar_3.1-7 spatstat.data_3.1-9  
[16] terra_1.9-27          gstat_2.1-6          

loaded via a namespace (and not attached):
 [1] Matrix_1.7-4         jsonlite_2.0.0       compiler_4.5.3      
 [4] maps_3.4.3           spatstat.utils_3.2-3 Rcpp_1.1.1-1.1      
 [7] FNN_1.1.4.1          yaml_2.3.12          fastmap_1.2.0       
[10] lattice_0.22-9       deldir_2.0-4         intervals_0.15.5    
[13] htmlwidgets_1.6.4    polyclip_1.10-7      dotCall64_1.2       
[16] rprojroot_2.1.1      rlang_1.2.0          spacetime_1.3-3     
[19] xfun_0.57            otel_0.2.0           cli_3.6.6           
[22] xts_0.14.2           digest_0.6.39        grid_4.5.3          
[25] rstudioapi_0.18.0    vctrs_0.7.2          evaluate_1.0.5      
[28] codetools_0.2-20     zoo_1.8-15           rmarkdown_2.31      
[31] tools_4.5.3          htmltools_0.5.9