1. Input data
1.1 Create lattice arena
Initially, we’ll establish the physical space designated for the mosquito introduction.
We’ll construct a squared lattice arena with:
- Dimensions: 5 km x 5 km
- Resolution: 250 m (which gives us 20 columns and 20 rows)
- Total Cells: 400
# Libraries
library(gstat)
library(terra)
library(eesim)
library(dynamAedes)
library(ggplot2)
Sys.setlocale("LC_TIME", "en_GB.UTF-8")# [1] "en_GB.UTF-8"gridDim <- 20 # 5000m/250 m = 20 columns and rows
xy <- expand.grid(x=1:gridDim, y=1:gridDim)After establishing the lattice arena, we’ll introduce a spatial pattern to it. This pattern becomes crucial when we aim to incorporate spatial autocorrelation (SAC) into the temperature time series later on.
To achieve this spatial autocorrelation:
- Semivariogram Model: We’ll employ a semivariogram model with specified parameters:
- Sill: This is the value that the semivariogram reaches at the defined range.
- Range: This represents the distance at which spatial correlation becomes zero.
- Unconditional Gaussian Simulation: Using this approach, we will predict the semivariogram model across our lattice grid, producing the desired spatial autocorrelated pattern.
varioMod <- gstat::vgm(psill=0.005, range=100, model='Exp') # psill = partial sill = (sill-nugget)
# Set up an additional variable from simple kriging
zDummy <- gstat::gstat(formula=z~1,
locations = ~x+y,
dummy=TRUE,
beta=1,
model=varioMod,
nmax=1)
# Generate a randomly autocorrelated predictor data field
set.seed(123)
xyz <- predict(zDummy, newdata=xy, nsim=1)# [using unconditional Gaussian simulation]Having defined our spatial pattern, we can now generate a spatially autocorrelated raster. We’ll achieve this by incorporating the SA variable, denoted as xyz$sim1, into the SpatRaster object.
It’s worth noting that such a spatially autocorrelated surface might, for instance, depict the distribution of vegetation cover within an urban landscape, providing valuable context and depth to our model.
utm32N <- "+proj=utm +zone=32 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
r <- terra::rast(nrow=gridDim, ncol=gridDim, crs=utm32N, ext=terra::ext(1220000,1225000, 5700000,5705000))
terra::values(r) <- xyz$sim1
plot(r, main="SAC landscape")
# convert to a data.frame
df <- data.frame("id"=1:nrow(xyz), terra::crds(r))
bbox <- terra::as.polygons(terra::ext(r), crs=utm32N)
# Store Parameters for autocorrelation
autocorr_factor <- terra::values(r)1.2 Simulate temperature data with seasonal trend
We simulate a 1-year temperature time series with seasonal trend. For the time series we consider a mean value of 16°C and standard deviation of 2°C.
ndays <- 365*1 #length of the time series in days
set.seed(123)
sim_temp <- eesim::create_sims(n_reps = 1,
n = ndays,
central = 16,
sd = 2,
exposure_type = "continuous",
exposure_trend = "cos1", exposure_amp = -1.0,
average_outcome = 12,
outcome_trend = "cos1",
outcome_amp = 0.8,
rr = 1.0055)A visualisation of the distribution of temperature values and temporal trend.
hist(sim_temp[[1]]$x,
xlab="Temperature (°C)",
main="Histogram of simulated temperatures")
plot(sim_temp[[1]]$date,
sim_temp[[1]]$x,
main="Simulated temperatures seasonal trend",
xlab="Date", ylab="Temperature (°C)"
)
Once we have our spatially autocorrelated raster, the next step is to spatially distribute our temperature time series. We achieve this by multiplying the time series with the autocorrelated surface we generated in the previous step.
By doing so, we incorporate spatial variability into the temperature data, ensuring that the temperature values differ across the landscape in a manner consistent with the underlying spatial pattern (e.g., vegetation cover). This offers a more realistic representation of how temperatures might vary across different areas due to various environmental factors.
mat <- do.call(rbind, lapply(1:ncell(r), function(x) {
d_t <- sim_temp[[1]]$x*autocorr_factor[[x]]
return(d_t)
}))oldpar <- par(mfrow = c(1,2)) A comparison between the distribution of the initial temperature time series and autocorrelated temperature surface
par(mfrow=c(2,1))
hist(mat, xlab="Temperature (°C)", main="Histogram of simulated spatial autocorreled temperature")
hist(sim_temp[[1]]$x, xlab="Temperature (°C)", main="Histogram of simulated temperatures", col="red")
par(mfrow=c(1,1))par(oldpar) 1.2.1 Format temperature data
names(mat) <- paste0("d_", 1:ndays)
df_temp <- cbind(df, mat)1.3 Incorporating Medium-Range Passive Dispersal
In many real-world scenarios, mosquitoes do not just disperse naturally – they often hitch rides and spread to new areas, especially through human transportation methods. In the context of our model, this concept is termed as “medium-range passive dispersal”.
To account for this mechanism of dispersal, we introduce an arbitrary road segment into our spatial domain. This road segment will represent pathways along which adult mosquitoes can disperse passively, most notably through car traffic.
set.seed(123)
roads <- vect(terra::crds(terra::spatSample(bbox, 5, method="random")),
type="lines",
crs=utm32N)
# Check simulated segment
plot(r)
terra::plot(roads, add=TRUE)
After defining the road segment we add a “buffer” of 250 m around the road segment. Adult mosquitoes that reach or develop into cells comprised in the 250 buffer around roads are thus able to undergo passive dispersal.
# Create a buffer around the converted polygon
buff <- terra::buffer(roads, width=250)
# Check grid, road segment and buffer
plot(r)
plot(buff, add=T)
plot(roads, add=T, col="red")
Next, we derive a distance matrix between cells comprised in the spatial buffer along the road network. First, we select the cells.
df_sp <- vect(x=df, geom=c("x", "y"), crs=utm32N)
df_sp <- terra::intersect(df_sp, buff)
# Check selected cells
plot(r)
plot(buff, add=T)
plot(df_sp, add=T, col="red")
Then, we compute the Euclidean distance (in meters) between each selected cell.
dist_matrix <- as.matrix(stats::dist(terra::crds(df_sp)))
At local scale, the interpretation of this output is more nuanced than for the other scales: a pixel having psi=0 can be a pixel where all the simulations resulted in an extinction or where the species has not yet arrived through dispersal. To understand the value of the pixels may be useful to iterate the psi_sp function along the vector of simulated days.
Note that if only a small number of mosquitoes are present in a pixel over many iterations, quantiles may be 0 (especially for low quantiles) and you may see a series of empty rasters!