From c912735d580a8b7ec4991a9ecb53b7956cee24b7 Mon Sep 17 00:00:00 2001
From: Marco Matthies <71844+marcom@users.noreply.github.com>
Date: Mon, 10 Mar 2025 08:07:23 +0100
Subject: [PATCH] Allow for using multiple crop models at the same time in one
simulation
---
data/crops/aquacrop/crop_data.csv | 4 -
src/core/input.jl | 20 +++--
src/core/simulation.jl | 4 +-
src/crop/almass.jl | 2 +-
src/crop/aquacrop.jl | 129 ++++++++++++++++-----------
src/crop/cropmodels.jl | 45 ++++++----
src/crop/farmplot.jl | 7 +-
src/crop/simplecrop.jl | 4 +-
src/parameters.toml | 4 +-
test/paramscan.toml | 2 +-
test/test_parameters.toml | 2 +-
test/test_parameters_almass.toml | 2 +-
test/test_parameters_aquacrop.toml | 4 +-
test/test_parameters_simplecrop.toml | 2 +-
14 files changed, 135 insertions(+), 96 deletions(-)
diff --git a/data/crops/aquacrop/crop_data.csv b/data/crops/aquacrop/crop_data.csv
index c9a871f..89ed915 100644
--- a/data/crops/aquacrop/crop_data.csv
+++ b/data/crops/aquacrop/crop_data.csv
@@ -1,8 +1,4 @@
persefone_name,aquacrop_name,crop_group,min_sow_date,max_sow_date
-natural grass,alfalfaGDD,semi-natural,NA,NA
-permanent grassland (grazed),alfalfaGDD,grass,NA,NA
-permanent set-aside,alfalfaGDD,grass,NA,NA
-no growth,alfalfaGDD,semi-natural,NA,NA
winter rape,rapeseed,grain,20 August,25 August
winter wheat,wheatGDD,grain,15 October,31 October
maize,maize,grain,15 April,30 April
diff --git a/src/core/input.jl b/src/core/input.jl
index b149699..9437f0e 100644
--- a/src/core/input.jl
+++ b/src/core/input.jl
@@ -104,21 +104,23 @@ function preprocessparameters(settings::Dict{String,Any}, defaultoutdir::String)
Base.error("Couldn't find map directory $(settings["world.mapdirectory"]).")
end
end
- if !isdir(settings["crop.cropdirectory"])
- if isdir(abspath(pkgdir(@__MODULE__), settings["crop.cropdirectory"]))
- settings["crop.cropdirectory"] = abspath(pkgdir(@__MODULE__), settings["crop.cropdirectory"])
- @debug "Using package directory to load crop data: $(settings["crop.cropdirectory"])."
- else
- Base.error("Couldn't find crop directory $(settings["crop.cropdirectory"]).")
+ # check cropmodel settings
+ cropmodels = string.(split(settings["crop.cropmodel"], ","))
+ for cropmodel in cropmodels
+ if !(cropmodel in AVAILABLE_CROPMODELS)
+ error("cropmodel = \"$cropmodel\", but has to be one of: $AVAILABLE_CROPMODELS")
+ end
+ end
+ cropdirs = string.(split(settings["crop.cropdirectory"], ","))
+ for cropdir in cropdirs
+ if !isdir(cropdir)
+ error("Can't access crop model data dir $cropdir")
end
end
# sanity checks
if settings["core.startdate"] > settings["core.enddate"]
Base.error("Enddate is earlier than startdate.")
end
- if !(settings["crop.cropmodel"] in AVAILABLE_CROPMODELS)
- Base.error("crop.cropmodel = \"$(settings["crop.cropmodel"])\", but has to be one of: $AVAILABLE_CROPMODELS")
- end
#FIXME enable parallelisation
# if !isempty(settings["internal.scanparams"]) && (nprocs() < 2)
# @warn "To parallelise multiple simulation runs, use `julia -p <n>`."
diff --git a/src/core/simulation.jl b/src/core/simulation.jl
index 667f766..2aa11ce 100644
--- a/src/core/simulation.jl
+++ b/src/core/simulation.jl
@@ -140,8 +140,8 @@ function initmodel(settings::Dict{String, Any})
settings["world.weatherfile"]),
settings["core.startdate"],
settings["core.enddate"])
- crops = initcropmodel(settings["crop.cropmodel"],
- settings["crop.cropdirectory"])
+ crops = initcropmodels(settings["crop.cropmodel"],
+ settings["crop.cropdirectory"])
model = AgricultureModel(;
settings,
rng = StableRNG(settings["core.seed"]),
diff --git a/src/crop/almass.jl b/src/crop/almass.jl
index ecd155c..de804b3 100644
--- a/src/crop/almass.jl
+++ b/src/crop/almass.jl
@@ -417,7 +417,7 @@ function readcropparameters(cropdirectory::String)
growthdata = CSV.File(joinpath(cropdirectory, GROWTHFILE), missingstring="NA",
types=[Int,String,String,GrowthPhase,String,
Float64,Float64,Float64,Float64])
- croptypes = Dict{String,CropType}()
+ croptypes = Dict{String,AbstractCropType}()
for crop in cropdata
cropgrowthdata = growthdata |> filter(x -> !ismissing(x.crop_name) &&
x.crop_name == crop.name)
diff --git a/src/crop/aquacrop.jl b/src/crop/aquacrop.jl
index f95e516..d9c9f1c 100644
--- a/src/crop/aquacrop.jl
+++ b/src/crop/aquacrop.jl
@@ -91,7 +91,7 @@ needed for Persefone (cropgroup, minsowdate, maxsowdate).
"""
function readcropparameters(cropdirectory::String)
@debug "Reading extra crop parameters for AquaCrop crop model"
- croptypes = Dict{String,CropType}()
+ croptypes = Dict{String, AbstractCropType}()
df = CSV.read(joinpath(cropdirectory, CROPFILE), DataFrame;
missingstring="NA", types=[String, String, String, AnnualDate, AnnualDate])
for row in eachrow(df)
@@ -105,61 +105,68 @@ function readcropparameters(cropdirectory::String)
return croptypes
end
+function make_aquacrop_cropfield(croptype::CropType,
+ model::SimulationModel,
+ soiltype::SoilType)
+ # TODO: get this mapping for soil types from a CSV file?
+ soiltype_to_aquacrop = Dict(soiltype => replace(string(soiltype), "_" => " ")
+ for soiltype in instances(SoilType))
+ aquacrop_soiltype = soiltype_to_aquacrop[soiltype]
+
+ start_date = model.date # TODO: maybe `today(model)`
+ end_date = model.weather.lastdate # TODO: maybe `lastdate(model)`
+ start_daynum = daynumber(model.weather, start_date)
+ end_daynum = daynumber(model.weather, end_date)
+ dayrange = start_daynum:end_daynum
+ input_Tmin = @view model.weather.mintemp[dayrange]
+ input_Tmax = @view model.weather.maxtemp[dayrange]
+ input_ETo = @view model.weather.evapotranspiration[dayrange]
+ input_Rain = @view model.weather.precipitation[dayrange]
+
+ # Generate the keyword object for the AquaCrop simulation
+ kwargs = (
+ ## Necessary keywords
+ # runtype
+ runtype = AquaCrop.NoFileRun(),
+ # project input
+ Simulation_DayNr1 = start_date,
+ Simulation_DayNrN = end_date,
+ Crop_Day1 = start_date, # TODO: originally start_date + Week(1)
+ Crop_DayN = end_date,
+ # soil
+ soil_type = aquacrop_soiltype,
+ # crop
+ crop_type = croptype.aquacrop_name,
+ # Climate
+ InitialClimDate = start_date,
+ ## Optional keyworkds
+ # Climate
+ Tmin = input_Tmin,
+ Tmax = input_Tmax,
+ ETo = input_ETo,
+ Rain = input_Rain,
+ # change soil properties
+ soil_layers = Dict("Thickness" => 5.0)
+ )
+
+ aquacrop_cropfield, all_ok = AquaCrop.start_cropfield(; kwargs...)
+ if ! all_ok.logi
+ @error "Failed calling AquaCrop.start_cropfield()\nAquaCrop error: $(all_ok.msg)"
+ end
+
+ return aquacrop_cropfield
+end
+
mutable struct CropState <: AbstractCropState
croptype::CropType
height::Tlength # TODO: remove height field, supply from cropstate
soiltype::SoilType
cropstate::AquaCrop.AquaCropField
- function CropState(crop_type::CropType, soiltype::SoilType,
+ function CropState(croptype::CropType, soiltype::SoilType,
model::SimulationModel, height::Tlength=0.0cm)
- # TODO: get this mapping for soil types from a CSV file?
- soiltype_to_aquacrop = Dict(soiltype => replace(string(soiltype), "_" => " ")
- for soiltype in instances(SoilType))
- aquacrop_soiltype = soiltype_to_aquacrop[soiltype]
-
- start_date = model.date # TODO: maybe `today(model)`
- end_date = model.weather.lastdate # TODO: maybe `lastdate(model)`
- start_daynum = daynumber(model.weather, start_date)
- end_daynum = daynumber(model.weather, end_date)
- dayrange = start_daynum:end_daynum
- input_Tmin = @view model.weather.mintemp[dayrange]
- input_Tmax = @view model.weather.maxtemp[dayrange]
- input_ETo = @view model.weather.evapotranspiration[dayrange]
- input_Rain = @view model.weather.precipitation[dayrange]
-
- # Generate the keyword object for the AquaCrop simulation
- kwargs = (
- ## Necessary keywords
- # runtype
- runtype = AquaCrop.NoFileRun(),
- # project input
- Simulation_DayNr1 = start_date,
- Simulation_DayNrN = end_date,
- Crop_Day1 = start_date, # TODO: originally start_date + Week(1)
- Crop_DayN = end_date,
- # soil
- soil_type = aquacrop_soiltype,
- # crop
- crop_type = crop_type.aquacrop_name,
- # Climate
- InitialClimDate = start_date,
- ## Optional keyworkds
- # Climate
- Tmin = input_Tmin,
- Tmax = input_Tmax,
- ETo = input_ETo,
- Rain = input_Rain,
- # change soil properties
- soil_layers = Dict("Thickness" => 5.0)
- )
-
- aquacrop_cropfield, all_ok = AquaCrop.start_cropfield(; kwargs...)
- if ! all_ok.logi
- @error "Failed calling AquaCrop.start_cropfield()\nAquaCrop error: $(all_ok.msg)"
- end
-
- return new(crop_type, height, soiltype, aquacrop_cropfield)
+ aquacrop_cropfield = make_aquacrop_cropfield(croptype, model, soiltype)
+ return new(croptype, height, soiltype, aquacrop_cropfield)
end
end
@@ -168,9 +175,12 @@ cropname(cs::CropState) = cropname(croptype(cs))
cropheight(cs::CropState) = cs.height # TODO: calculate from AquaCrop state info
cropcover(cs::CropState) = AquaCrop.canopycover(cs.cropstate)
cropyield(cs::CropState) = AquaCrop.dryyield(cs.cropstate) # TODO: there is also freshyield
-isharvestable(cs::CropState) = true # TODO: implement this correctly
-makecropstate(crop_type::CropType, model::SimulationModel, soiltype::SoilType) =
- CropState(crop_type, soiltype, model)
+function isharvestable(cs::CropState)
+ stage = get_aquacrop_stage(cs)
+ return !ismissing(stage) && stage == 0
+end
+makecropstate(croptype::CropType, model::SimulationModel, soiltype::SoilType) =
+ CropState(croptype, soiltype, model)
function setsoiltype!(cs::CropState, soiltype::SoilType)
# TODO: this does not affect the actual soil type of the state of
# the aquacrop simulation
@@ -178,6 +188,11 @@ function setsoiltype!(cs::CropState, soiltype::SoilType)
return nothing
end
+function get_aquacrop_stage(cs::CropState)
+ stages = cs.cropstate.dayout.Stage
+ length(stages) > 0 ? last(stages) : missing
+end
+
"""
stepagent!(cropstate, model)
@@ -196,9 +211,15 @@ function sow!(cs::CropState, model::SimulationModel, cropname::String)
if cropname ∉ keys(model.crops)
@error "cropname \"$cropname\" not found"
end
- cs.croptype = model.crops[cropname]
+ new_croptype = model.crops[cropname]
+ if typeof(new_croptype) !== CropType
+ error("Cannot change crop model of AquaCrop crop state.")
+ end
+
+ cs.croptype = new_croptype
cs.height = 0.0cm
- # TODO: other things needed for AquaCrop?
+ # cs.soiltype stays the way it is
+ cs.cropstate = make_aquacrop_cropfield(cs.croptype, model, soiltype)
end
"""
diff --git a/src/crop/cropmodels.jl b/src/crop/cropmodels.jl
index 153eb39..5aafe24 100644
--- a/src/crop/cropmodels.jl
+++ b/src/crop/cropmodels.jl
@@ -6,24 +6,39 @@
using StatsBase: mode
"""
- initcropmodel(cropmodel, cropdirectory)
+ initcropmodels(cropmodels, cropdirectory)
-Initialise the crop model. Returns the crop types available in the
-simulation, as well as the types `Tcroptype` and `Tcropstate`.
+Initialise the crop models given as a comma-delimited string
+(e.g. "almass,aquacrop"). The crop model parameters are read from the
+`cropdirectories`, also a comma-delimited string. Returns the crop
+types available in the simulation.
"""
-function initcropmodel(cropmodel::AbstractString, cropdirectory::AbstractString)
- if cropmodel == "almass"
- crops = ALMaSS.readcropparameters(cropdirectory)
- elseif cropmodel == "simple"
- crops_almass = ALMaSS.readcropparameters(cropdirectory)
- crops = Dict(name => SimpleCrop.CropType(ct.name, ct.group, ct.minsowdate, ct.maxsowdate)
- for (name, ct) in crops_almass)
- elseif cropmodel == "aquacrop"
- crops = AquaCropWrapper.readcropparameters(cropdirectory)
- else
- error("initcropmodel: no implementation for crop model '$cropmodel'")
+function initcropmodels(cropmodels::AbstractString,
+ cropdirectories::AbstractString)
+ models = string.(split(cropmodels, ","))
+ dirs = string.(split(cropdirectories, ","))
+ if length(models) != length(dirs)
+ error("Different number of cropmodels and cropdirectories.\n" *
+ "cropmodels = $cropmodels\ncropdirectories = $cropdirectories")
end
- return crops
+
+ allcrops = Dict{String, AbstractCropType}()
+ for (i, cropmodel) in enumerate(models)
+ cropdirectory = dirs[i]
+ if cropmodel == "almass"
+ crops = ALMaSS.readcropparameters(cropdirectory)
+ elseif cropmodel == "simple"
+ crops_almass = ALMaSS.readcropparameters(cropdirectory)
+ crops = Dict(name => SimpleCrop.CropType(ct.name, ct.group, ct.minsowdate, ct.maxsowdate)
+ for (name, ct) in crops_almass)
+ elseif cropmodel == "aquacrop"
+ crops = AquaCropWrapper.readcropparameters(cropdirectory)
+ else
+ error("initcropmodel: no implementation for crop model '$cropmodel'")
+ end
+ merge!(allcrops, crops)
+ end
+ return allcrops
end
"""
diff --git a/src/crop/farmplot.jl b/src/crop/farmplot.jl
index 40a3917..b39e5bd 100644
--- a/src/crop/farmplot.jl
+++ b/src/crop/farmplot.jl
@@ -44,8 +44,13 @@ end
Sow the specified crop on the farmplot.
"""
function sow!(farmplot::FarmPlot, model::SimulationModel, cropname::String)
+ new_croptype = model.crops[cropname]
createevent!(model, farmplot.pixels, sowing)
- sow!(farmplot.cropstate, model, cropname)
+ if typeof(new_croptype) === typeof(croptype(farmplot))
+ sow!(farmplot.cropstate, model, cropname)
+ else
+ farmplot.cropstate = makecropstate(new_croptype, model, farmplot.soiltype)
+ end
@debug "Farmer $(farmplot.farmer) sowed $(cropname) on farmplot $(farmplot.id)."
end
diff --git a/src/crop/simplecrop.jl b/src/crop/simplecrop.jl
index 7f74ca4..2fd9c29 100644
--- a/src/crop/simplecrop.jl
+++ b/src/crop/simplecrop.jl
@@ -44,8 +44,8 @@ cropheight(cs::CropState) = cs.height
cropcover(cs::CropState) = 0.0
cropyield(cs::CropState) = 0.0 # TODO: units?
isharvestable(cs::CropState) = true
-makecropstate(crop_type::CropType, model::SimulationModel, soiltype::SoilType) =
- CropState(crop_type, 0.0cm)
+makecropstate(croptype::CropType, model::SimulationModel, soiltype::SoilType) =
+ CropState(croptype, 0.0cm)
"""
diff --git a/src/parameters.toml b/src/parameters.toml
index c3d8aff..1b7d420 100644
--- a/src/parameters.toml
+++ b/src/parameters.toml
@@ -46,5 +46,5 @@ indoutfreq = "end" # output frequency individual-level data, daily/monthly/yearl
insectmodel = ["season", "habitat", "pesticides", "weather"] # factors affecting insect growth
[crop]
-cropmodel = "aquacrop" # crop growth model to use, "almass", "aquacrop", or "simple"
-cropdirectory = "data/crops/aquacrop/" # the directory storing all data files for the selected crop model
+cropmodel = "almass,aquacrop" # crop growth model to use, "almass", "aquacrop", or "simple"
+cropdirectory = "data/crops/almass/,data/crops/aquacrop/" # the directory storing all data files for the selected crop model
diff --git a/test/paramscan.toml b/test/paramscan.toml
index 565e685..fecac92 100644
--- a/test/paramscan.toml
+++ b/test/paramscan.toml
@@ -29,4 +29,4 @@ indoutfreq = "end" # output frequency individual-level data, daily/monthly/yearl
[crop]
cropmodel = "simple" # crop growth model to use, "simple" or "almass"
-
+cropdirectory = "../data/crops/almass/"
diff --git a/test/test_parameters.toml b/test/test_parameters.toml
index fa539d9..a9d6b9b 100644
--- a/test/test_parameters.toml
+++ b/test/test_parameters.toml
@@ -42,4 +42,4 @@ insectmodel = ["season", "habitat", "pesticides"] # which factors affect insect
[crop]
cropmodel = "almass" # crop growth model to use, "almass", "aquacrop", or "simple"
-cropdirectory = "data/crops/almass/" # the directory storing all data files for the selected crop model
+cropdirectory = "../data/crops/almass/" # the directory storing all data files for the selected crop model
diff --git a/test/test_parameters_almass.toml b/test/test_parameters_almass.toml
index 6d8fa7e..c497679 100644
--- a/test/test_parameters_almass.toml
+++ b/test/test_parameters_almass.toml
@@ -42,4 +42,4 @@ insectmodel = ["season", "habitat", "pesticides"] # which factors affect insect
[crop]
cropmodel = "almass" # crop growth model to use, "almass", "aquacrop", or "simple"
-cropdirectory = "data/crops/almass/" # the directory storing all data files for the selected crop model
+cropdirectory = "../data/crops/almass/" # the directory storing all data files for the selected crop model
diff --git a/test/test_parameters_aquacrop.toml b/test/test_parameters_aquacrop.toml
index 9f29fdc..13d5935 100644
--- a/test/test_parameters_aquacrop.toml
+++ b/test/test_parameters_aquacrop.toml
@@ -41,5 +41,5 @@ indoutfreq = "daily" # output frequency individual-level data, daily/monthly/yea
insectmodel = ["season", "habitat", "pesticides"] # which factors affect insect growth ("weather" is not yet implemented)
[crop]
-cropmodel = "aquacrop" # crop growth model to use, "almass", "aquacrop", or "simple"
-cropdirectory = "data/crops/aquacrop/"
+cropmodel = "almass,aquacrop" # crop growth model to use, "almass", "aquacrop", or "simple"
+cropdirectory = "../data/crops/almass/,../data/crops/aquacrop/"
diff --git a/test/test_parameters_simplecrop.toml b/test/test_parameters_simplecrop.toml
index 6540f38..f2fd15c 100644
--- a/test/test_parameters_simplecrop.toml
+++ b/test/test_parameters_simplecrop.toml
@@ -42,4 +42,4 @@ insectmodel = ["season", "habitat", "pesticides"] # which factors affect insect
[crop]
cropmodel = "simple" # crop growth model to use, "almass", "aquacrop", or "simple"
-cropdirectory = "data/crops/almass/" # the directory storing all data files for the selected crop model
+cropdirectory = "../data/crops/almass/" # the directory storing all data files for the selected crop model
--
GitLab