Local scale model
The “local scale” means that the model accounts for both active and passive dispersal of the mosquitoes. With this setting, the model requires three input datasets:
a numerical temperature matrix (in degree Celsius) defined in space and time (space in the rows, time in the columns);
a two-column numerical matrix reporting the coordinates (in meters) of each space-unit (cell);
a numerical distance matrix which reports the distance in meters between the cells connected through a road network.
For the purpose of this tutorial, we will use the following simulated datasets:
- A 5 km lattice grid with 250 m cell size;
- A 1-year long spatially and temporally correlated temperature time series;
- A matrix of distances between cells connected through a simulated road network;
Input data
Create lattice arena
First, we define the physical space into which the introduction of our mosquitoes will happen. We define a squared lattice arena having 5 km side and 250 m resolution (20 columns and 20 rows, 400 total cells).
gridDim <- 20 # 5000m/250 m = 20 columns and rows
xy <- expand.grid(x=1:gridDim, y=1:gridDim)We then add a spatial pattern into the lattice area. This spatial pattern will be used later to add spatial correllation (SAC) to the temperature time series. The spatial autocorrelated pattern will be obtained using a semivariogram model with defined sill (value that the semivarion attains at the range) and range (distance of 0 spatial correlation) and then predicting the semivariogram model over the lattice grid using unconditional Gaussian simulation.
varioMod <- vgm(psill=0.005, range=100, model='Exp') # psill = partial sill = (sill-nugget)
# Set up an additional variable from simple kriging
zDummy <- 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]We generate a spatially autocorrelated raster adding the SA variable (xyz$sim1) to the RasterLayer object. The autocorrelated surface could for example represent the distribution of vegetation cover in a urban landscape.
utm32N <- "+proj=utm +zone=32 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
r <- raster(nrow=gridDim, ncol=gridDim, crs=utm32N, ext=extent(1220000,1225000, 5700000,5705000))
values(r)=xyz$sim1
plot(r, main="SAC landscape")
df <- data.frame("id"=1:nrow(xyz), raster::coordinates(r))
bbox <- as(extent(r), "SpatialPolygons")
projection(bbox) <- projection(utm32N)
# Store Parameters for autocorrelation
autocorr_factor <- values(r)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 <- 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)"
)
We can then “expand onto space” the temperature time series by multiplying it with the autocorrelated surface simulated above.
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) 
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.
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!