O3-landscape_dispersal.R
In landsepi: Landscape Epidemiology and Evolution

## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----eval=FALSE---------------------------------------------------------------
# library(sf)
# mylandscape <- st_read(dsn = "myshapefile.shp")
# library(landsepi)
# simul_params <- createSimulParams(outputDir = getwd())
# simul_params <- setLandscape(simul_params, mylandscape)
# simul_params@Landscape

## ----eval=FALSE---------------------------------------------------------------
# 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)

## ----eval=FALSE---------------------------------------------------------------
# simul_params <- setLandscape(simul_params, mylandscape)
# simul_params@Landscape

## ----eval=FALSE---------------------------------------------------------------
# install.packages("RCALI")
# library(RCALI)
# library(landsepi)
# landscape <- landscapeTEST1
# Npoly <- length(landscape)
# Npoly
# plot(landscape)

## ----eval=FALSE---------------------------------------------------------------
# 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")

## ----eval=FALSE---------------------------------------------------------------
# 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)

## ----eval=FALSE---------------------------------------------------------------
# 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)

## ----eval=FALSE---------------------------------------------------------------
# ## 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=",")

## ----eval=FALSE---------------------------------------------------------------
# disp_patho <- scan("dispersal.txt", sep=",")

## ----eval=FALSE---------------------------------------------------------------
# landscape <- landscapeTEST1
# plot(landscape)
# plotland(landscape)

## ----eval=FALSE---------------------------------------------------------------
# poly <- 10
# colFields <- rep("white", length(landscape))
# colFields[poly] <- "red"
# plot(landscape, col = colFields)

## ----eval=FALSE---------------------------------------------------------------
# ## 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))

## ----eval=FALSE---------------------------------------------------------------
# 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())