Commit eb99ff23 authored by Jean-Francois Rey's avatar Jean-Francois Rey
Browse files

wgpkgc++ to intR

parent c8556abb
......@@ -37,7 +37,7 @@ r-Mac-build_test:
- "Mac"
- "R"
script:
- "Rscript -e \"pkgbuild::compile_dll()\""
- 'TMPDIR="/tmp" && Rscript -e "pkgbuild::compile_dll()"'
- "Rscript -e \"roxygen2::roxygenize('.', roclets=c('rd', 'collate', 'namespace'))\""
- "R CMD build . --resave-data"
- "R CMD check --as-cran $(ls -rt landsepiDev_* |tail -1)"
......
......@@ -11,6 +11,8 @@ export(createSimulParams)
export(demo_landsepi)
export(epid_output)
export(evol_output)
export(getGPKGArea)
export(getGPKGRotation)
export(invlogit)
export(loadCroptypes)
export(loadCultivar)
......
......@@ -47,6 +47,45 @@ createLandscapeGPKG <- function(landscape, outputfile) {
return(outputfile)
}
#' getGPKGArea
#' @description get Area layer as a vector
#' @param gpkgfile a GPKG file
#' @return a vector of the area of each polygons
#' @export
getGPKGArea <- function(gpkgfile) {
area <- st_read(
gpkgfile,
"area"
)
st_geometry(area) <- NULL
areadf <- as.data.frame(area)
return(as.vector(areadf[,1]))
}
#' getGPKGRotation
#' @description get Croptypes ID rotation years as a matrix
#' @param gpkgfile a GPKG file
#' @return a matrix as rows for polygons and cols for years
#' @export
getGPKGRotation <- function(gpkgfile) {
croptypeID <- st_read(
gpkgfile,
"CroptypeID"
)
st_geometry(croptypeID) <- NULL
cdf <- as.data.frame(croptypeID)
ncol <- length(grep("^year_", colnames(cdf)) %in% colnames(cdf))
#nrow <- nrow(cdf)
croptype_matrix <- as.matrix(cdf[,grep("^year_", colnames(cdf))], ncol=ncol)
return(croptype_matrix)
}
#' @title GPKGAddTables
#' @description Add non spatial data table definitions (sqlite) into a GPKG file.
#' It add Cultivar, CultivarList, Gene, GeneList table
......@@ -196,6 +235,22 @@ getGPKGCroptypes <- function(inputGPKG) {
return(dfcroptypes)
}
#' getGPKGCroptypesRaw
#' @description Get Croptypes and cultivars proportions from a GPKG file
#' @param inputGPKG a GPKG filename
#' @return a data.frame with croptypeID, CultivarID, and proportions
getGPKGCroptypesRaw <- function(inputGPKG) {
inputDB <- DBI::dbConnect(RSQLite::SQLite(), inputGPKG)
# Select croptypes and cultivars name
sql <- "SELECT croptypeID, cultivarID, proportion FROM CultivarList;"
croptypes <- DBI::dbGetQuery(inputDB, sql)
DBI::dbDisconnect(inputDB)
return(croptypes)
}
#' getGPKGCultivars
#' @description Get Cultivars from a GPKG file
#' @param inputGPKG a GPKG filename
......
......@@ -198,7 +198,7 @@ addDefaultParameters <- function(params) {
message("no Host dispersal matrix added")
message("loading default identity matrix (absence of host dispersal)")
message("Use methods setDispersalHost or loadDispersalHost to modify value")
disp_host <- loadDispersalHost(type = "no")
disp_host <- loadDispersalHost(params, type = "no")
params <- setDispersalHost(params, disp_host)
}
......@@ -313,7 +313,9 @@ saveDeploymentStrategy <- function(params, overwrite = FALSE) {
}
## create gpkg file
message("Create ", paste0(params@OutputDir, "/", params@OutputGPKG), " file")
gpkgFile <- createLandscapeGPKG(params@Landscape, paste0(params@OutputDir, "/", params@OutputGPKG))
## save only years_ and area ?
land <- params@Landscape[,grep("^year_", colnames(params@Landscape))]
gpkgFile <- createLandscapeGPKG(land, paste0(params@OutputDir, "/", params@OutputGPKG))
## add data tables
GPKGAddTables(gpkgFile)
## Fill data tables
......@@ -340,14 +342,34 @@ runSimul <- function(params, outputs=list(epid_outputs="all", evol_outputs="all"
, GLAnoDis=1.48315, audpc100S=0.38)
, graphic = TRUE, write = TRUE, video = FALSE, keepRawResults = FALSE) {
### !!!!! BE CAREFUL !!!!! ###
### croptypes, cultivars, genes and cultivarsGenes have to be ordored all in the same way
### ID have to match row index and col index
## TODO move this method somewhere else
params <- addDefaultParameters(params)
cultivars_genes_list <- lapply(1:nrow(params@Cultivars), FUN = function(i){
return(which(params@CultivarsGenes[i,] == 1)-1)
})
cdf <- as.data.frame(params@Landscape)
ncol <- length(grep("^year_", colnames(cdf)) %in% colnames(cdf)) ## Utiliser la valeur de time Nyears
rotation <- as.matrix(cdf[,grep("^year_", colnames(cdf))], ncol=ncol)
croptypes_cultivars <- params2CroptypeBDD(params)[,c(2,3,4)]
## Run the simulation
simul_landsepi(
outputs <- simul_landsepi(
seed = params@Seed,
time_param = params@TimeParam,
GPKGLandscape = paste0(params@OutputDir, "/", params@OutputGPKG),
croptype_names = params@Croptypes$croptypeName,
cultivars = params@Cultivars,
cultivars_genes_list = cultivars_genes_list,
genes = params@Genes,
landscape = as_Spatial(st_geometry(params@Landscape)),
area = as.vector(params@Landscape$area[,1]),
rotation = rotation,
croptypes_cultivars_prop = croptypes_cultivars,
basic_patho_param = params@Pathogen,
disp_patho = params@DispPatho,
disp_host = params@DispHost,
......@@ -358,6 +380,8 @@ runSimul <- function(params, outputs=list(epid_outputs="all", evol_outputs="all"
video = video,
keepRawResults = keepRawResults
)
#return(outputs)
}
#' setSeed
......@@ -415,6 +439,8 @@ setLandscape <- function(params, land) {
} else {
params@Landscape <- land
}
params@Landscape$area <- data.frame(area = st_area(params@Landscape))
return(params)
}
......@@ -450,7 +476,7 @@ checkLandscape <- function(params) {
ldf <- as.data.frame(land)
# check CroptypeID present in Landscape
if (sum(!unique(unlist(ldf)) %in% params@Croptypes$croptypeID) != 0) {
if (sum(!unique(as.integer(unlist(ldf[,grep("^year_", colnames(ldf))]))) %in% params@Croptypes$croptypeID) != 0) {
warning("Croptypes undef in Landscape")
warning("Croptypes ID : ", params@Croptypes$croptypeID)
warning("Croptypes ID in Landscape : ", unique(unlist(ldf)))
......@@ -631,6 +657,8 @@ allocateLandscapeCroptypes <- function(params, rotation_period, rotation_sequenc
, prop, aggreg, graphic=TRUE) {
### TODO Check validity of params slot before Agriland
orig_landscape = params@Landscape
croptypeSP <- AgriLand(as_Spatial(params@Landscape),
Nyears = params@TimeParam$Nyears,
rotation_period = rotation_period,
......@@ -643,6 +671,11 @@ allocateLandscapeCroptypes <- function(params, rotation_period, rotation_sequenc
outputDir = params@OutputDir
)
params@Landscape <- st_as_sf(croptypeSP)
params@Landscape$area <- data.frame(area = st_area(params@Landscape))
if( length(orig_landscape$ID) != 0 && length(orig_landscape$Name) != 0 ) {
params@Landscape$ID <- orig_landscape$ID
params@Landscape$NAME <- orig_landscape$NAME
}
return(params)
}
......@@ -723,14 +756,14 @@ checkPathogen <- function(params) {
#' loadCroptype
#' @description create a croptype data.frame parameters depending of cultivars composing the croptype
#' @param params a LandsepiParams object
#' @param name the croptype name
#' @param cultivarsInCroptype name of cultivars composing the croptype
#' @param prop vector of proportions of each cultivar in the croptype. Default to balanced proportions.
#' @return a data.frame with croptype parameters
#' @export
# loadCroptype
# @description create a croptype data.frame parameters depending of cultivars composing the croptype
# @param params a LandsepiParams object
# @param name the croptype name
# @param cultivarsInCroptype name of cultivars composing the croptype
# @param prop vector of proportions of each cultivar in the croptype. Default to balanced proportions.
# @return a data.frame with croptype parameters
# @export
# loadCroptype <- function(params, name, cultivarsInCroptype, prop=NULL) {
# cultivar_names <- params@Cultivars$cultivarName
# cropt <- list(croptypeName=name)
......
......@@ -10,11 +10,9 @@
#' \item Nyears = number cropping seasons,
#' \item nTSpY = number of time-steps per cropping season.
#' }
#' @param landscape list of landscape parameters:\itemize{
#' \item shapefilename = name of the landscape,
#' \item layername_croptype = name of the layer containing croptype indices,
#' \item layername_area = name of the layer containing polygon areas.
#' }
#' @param area_vector a Rcpp::NumericVector containing polygons area
#' @param rotation_matrix a Rcpp::NumericMatrix containing polygons (row) croptypes indices by years (col)
#' @param croptypes_cultivars_prop a Rcpp::NumericMatrix (data.frame) with CroptypesID, CultivarsID and Proportion
#' @param dispersal list of dispersal parameters:\itemize{
#' \item disp_patho = vectorised dispersal matrix of the pathogen,
#' \item disp_host = vectorised dispersal matrix of the host.
......@@ -101,12 +99,49 @@
#' Each file indicates for every time-step the number of individuals in each field, and when appropriate for
#' each host and pathogen genotypes).
#'
#' @examples
#' \dontrun{
#' ## Spatially-implicit simulation with 2 patches (S + R) during 3 years
#' cultivars <- as.list(rbind(loadCultivar(name="Susceptible", type="growingHost")
#' , loadCultivar(name="Resistant", type="growingHost")))
#' names(cultivars)[names(cultivars)=="cultivarName"] <- "name"
#' cultivars <- c(cultivars, list(Nhost=2, sigmoid_kappa_host=0.002, sigmoid_sigma_host=1.001, sigmoid_plateau_host=1
#' , cultivars_genes_list=list(numeric(0),0)))
#'
#' model_landsepi(seed=1,
#' time_param = list(Nyears=3, nTSpY=120),
#' basic_patho_param = loadPathogen(disease = "rust"),
#' inits = list(pI0=5E-4),
#' area_vector = c(10000, 10000),
#' dispersal = list(disp_patho=c(0.99,0.01,0.01,0.99), disp_host=c(1,0,0,1)),
#' rotation_matrix = as.matrix(data.frame(year_1=c(0,1), year_2=c(0,1), year_3=c(0,1), year_4=c(0,1))),
#' croptypes_cultivars_prop = as.matrix(data.frame(croptypeID=c(0,1), cultivarID=c(0,1), proportion=c(1,1))),
#' cultivars_param = cultivars,
#' genes_param = as.list(loadGene(name="MG", type="majorGene")))
#'
#'
#' ## 1-year simulation of a rust epidemic in pure susceptible crop in a single patch of 1 ha
#' cultivars <- as.list(rbind(loadCultivar(name="Susceptible", type="growingHost")))
#' names(cultivars)[names(cultivars)=="cultivarName"] <- "name"
#' cultivars <- c(cultivars, list(Nhost=1, sigmoid_kappa_host=0.002, sigmoid_sigma_host=1.001, sigmoid_plateau_host=1
#' , cultivars_genes_list=list(numeric(0))))
#' model_landsepi(seed=1,
#' time_param = list(Nyears=1, nTSpY=120),
#' basic_patho_param = loadPathogen(disease = "rust"),
#' inits = list(pI0=5E-4),
#' area_vector = c(10000),
#' dispersal = list(disp_patho=c(1), disp_host=c(1)),
#' rotation_matrix = as.matrix(data.frame(year_1=c(0), year_2=c(0))),
#' croptypes_cultivars_prop = as.matrix(data.frame(croptypeID=c(0), cultivarID=c(0), proportion=c(1))),
#' cultivars_param = cultivars,
#' genes_param = as.list(loadGene(name="MG", type="majorGene")))
#' }
#' @references Rimbaud L., Papaïx J., Rey J.-F., Barrett L. G. and Thrall P. H. (2018).
#' Assessing the durability andefficiency of landscape-based strategies to deploy plant resistance to pathogens.
#' \emph{PLoS Computational Biology} 14(4):e1006067.
#'
#' @export
model_landsepi <- function(time_param, landscape, dispersal, inits, seed, cultivars_param, basic_patho_param, genes_param) {
invisible(.Call(`_landsepiDev_model_landsepi`, time_param, landscape, dispersal, inits, seed, cultivars_param, basic_patho_param, genes_param))
model_landsepi <- function(time_param, area_vector, rotation_matrix, croptypes_cultivars_prop, dispersal, inits, seed, cultivars_param, basic_patho_param, genes_param) {
invisible(.Call(`_landsepiDev_model_landsepi`, time_param, area_vector, rotation_matrix, croptypes_cultivars_prop, dispersal, inits, seed, cultivars_param, basic_patho_param, genes_param))
}
......@@ -38,11 +38,9 @@
#' \item nTSpY = number of time-steps per cropping season.
#' }
#' @param Npatho number of pathogen genotypes.
#' @param landscape list of landscape parameters:\itemize{
#' \item shapefilename = name of the landscape,
#' \item layername_croptype = name of the layer containing croptype indices,
#' \item layername_area = name of the layer containing polygon areas.
#' }
#' @param area a vector containing polygons areas
#' @param rotation a data.frame containing croptypes indices by years for polygons indices
#' @param croptypes a data.frame containing croptypeID, CulitvarsID and proportion
#' @param cultivars_param list of parameters associated with each host genotype (i.e. cultivars)\itemize{
#' \item name = vector of cultivar names,
#' \item initial_density = vector of host densities (per square meter) at the beginning of the cropping season,
......@@ -99,7 +97,7 @@
#' @importFrom grDevices colorRampPalette dev.off graphics.off gray png tiff
#' @importFrom utils write.table
#' @export
epid_output <- function(types="all", time_param, Npatho, landscape, cultivars_param, eco_param
epid_output <- function(types="all", time_param, Npatho, area, rotation, croptypes, cultivars_param, eco_param
, GLAnoDis=cultivars_param$max_density[1] ## true for non-growing host
, ylim_param = list(audpc=c(0,0.38)
, gla_abs=c(0,1.48)
......@@ -124,11 +122,12 @@ epid_output <- function(types="all", time_param, Npatho, landscape, cultivars_pa
nTS <- Nyears * nTSpY ## Total number of time-steps
## Landscape (Need package sf)
landscape_tmp <- st_read(landscape$shapefilename, layer = landscape$layername_croptype)
landscape_spatial <- as(landscape_tmp, "Spatial")
rotation <- data.frame(landscape_spatial)
croptypes <- st_read(landscape$shapefilename, layer = "CultivarList")
area <- st_read(landscape$shapefilename, layer = landscape$layername_area)$area
#landscape_tmp <- st_read(landscape$shapefilename, layer = landscape$layername_croptype)
#landscape_spatial <- as(landscape_tmp, "Spatial")
#rotation <- data.frame(landscape_spatial)
#croptypes <- st_read(landscape$shapefilename, layer = "CultivarList")
#area <- st_read(landscape$shapefilename, layer = landscape$layername_area)$area
rotation <- data.frame(rotation)
areaTot <- sum(area)
Npoly <- length(area)
......@@ -798,11 +797,15 @@ evol_output <- function(types="all", time_param, Npoly, cultivars_param, genes_p
#' \item nTSpY = number of time-steps per cropping season.
#' }
#' @param Npatho number of pathogen genotypes.
#' @param landscape list of landscape parameters:\itemize{
#' \item shapefilename = name of the landscape,
#' \item layername_croptype = name of the layer containing croptype indices,
#' \item layername_area = name of the layer containing polygon areas.
#' }
#' @param landscape a sp object
#' @param area a vector containing polygons area
#' @param rotation a data.frame containing polygons (row) CroptypesID by years (col)
#' @param croptypes a data.frame containing CroptypeID, CultivarsID and proportion
# @param landscape list of landscape parameters:\itemize{
# \item shapefilename = name of the landscape,
# \item layername_croptype = name of the layer containing croptype indices,
# \item layername_area = name of the layer containing polygon areas.
# }
#' @param croptype_names a vector of croptype names (for legend).
#' @param cultivars_param list of parameters associated with each host genotype (i.e. cultivars)
#' when cultivated in pure crops:\itemize{
......@@ -831,7 +834,7 @@ evol_output <- function(types="all", time_param, Npoly, cultivars_param, genes_p
#' @import parallel
#' @importFrom grDevices colorRampPalette dev.off graphics.off gray png tiff
#' @export
video <- function(audpc, time_param, Npatho, landscape, croptype_names=c(), cultivars_param, audpc100S, keyDates = NULL
video <- function(audpc, time_param, Npatho, landscape, area, rotation, croptypes, croptype_names=c(), cultivars_param, audpc100S, keyDates = NULL
, nMapPY = 5, path = getwd()) {
## TODO: update @include and @importFrom --> what is "gray" ?
......@@ -845,11 +848,12 @@ video <- function(audpc, time_param, Npatho, landscape, croptype_names=c(), cult
nTS <- Nyears * nTSpY ## Total number of time-steps
## Landscape (Need package sf)
landscape_tmp <- st_read(landscape$shapefilename, layer = landscape$layername_croptype)
landscape_spatial <- as(landscape_tmp, "Spatial")
rotation <- data.frame(landscape_spatial)
croptypes <- st_read(landscape$shapefilename, layer = "CultivarList")
area <- st_read(landscape$shapefilename, layer = landscape$layername_area)$area
# landscape_tmp <- st_read(landscape$shapefilename, layer = landscape$layername_croptype)
# landscape_spatial <- as(landscape_tmp, "Spatial")
# rotation <- data.frame(landscape_spatial)
# croptypes <- st_read(landscape$shapefilename, layer = "CultivarList")
# area <- st_read(landscape$shapefilename, layer = landscape$layername_area)$area
rotation <- data.frame(rotation)
Npoly <- length(area)
## croptype parameters (for legend)
......@@ -996,7 +1000,7 @@ video <- function(audpc, time_param, Npatho, landscape, croptype_names=c(), cult
)
## Map dynamic of diseased hosts
plotland(landscape_spatial, col_prev[intvlsIR], density_hatch[rotation[, index_year] + 1], angle_hatch[rotation[, index_year] + 1]
plotland(landscape, col_prev[intvlsIR], density_hatch[rotation[, index_year] + 1], angle_hatch[rotation[, index_year] + 1]
, col_prev, density_hatch, angle_hatch, title, subtitle
, croptype_names
, legend_prev, " ")
......
......@@ -28,8 +28,14 @@
#' \item Nyears = number cropping seasons,
#' \item nTSpY = number of time-steps per cropping season.
#' }
#' @param GPKGLandscape a GPKG file generated by landsepi method saveDeploymentStrategy. Areas must be
#' expressed in square meters.
#' @param croptype_names Croptypes Name
#' @param cultivars a data.frame of cultivars parameters
#' @param cultivars_genes_list Genes list indices by cultivars order
#' @param genes a data.frame of Genes parameters
#' @param landscape a sp object
#' @param area a vector containing polygons area (square meters)
#' @param rotation a data.frame containing croptypes indices for years
#' @param croptypes_cultivars_prop a data.frame containing CroptypeID, CultivarsID and proportions
#' @param basic_patho_param a list of pathogen aggressiveness parameters on a susceptible host
#' for a pathogen genotype not adapted to resistance: \itemize{
#' \item infection_rate = maximal expected infection rate of a propagule on a healthy host,
......@@ -134,7 +140,8 @@
#' @include RcppExports.R
#' @export
simul_landsepi <- function(seed = 12345, time_param = list(Nyears = 20, nTSpY = 120)
, GPKGLandscape, basic_patho_param, disp_patho, disp_host, pI0 = 5e-4
, croptype_names, cultivars, cultivars_genes_list, genes, landscape=NULL, area,
rotation, croptypes_cultivars_prop, basic_patho_param, disp_patho, disp_host, pI0 = 5e-4
, outputs=list(epid_outputs="all", evol_outputs="all", thres_breakdown=50000
, GLAnoDis=1.48315, audpc100S=0.38)
, write = TRUE, graphic = TRUE, video = FALSE, keepRawResults = FALSE) {
......@@ -144,20 +151,20 @@ simul_landsepi <- function(seed = 12345, time_param = list(Nyears = 20, nTSpY =
# croptype parameters
## TODO: extract croptype_names (doesn't work currently from GPKG)
croptypes <- getGPKGCroptypes(GPKGLandscape)
croptype_names <- croptypes$croptypeName
#croptypes <- getGPKGCroptypes(GPKGLandscape)
#croptype_names <- croptypes$croptypeName
# croptype_names <- c()
# Host parameters
cultivar <- getGPKGCultivars(GPKGLandscape)
Nhost <- nrow(cultivar)
#cultivar <- getGPKGCultivars(GPKGLandscape)
Nhost <- nrow(cultivars)
# GeneList
cultivars_genes_list <- list()
for (i in 1:Nhost) {
GeneList <- getGPKGGeneIDForCultivar(GPKGLandscape, i-1)
cultivars_genes_list[i] <- as.list(GeneList)
}
# cultivars_genes_list <- list()
# for (i in 1:Nhost) {
# GeneList <- getGPKGGeneIDForCultivar(GPKGLandscape, i-1)
# cultivars_genes_list[i] <- as.list(GeneList)
# }
# rotation <- as.data.frame(st_read(landscape$shapefilename, layer = landscape$layername_croptype))
# rotation$geom <- NULL
......@@ -166,12 +173,12 @@ simul_landsepi <- function(seed = 12345, time_param = list(Nyears = 20, nTSpY =
#cultivar_list <- DBI::dbGetQuery(inputDB, "SELECT * FROM CultivarList")
cultivars_param <- list(
Nhost = Nhost, ## TODO: we don't need this one if we use Nhost = cultivars.size() in code C
name = cultivar$cultivarName,
initial_density = cultivar$initial_density,
max_density = cultivar$max_density,
growth_rate = cultivar$growth_rate,
reproduction_rate = cultivar$reproduction_rate,
death_rate = cultivar$death_rate,
name = as.character(cultivars$cultivarName),
initial_density = as.numeric(cultivars$initial_density),
max_density = as.numeric(cultivars$max_density),
growth_rate = as.numeric(cultivars$growth_rate),
reproduction_rate = as.numeric(cultivars$reproduction_rate),
death_rate = as.numeric(cultivars$death_rate),
sigmoid_kappa_host = 0.002,
sigmoid_sigma_host = 1.001,
sigmoid_plateau_host = 1.0,
......@@ -179,20 +186,20 @@ simul_landsepi <- function(seed = 12345, time_param = list(Nyears = 20, nTSpY =
)
# Evolution parameters
gene <- getGPKGGenes(GPKGLandscape)
genes_param <- list(name = gene$geneName
, fitness_cost = gene$fitness_cost
, mutation_prob = gene$mutation_prob
, efficiency = gene$efficiency
, tradeoff_strength = gene$tradeoff_strength
, Nlevels_aggressiveness = gene$Nlevels_aggressiveness
, time_to_activ_exp = gene$time_to_activ_exp
, time_to_activ_var = gene$time_to_activ_var
, target_trait = gene$target_trait)
## genes <- getGPKGGenes(GPKGLandscape)
genes_param <- list(name = as.character(genes$geneName)
, fitness_cost = as.numeric(genes$fitness_cost)
, mutation_prob = as.numeric(genes$mutation_prob)
, efficiency = as.numeric(genes$efficiency)
, tradeoff_strength = as.numeric(genes$tradeoff_strength)
, Nlevels_aggressiveness = as.numeric(genes$Nlevels_aggressiveness)
, time_to_activ_exp = as.numeric(genes$time_to_activ_exp)
, time_to_activ_var = as.numeric(genes$time_to_activ_var)
, target_trait = as.character(genes$target_trait))
# Economic parameters
eco_param <- list(yield_perHa = cbind(H = cultivar$yield_H, L = cultivar$yield_L, I = cultivar$yield_I, R = cultivar$yield_R)
, production_cost_perHa = cultivar$production_cost, market_value = cultivar$market_value)
eco_param <- list(yield_perHa = cbind(H = as.numeric(cultivars$yield_H), L = as.numeric(cultivars$yield_L), I = as.numeric(cultivars$yield_I), R = as.numeric(cultivars$yield_R))
, production_cost_perHa = as.numeric(cultivars$production_cost), market_value = as.numeric(cultivars$market_value))
####
......@@ -202,16 +209,29 @@ simul_landsepi <- function(seed = 12345, time_param = list(Nyears = 20, nTSpY =
dispersal <- list(disp_patho = disp_patho, disp_host = disp_host)
inits <- list(pI0 = pI0)
landscape <- list(shapefilename = GPKGLandscape, layername_croptype = "croptypeID", layername_area = "area")
model_landsepi(time_param = time_param, landscape = landscape, dispersal = dispersal, inits = inits, seed = seed
, cultivars_param = cultivars_param, basic_patho_param = basic_patho_param, genes_param = genes_param)
#area <- getGPKGArea(GPKGLandscape)
#rotation <- as.matrix(getGPKGRotation(GPKGLandscape))
#croptypes_cultivars_prop <- as.matrix(getGPKGCroptypesRaw(GPKGLandscape))
#landscape <- list(shapefilename = GPKGLandscape, layername_croptype = "croptypeID", layername_area = "area")
model_landsepi(time_param = time_param,
area_vector = area,
rotation_matrix = rotation,
croptypes_cultivars_prop = as.matrix(croptypes_cultivars_prop),
dispersal = dispersal,
inits = inits,
seed = seed,
cultivars_param = cultivars_param,
basic_patho_param = basic_patho_param,
genes_param = genes_param)
### TODO Move outputs into a new method
####
#### Generate outputs
####
print("Compute model outputs")
## TODO: find another way to compute Npoly?
Npoly <- length(st_read(landscape$shapefilename, layer = landscape$layername_area)$area)
Npoly <- length(area)
#Npoly <- length(st_read(landscape$shapefilename, layer = landscape$layername_area)$area)
## Evolutionary output
Npatho <- prod(genes_param$Nlevels_aggressiveness)
......@@ -232,7 +252,7 @@ simul_landsepi <- function(seed = 12345, time_param = list(Nyears = 20, nTSpY =
, eco_product=c(0,NA)
, eco_benefit=c(0,NA)
, eco_grossmargin=c(NA,NA))
epid_res <- epid_output(outputs[["epid_outputs"]], time_param, Npatho, landscape, cultivars_param, eco_param
epid_res <- epid_output(outputs[["epid_outputs"]], time_param, Npatho, area, rotation, croptypes_cultivars_prop, cultivars_param, eco_param
, outputs[["GLAnoDis"]], ylim_param, write=write, graphic=graphic)
} else {
epid_res <- NULL
......@@ -240,7 +260,7 @@ simul_landsepi <- function(seed = 12345, time_param = list(Nyears = 20, nTSpY =
## Video
if (video & !is.null(epid_res[["audpc"]])) {
video(epid_res[["audpc"]], time_param, Npatho, landscape, croptype_names, cultivars_param, outputs[["audpc100S"]], evol_res$durability)
video(epid_res[["audpc"]], time_param, Npatho, landscape, area, rotation, croptypes_cultivars_prop, croptype_names, cultivars_param, outputs[["audpc100S"]], evol_res$durability)
}
if (!keepRawResults){
......
......@@ -8,13 +8,15 @@ epid_output(
types = "all",
time_param,
Npatho,
landscape,
area,
rotation,
croptypes,
cultivars_param,
eco_param,
GLAnoDis = cultivars_param$max_density[1],
ylim_param = list(audpc = c(0, 0.38), gla_abs = c(0, 1.48), gla_rel = c(0, 1),
eco_cost = c(0, NA), eco_product = c(0, NA), eco_benefit = c(0, NA), eco_grossmargin =
c(NA, NA)),
eco_cost = c(0, NA), eco_product = c(0, NA), eco_benefit = c(0, NA), eco_grossmargin
= c(NA, NA)),
write = TRUE,
graphic = TRUE,
path = getwd()
......@@ -42,11 +44,11 @@ related to the specified sanitary status (H, L, I or R and all their combination
\item{Npatho}{number of pathogen genotypes.}
\item{landscape}{list of landscape parameters:\itemize{
\item shapefilename = name of the landscape,
\item layername_croptype = name of the layer containing croptype indices,
\item layername_area = name of the layer containing polygon areas.
}}
\item{area}{a vector containing polygons areas}
\item{rotation}{a data.frame containing croptypes indices by years for polygons indices}
\item{croptypes}{a data.frame containing croptypeID, CulitvarsID and proportion}
\item{cultivars_param}{list of parameters associated with each host genotype (i.e. cultivars)\itemize{
\item name = vector of cultivar names,
......
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/GPKGTools.R
\name{getGPKGArea}
\alias{getGPKGArea}
\title{getGPKGArea}
\usage{
getGPKGArea(gpkgfile)
}
\arguments{
\item{gpkgfile}{a GPKG file}
}
\value{
a vector of the area of each polygons
}
\description{
get Area layer as a vector
}
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/GPKGTools.R
\name{getGPKGCroptypesRaw}
\alias{getGPKGCroptypesRaw}
\title{getGPKGCroptypesRaw}
\usage{
getGPKGCroptypesRaw(inputGPKG)
}
\arguments{
\item{inputGPKG}{a GPKG filename}
}
\value{