knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
To simulate epidemics in a heterogeneous landscape, landsepi needs (among others) these three elements which are
related one each other:
- the spatial coordinates of fields composing the landscape (represented as polygons),
- the allocation of croptypes in the different fields,
- a dispersal matrix for between-field pathogen migration.
landsepi includes built-in landscapes (and associated dispersal matrices for rust pathogens) and an algorithm to allocate croptypes, but is it possible to use your own landscape, dispersal matrix and croptype allocation.
Any landscape can be used to simulate epidemics in landsepi, provided that it is in sp or sf format and contains, at least, polygon coordinates.
library(sf) mylandscape <- st_read(dsn = "myshapefile.shp") library(landsepi) simul_params <- createSimulParams(outputDir = getwd()) simul_params <- setLandscape(simul_params, mylandscape) simul_params@Landscape
Then you can simply call the method allocateLandscapeCroptypes to allocate croptypes to the fields of the landscape with controlled proportions and spatio-temporal aggregation (see tutorial on how to run a simple simulation). Otherwise, you can use your own allocation (see below).
You must define for each year of simulation the index of the croptype ("croptypeID") cultivated in each feature (polygons). Each feature has a field identified by "year_XX" (XX <- seq(1:Nyears+1)) and containing the croptype ID. Note that the allocation must contain one more year than the real number of simulated years (this is only for simulation purpose, the content of the allocation in year Nyears+1 does not affect the result).
| Features/fields | year_1 | year_2 | ... year_Nyears+1 | |---------------- | ------ | ------ | ----------------- | | polygons1 | 13 | 10 | 13 | | polygonsX | 2 | 1 | 2 | | ... | | | |
An example for sf landscape:
mylandscape$year_1 <- c(13,2,4,1,1) # croptypes ID allocated to the different polygons mylandscape$year_2 <- c(2,2,13,1,1)
Then simply add your landscape to the simulation parameters:
simul_params <- setLandscape(simul_params, mylandscape) simul_params@Landscape
To simulate pathogen dispersal, landsepi needs a vectorized matrix giving the probability
of propagule dispersal from any field of the landscape to any other field.
This matrix must be computed before running any simulation with landsepi.
It is a square matrix whose size is the number of fields in the landscape and whose elements are,
for each line $i$ and each column $i'$ the probability $\mu_{ii'}$ that propagules migrate
from field $i$ (whose area is $A_i$) to field $i'$ (whose area is $A_{i'}$). This probability
is computed from:
$$\mu_{ii'} = \frac { \int_{A_i} \int_{A_{i'}} g(\mid\mid z'-z \mid\mid).dz.dz' } { A_i }$$
with $\mid\mid z'-z \mid\mid$ the Euclidian distance between locations $z$ and $z'$ in fields $i$ and $i'$,
respectively, and $g(.)$ the two-dimensional dispersal kernel of the propagules. Note that
$\sum_i \mu_{ii'} = 1$.
landsepi includes built-in dispersal matrices to represent rust dispersal in the five built-in landscapes. These have been computed from a power-law dispersal kernel: $$g(\mid\mid z'-z \mid\mid) = \frac {(b-2).(b-1)} {2.\pi.a^2} . (1+ \frac {\mid\mid z'-z \mid\mid} {a})^{-b}$$ with $a$ the scale parameter and $b$ a parameter related to the width of the dispersal kernel.
A new dispersal matrix must be computed to run simulations with a different landscape or a different dispersal kernel.
The computation of $\mu_{ii'}$ is performed using the CaliFloPP algorithm from the R package RCALI.
The RCALI package has a limited number of built-in dispersal kernels.
However, users can code for their own dispersal kernel.
See section "Details" in the documentation of the function califlopp
to learn how to implement
your own kernel.
Here is an example illustrating the computation of the dispersal matrix on the first landscape supplied in landsepi.
install.packages("RCALI") library(RCALI) library(landsepi) landscape <- landscapeTEST1 Npoly <- length(landscape) Npoly plot(landscape)
For compatibility with the function califlopp
, the landscape can be modified with specific functions
of package sf relative to geographic projection (st_transform
), polygon simplification (st_simplify
).
The function califlopp
needs a specific format for the coordinates of each polygon
(i.e. fields) composing the landscape.
file_land <- "land_rcali.txt" ## input for califlopp file_disp <- "disp_rcali.txt" ## output for califlopp (DO NOT WRITE A PATH) ## Formatting the polygons-file for califlopp cat(Npoly, file=file_land) for (k in 1:Npoly) { ## extract coordinates of polygon vertices coords <- landscape@polygons[[k]]@Polygons[[1]]@coords ## alternatively: # coords <- as.data.frame(landscape$geometry[[k]][[1]]) n <- nrow(coords) cat(NULL, file=file_land, append=T, sep="\n") cat(c(k,k,n), file=file_land, append=T, sep="\t") cat(NULL, file=file_land, append=T, sep="\n") cat(coords[1:n,1], file=file_land, append=T, sep="\t") cat(NULL,file=file_land,append=T,sep="\n") cat(coords[1:n,2], file=file_land, append=T, sep="\t") } cat(NULL, file=file_land, append=T, sep="\n")
Then the function califlopp
calculates the flow of particles between polygons
using an integration method. Here we use the dispersal kernel of oilseed rape pollen
(available in RCALI: use dispf=1
in the arguments of function califlopp
,
see ?califlopp
for details).
param <- list(input=2, output=0, method="cub", dp=6000, dz=6000 , warn.poly=FALSE, warn.conv=FALSE, verbose=FALSE) califlopp(file=file_land, dispf=1, param=param, resfile=file_disp)
The RCALI package has a limited number of built-in dispersal kernels (dispf = 1 in our example).
However, users can code for their own dispersal kernel (let say the name of your kernel is f
)
using dispf=f
in the function califlopp
:
my_df <-function(x, a=40, b=7) ((b-2)*(b-1)/(2*a^2*pi)*(1+(abs(x)/a))^(-b)) param <- list(input=2, output=0, method="cub", dp=6000, dz=6000, warn.poly=FALSE, warn.conv=FALSE, verbose=FALSE) califlopp(file=file_land, dispf=my_df, param=param, resfile=file_disp)
However, if there are many polygons in the landscape, computations may be long.
In this situation, we recommend to replace one of the built-in functions of RCALI by your own function
in the source code, and to recompile RCALI.
See paragraph "The individual dispersion functions" in the details of the documentation of the
califlopp function (?califlopp
).
The output of califlopp must then be reformatted to generate the dispersal matrix that will be further used in landsepi. The vector of field areas can also be generated.
## Import califlopp results disp_df <- getRes(file_disp) ## Double the table because only half of the flows have been calculated emitter <- c(disp_df$poly1, disp_df$poly2) receiver <- c(disp_df$poly2, disp_df$poly1) ## Write a text file containing a vector of areas of all polygons area_e <- c(disp_df$area1, disp_df$area2) area_r <- c(disp_df$area2, disp_df$area1) area <- as.vector(by(area_e, emitter, mean)) write(area, file="area.txt", sep=",") ## Generation of the dispersal matrix name_f <- "mean.flow" flow_mean <- c(disp_df[,name_f], disp_df[,name_f]) flow_f <- cbind(emitter, receiver, flow_mean, area_e, area_r) ## Remove the doublons (i.e. half the lines where emitter == receiver) flow_f[1:nrow(disp_df),][(disp_df$poly2 - disp_df$poly1) == 0,] <- NA flow_f <- flow_f[is.na(apply(flow_f, 1, sum)) == F,] flow_f <- as.data.frame(flow_f) colnames(flow_f) <- c("emitter", "receiver", "flow", "area_e", "area_r") flow_f <- flow_f[order(flow_f$emitter),] ## lines: emitter ## columns: receiver matrix_f <- NULL for(k in 1:Npoly){ ## flow divided by the emitter area matrix_f <- cbind(matrix_f, flow_f$flow[flow_f$receiver==k] / area) } ## Normalisation of the matrix (reflecting boundaries) ## (do not normalise for absorbing boundaries) flowtot_f <- apply(matrix_f,1,sum) for(k in 1:Npoly){ matrix_f[k,] <- (matrix_f[k,] / flowtot_f[k]) ## In order to have sum == 1 } write(as.vector(matrix_f), file="dispersal.txt", sep=",")
Then, to read the file, use:
disp_patho <- scan("dispersal.txt", sep=",")
Landscape structure can be plotted using the basic function plot()
, or using the landsepi function
plotland()
:
landscape <- landscapeTEST1 plot(landscape) plotland(landscape)
To highlight a specific field:
poly <- 10 colFields <- rep("white", length(landscape)) colFields[poly] <- "red" plot(landscape, col = colFields)
To check the dispersal matrix and represent in a graphic the flow emitted by a specific polygon, use:
## convert dispersal in matrix mat <- matrix(disp_patho, nrow=sqrt(length(disp_patho))) poly <- 1 dispToPlot <- log10(mat[poly,] +1E-20) ## 1E-20 to avoid log(0) ## Colour palette nCol <- 11 whiteYellowRed <- colorRampPalette(c("white", "#FFFF99", "#990000")) col_disp <- whiteYellowRed(nCol) intvls <- seq(min(dispToPlot) - 1, max(dispToPlot) + 1, length.out=nCol) intvls_disp <- findInterval(dispToPlot, intvls) ## Plot plot(landscape, col = col_disp[intvls_disp], main=paste("Dispersal from polygon", poly))
With package ggplot2:
library(ggplot2) ggplot(landscape) + ggtitle(paste("Dispersal from polygon", poly)) + geom_sf(colour="black", aes(fill = dispToPlot)) + scale_fill_gradientn(name="Prob. of\ndispersal", colours=rev(heat.colors(10)), breaks=-1:-10, labels=10^(-1:-10)) + # theme_classic() + theme(axis.line=element_blank(),axis.text.x=element_blank(), axis.text.y=element_blank(),axis.ticks=element_blank(), axis.title.x=element_blank(), axis.title.y=element_blank(), panel.background=element_blank(),panel.border=element_blank(),panel.grid.major=element_blank(), panel.grid.minor=element_blank(),plot.background=element_blank())
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.