diff --git a/src/Persefone.jl b/src/Persefone.jl
index ec269bfd733887472526f65acaad82a1e779f3e5..ad32aef3a8e66814209fc56b7ab69397b74c40c8 100644
--- a/src/Persefone.jl
+++ b/src/Persefone.jl
@@ -100,7 +100,12 @@ export
populationtrends,
visualiseoutput,
savemodelobject,
- loadmodelobject
+ loadmodelobject,
+ croptype,
+ cropname,
+ cropheight,
+ cropcover,
+ cropyield
## Import and define units and dimensions
import Unitful: cm, m, km, ha, mg, g, kg, Length, Area, Mass
@@ -127,6 +132,8 @@ The supertype of all agents in the model (animal species, farmer types, farmplot
"""
abstract type ModelAgent end
+function stepagent! end
+
## include all module files (note that the order matters - if file
## b references something from file a, it must be included later)
include("core/input.jl")
@@ -136,8 +143,11 @@ include("analysis/makieplots.jl")
include("world/landscape.jl")
include("world/weather.jl")
-include("crop/crops.jl")
include("crop/farmplot.jl")
+include("crop/cropmodels.jl")
+include("crop/almass.jl")
+include("crop/simplecrop.jl")
+
include("farm/farm.jl")
include("nature/insects.jl")
diff --git a/src/core/input.jl b/src/core/input.jl
index a8bc2a42d392c2ab220590867d791a53644fd9fe..46460cb92af1701ff1109771a439809722bb44a5 100644
--- a/src/core/input.jl
+++ b/src/core/input.jl
@@ -17,6 +17,11 @@ const PARAMFILE = joinpath(pkgdir(Persefone), "src/parameters.toml")
#XXX do I need to use absolute paths for all input files in case working dir is changed?
# (can be done with `joinpath(dirname(@__FILE__), <filename>)` )
+"""
+The crop models that can be used in the simulation.
+"""
+const AVAILABLE_CROPMODELS = ["almass", "simple"]
+
"""
@param(domainparam)
@@ -89,6 +94,9 @@ function preprocessparameters(settings::Dict{String,Any}, defaultoutdir::String)
Base.error("Enddate is earlier than startdate.") #TODO replace with exception
end
settings["world.mapresolution"] = settings["world.mapresolution"] * 1m
+ if ! (settings["crop.cropmodel"] in AVAILABLE_CROPMODELS)
+ 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 05158aa6590a6e4bcc38941ca6332d45f734b778..ac9a4570dfe97298e33f1cf582162cb4236daccd 100644
--- a/src/core/simulation.jl
+++ b/src/core/simulation.jl
@@ -12,7 +12,7 @@ This is the heart of the model - a struct that holds all data and state
for one simulation run. It is created by [`initialise`](@ref) and passed
as input to most model functions.
"""
-mutable struct AgricultureModel <: SimulationModel
+mutable struct AgricultureModel{Tcroptype,Tcropstate} <: SimulationModel
settings::Dict{String,Any}
rng::AbstractRNG
logger::AbstractLogger
@@ -21,9 +21,9 @@ mutable struct AgricultureModel <: SimulationModel
date::Date
landscape::Matrix{Pixel}
weather::Dict{Date,Weather}
- crops::Dict{String,CropType}
+ crops::Dict{String,Tcroptype}
farmers::Vector{Farmer}
- farmplots::Vector{FarmPlot}
+ farmplots::Vector{FarmPlot{Tcropstate}}
animals::Vector{Union{Animal,Nothing}}
migrants::Vector{Pair{Animal,Date}}
events::Vector{FarmEvent}
@@ -117,23 +117,29 @@ function initmodel(settings::Dict{String, Any})
settings["world.weatherfile"]),
settings["core.startdate"],
settings["core.enddate"])
- crops = readcropparameters(settings["crop.cropfile"],
- settings["crop.growthfile"])
- model = AgricultureModel(settings,
- StableRNG(settings["core.seed"]),
- logger,
- Vector{DataOutput}(),
- Dict{String, DataFrame}(),
- settings["core.startdate"],
- landscape,
- weather,
- crops,
- Vector{Farmer}(),
- Vector{FarmPlot}(),
- Vector{Union{Animal,Nothing}}(),
- Vector{Pair{Animal, Date}}(),
- Vector{FarmEvent}())
+ crops, Tcroptype, Tcropstate = initcropmodel(settings["crop.cropmodel"],
+ settings["crop.cropfile"],
+ settings["crop.growthfile"])
+ farmers = Vector{Farmer}()
+ farmplots = Vector{FarmPlot{Tcropstate}}()
+ model = AgricultureModel{Tcroptype,Tcropstate}(
+ settings,
+ StableRNG(settings["core.seed"]),
+ logger,
+ Vector{DataOutput}(),
+ Dict{String, DataFrame}(),
+ settings["core.startdate"],
+ landscape,
+ weather,
+ crops,
+ farmers,
+ farmplots,
+ Vector{Union{Animal,Nothing}}(),
+ Vector{Pair{Animal, Date}}(),
+ Vector{FarmEvent}()
+ )
saveinputfiles(model)
+
initfields!(model)
initfarms!(model)
initnature!(model)
diff --git a/src/crop/almass.jl b/src/crop/almass.jl
new file mode 100644
index 0000000000000000000000000000000000000000..5ac4c53ed22c13ab66fe9125577659fd3879b838
--- /dev/null
+++ b/src/crop/almass.jl
@@ -0,0 +1,258 @@
+### Persefone.jl - a model of agricultural landscapes and ecosystems in Europe.
+###
+### This file is responsible for managing the crop growth modules.
+###
+
+#TODO write tests for input functions
+
+module ALMaSS
+
+using Persefone:
+ Management,
+ Length,
+ cm,
+ SimulationModel,
+ fertiliser,
+ maxtemp,
+ mintemp
+
+import Persefone:
+ stepagent!,
+ croptype,
+ cropname,
+ cropheight,
+ cropcover,
+ cropyield
+
+using Dates: Date, month, monthday
+using CSV: CSV
+
+"""
+ GrowthPhase
+
+ALMaSS crop growth curves are split into five phases, triggered by
+seasonal dates or agricultural events.
+"""
+@enum GrowthPhase janfirst sow marchfirst harvest1 harvest2
+
+"""
+ CropCurveParams
+
+The values in this struct define one crop growth curve.
+"""
+struct CropCurveParams
+ #TODO add Unitful
+ curveID::Int
+ highnutrients::Bool
+ GDD::Dict{GrowthPhase,Vector{Float64}}
+ LAItotal::Dict{GrowthPhase,Vector{Float64}}
+ LAIgreen::Dict{GrowthPhase,Vector{Float64}}
+ height::Dict{GrowthPhase,Vector{Length{Float64}}}
+end
+
+"""
+ CropType
+
+The type struct for all crops. Currently follows the crop growth model as
+implemented in ALMaSS.
+"""
+struct CropType
+ name::String
+ minsowdate::Union{Missing,Date}
+ maxsowdate::Union{Missing,Date}
+ minharvestdate::Union{Missing,Date}
+ maxharvestdate::Union{Missing,Date}
+ mingrowthtemp::Union{Missing,Float64}
+ highnutrientgrowth::Union{Missing,CropCurveParams}
+ lownutrientgrowth::Union{Missing,CropCurveParams}
+ #issowable::Union{Function,Bool}
+end
+
+cropname(ct::CropType) = ct.name
+
+"""
+ CropState
+
+The state data for an ALMaSS vegetation point calculation. Usually
+part of a `FarmPlot`.
+"""
+mutable struct CropState
+ #TODO add Unitful
+ croptype::CropType
+ phase::GrowthPhase
+ growingdegreedays::Float64
+ height::Length{Float64}
+ LAItotal::Float64
+ LAIgreen::Float64
+ #biomass::Float64 #XXX I need to figure out how to calculate this
+ events::Vector{Management}
+end
+
+croptype(cs::CropState) = cs.croptype
+cropname(cs::CropState) = cropname(croptype(cs))
+cropheight(cs::CropState) = cs.height
+cropcover(cs::CropState) = 0.0 # TODO: related to LAItotal, LAIgreen?
+cropyield(cs::CropState) = 0.0 # TODO: units? needs biomass?
+
+"""
+ Base.tryparse(type, str)
+
+Extend `tryparse` to allow parsing GrowthPhase values.
+(Needed to read in the CSV parameter file.)
+"""
+function Base.tryparse(type::Type{GrowthPhase}, str::String)
+ str == "janfirst" ? janfirst :
+ str == "sow" ? sow :
+ str == "marchfirst" ? marchfirst :
+ str == "harvest1" ? harvest1 :
+ str == "harvest2" ? harvest2 :
+ nothing
+end
+
+"""
+ buildgrowthcurve(data)
+
+Convert a list of rows from the crop growth data into a CropCurveParams object.
+"""
+function buildgrowthcurve(data::Vector{CSV.Row})
+ isempty(data) && return missing
+ GDD = Dict(janfirst=>Vector{Float64}(), sow=>Vector{Float64}(),
+ marchfirst=>Vector{Float64}(), harvest1=>Vector{Float64}(),
+ harvest2=>Vector{Float64}())
+ LAItotal = Dict(janfirst=>Vector{Float64}(), sow=>Vector{Float64}(),
+ marchfirst=>Vector{Float64}(), harvest1=>Vector{Float64}(),
+ harvest2=>Vector{Float64}())
+ LAIgreen = Dict(janfirst=>Vector{Float64}(), sow=>Vector{Float64}(),
+ marchfirst=>Vector{Float64}(), harvest1=>Vector{Float64}(),
+ harvest2=>Vector{Float64}())
+ height = Dict(janfirst=>Vector{Length{Float64}}(), sow=>Vector{Length{Float64}}(),
+ marchfirst=>Vector{Length{Float64}}(), harvest1=>Vector{Length{Float64}}(),
+ harvest2=>Vector{Length{Float64}}())
+ for e in data
+ append!(GDD[e.growth_phase], e.GDD)
+ append!(LAItotal[e.growth_phase], e.LAI_total)
+ append!(LAIgreen[e.growth_phase], e.LAI_green)
+ append!(height[e.growth_phase], e.height * cm) # assumes `height` is in units of `cm`
+ end
+ CropCurveParams(data[1].curve_id, data[1].nutrient_status=="high",
+ GDD, LAItotal, LAIgreen, height)
+end
+
+"""
+ readcropparameters(generalcropfile, cropgrowthfile)
+
+Parse a CSV file containing the required parameter values for each crop
+(as produced from the original ALMaSS file by `convert_almass_data.py`).
+"""
+function readcropparameters(generalcropfile::String, growthfile::String)
+ @debug "Reading crop parameters"
+ cropdata = CSV.File(generalcropfile, missingstring="NA", dateformat="d U",
+ types=[String,Date,Date,Date,Date,Float64])
+ growthdata = CSV.File(growthfile, missingstring="NA",
+ types=[Int,String,String,GrowthPhase,String,
+ Float64,Float64,Float64,Float64])
+ croptypes = Dict{String,CropType}()
+ for crop in cropdata
+ cropgrowthdata = growthdata |> filter(x -> !ismissing(x.crop_name) &&
+ x.crop_name == crop.name)
+ highnuts = buildgrowthcurve(cropgrowthdata |>
+ filter(x -> x.nutrient_status=="high"))
+ lownuts = buildgrowthcurve(cropgrowthdata |>
+ filter(x -> x.nutrient_status=="low"))
+ croptypes[crop.name] = CropType(crop.name, crop.minsowdate, crop.maxsowdate,
+ crop.minharvestdate, crop.maxharvestdate,
+ crop.mingrowthtemp, highnuts, lownuts)
+ end
+ croptypes
+end
+
+"""
+ stepagent!(cropstate, model)
+
+Update a farm plot by one day.
+"""
+function stepagent!(cs::CropState, model::SimulationModel)
+ # update growing degree days
+ # if no crop-specific base temperature is given, default to 5°C
+ # (https://www.eea.europa.eu/publications/europes-changing-climate-hazards-1/heat-and-cold/heat-and-cold-2014-mean)
+ basetemp = cs.croptype.mingrowthtemp
+ ismissing(basetemp) && (basetemp = 5.0)
+ gdd = (maxtemp(model)+mintemp(model))/2 - basetemp
+ gdd > 0 && (cs.growingdegreedays += gdd)
+ # update the phase on key dates
+ monthday(model.date) == (1,1) && (cs.phase = ALMaSS.janfirst)
+ monthday(model.date) == (3,1) && (cs.phase = ALMaSS.marchfirst)
+ # update crop growth
+ growcrop!(cs, model)
+end
+
+
+## CROP MANAGEMENT AND GROWTH FUNCTIONS
+
+"""
+ sow!(cropstate, model, cropname)
+
+Change the cropstate to sow the specified crop.
+"""
+function sow!(cs::CropState, model::SimulationModel, cropname::String)
+ #XXX test if the crop is sowable?
+ cs.croptype = model.crops[cropname]
+ cs.phase = ALMaSS.sow
+end
+
+"""
+ harvest!(cropstate, model)
+
+Harvest the crop of this cropstate.
+"""
+function harvest!(cs::CropState, model::SimulationModel)
+ cs.phase in [ALMaSS.harvest1, ALMaSS.harvest2] ?
+ cs.phase = ALMaSS.harvest2 :
+ cs.phase = ALMaSS.harvest1
+ # height & LAI will be automatically adjusted by the growth function
+ #TODO calculate and return yield
+end
+
+#TODO fertilise!()
+#TODO spray!()
+#TODO till!()
+
+"""
+ growcrop!(cropstate, model)
+
+Apply the relevant crop growth model to update the plants crop state
+on this farm plot. Implements the ALMaSS crop growth model by Topping
+et al.
+"""
+function growcrop!(cs::CropState, model::SimulationModel)
+ fertiliser in cs.events ?
+ curve = cs.croptype.lownutrientgrowth :
+ curve = cs.croptype.highnutrientgrowth
+ points = curve.GDD[cs.phase]
+ for p in 1:length(points)
+ if points[p] == 99999
+ return # the marker that there is no further growth this phase
+ elseif points[p] == -1 # the marker to set all variables to specified values
+ cs.height = curve.height[cs.phase][p]
+ cs.LAItotal = curve.LAItotal[cs.phase][p]
+ cs.LAIgreen = curve.LAIgreen[cs.phase][p]
+ return
+ else
+ gdd = cs.growingdegreedays
+ # figure out which is the correct slope value to use for growth
+ if p == length(points) || gdd < points[p+1]
+ cs.height += curve.height[cs.phase][p]
+ cs.LAItotal += curve.LAItotal[cs.phase][p]
+ cs.LAIgreen += curve.LAIgreen[cs.phase][p]
+ return
+ end
+ #XXX To be precise, we ought to figure out if one or more inflection
+ # points have been passed between yesterday and today, and calculate the
+ # growth exactly up to the inflection point with the old slope, and onward
+ # with the new slope. Not doing so will introduce a small amount of error,
+ # although I think this is acceptable.
+ end
+ end
+end
+
+end # module ALMaSS
diff --git a/src/crop/cropmodels.jl b/src/crop/cropmodels.jl
new file mode 100644
index 0000000000000000000000000000000000000000..bbbd92529424792d327ff7f5e3d55cfd0dfc5b6b
--- /dev/null
+++ b/src/crop/cropmodels.jl
@@ -0,0 +1,83 @@
+### Persefone.jl - a model of agricultural landscapes and ecosystems in Europe.
+###
+### Functions for switchable crop models.
+###
+
+"""
+ initfields!(cropmodel, cropfile, growthfile)
+
+Initialise the crop model. Returns the crop types available in the
+simulation, as well as the types `Tcroptype` and `Tcropstate`.
+"""
+function initcropmodel(cropmodel::AbstractString, cropfile::AbstractString, growthfile::AbstractString)
+ if cropmodel == "almass"
+ Tcroptype = ALMaSS.CropType
+ Tcropstate = ALMaSS.CropState
+ crops = ALMaSS.readcropparameters(cropfile, growthfile)
+ elseif cropmodel == "simple"
+ Tcroptype = SimpleCrop.CropType
+ Tcropstate = SimpleCrop.CropState
+ crops_almass = ALMaSS.readcropparameters(cropfile, growthfile)
+ crops = Dict(name => SimpleCrop.CropType(ct.name) for (name, ct) in crops_almass)
+ else
+ error("initcropmodel: no implementation for crop model '$cropmodel'")
+ end
+ return crops, Tcroptype, Tcropstate
+end
+
+"""
+ initfields!(model)
+
+Initialise the farm plots in the simulation model.
+"""
+function initfields!(model::SimulationModel)
+ n = 0
+ convertid = Dict{Int64,Int64}()
+ width, height = size(model.landscape)
+ for x in 1:width
+ for y in 1:height
+ # for each pixel, we need to extract the field ID given by the map input
+ # file, and convert it into the internal object ID used by Agents.jl,
+ # creating a new agent object if necessary
+ rawid = model.landscape[x,y].fieldid
+ (ismissing(rawid)) && continue
+ if rawid in keys(convertid)
+ objectid = convertid[rawid]
+ model.landscape[x,y].fieldid = objectid
+ push!(model.farmplots[objectid].pixels, (x,y))
+ else
+ cropstate = make_cropstate(model, @param(crop.cropmodel))
+ fp = FarmPlot(
+ length(model.farmplots) + 1,
+ [(x, y)],
+ cropstate
+ )
+ push!(model.farmplots, fp)
+ model.landscape[x,y].fieldid = fp.id
+ convertid[rawid] = fp.id
+ n += 1
+ end
+ end
+ end
+ @info "Initialised $n farm plots."
+end
+
+# internal utility function
+function make_cropstate(model::SimulationModel, cropmodel::AbstractString)
+ if cropmodel == "almass"
+ phase = (month(model.date) < 3 ? ALMaSS.janfirst : ALMaSS.marchfirst)
+ cs = ALMaSS.CropState(
+ model.crops["natural grass"],
+ phase,
+ 0.0, 0.0m, 0.0, 0.0, Vector{Management}()
+ )
+ elseif cropmodel == "simple"
+ cs = SimpleCrop.CropState(
+ model.crops["natural grass"],
+ 0.0m
+ )
+ else
+ error("Unhandled crop model '$cropmodel' in make_cropstate")
+ end
+ return cs
+end
diff --git a/src/crop/crops.jl b/src/crop/crops.jl
deleted file mode 100644
index 996c6dabd0187189e6dd701e62669598003f4f63..0000000000000000000000000000000000000000
--- a/src/crop/crops.jl
+++ /dev/null
@@ -1,121 +0,0 @@
-### Persefone.jl - a model of agricultural landscapes and ecosystems in Europe.
-###
-### This file is responsible for managing the crop growth modules.
-###
-
-#TODO write tests for input functions
-
-"""
- GrowthPhase
-
-ALMaSS crop growth curves are split into five phases, triggered by
-seasonal dates or agricultural events.
-"""
-@enum GrowthPhase janfirst sow marchfirst harvest1 harvest2
-
-"""
- CropCurveParams
-
-The values in this struct define one crop growth curve.
-"""
-struct CropCurveParams
- #TODO add Unitful
- curveID::Int
- highnutrients::Bool
- GDD::Dict{GrowthPhase,Vector{Float64}}
- LAItotal::Dict{GrowthPhase,Vector{Float64}}
- LAIgreen::Dict{GrowthPhase,Vector{Float64}}
- height::Dict{GrowthPhase,Vector{Float64}}
-end
-
-"""
- CropType
-
-The type struct for all crops. Currently follows the crop growth model as
-implemented in ALMaSS.
-"""
-struct CropType
- #TODO make this into an abstract type and create subtypes for different
- # crop submodels (#70)
- name::String
- minsowdate::Union{Missing,Date}
- maxsowdate::Union{Missing,Date}
- minharvestdate::Union{Missing,Date}
- maxharvestdate::Union{Missing,Date}
- mingrowthtemp::Union{Missing,Float64}
- highnutrientgrowth::Union{Missing,CropCurveParams}
- lownutrientgrowth::Union{Missing,CropCurveParams}
- #issowable::Union{Function,Bool}
-end
-
-"""
- Base.tryparse(type, str)
-
-Extend `tryparse` to allow parsing GrowthPhase values.
-(Needed to read in the CSV parameter file.)
-"""
-function Base.tryparse(type::Type{GrowthPhase}, str::String)
- str == "janfirst" ? janfirst :
- str == "sow" ? sow :
- str == "marchfirst" ? marchfirst :
- str == "harvest1" ? harvest1 :
- str == "harvest2" ? harvest2 :
- nothing
-end
-
-"""
- buildgrowthcurve(data)
-
-Convert a list of rows from the crop growth data into a CropCurveParams object.
-"""
-function buildgrowthcurve(data::Vector{CSV.Row})
- isempty(data) && return missing
- GDD = Dict(janfirst=>Vector{Float64}(), sow=>Vector{Float64}(),
- marchfirst=>Vector{Float64}(), harvest1=>Vector{Float64}(),
- harvest2=>Vector{Float64}())
- LAItotal = Dict(janfirst=>Vector{Float64}(), sow=>Vector{Float64}(),
- marchfirst=>Vector{Float64}(), harvest1=>Vector{Float64}(),
- harvest2=>Vector{Float64}())
- LAIgreen = Dict(janfirst=>Vector{Float64}(), sow=>Vector{Float64}(),
- marchfirst=>Vector{Float64}(), harvest1=>Vector{Float64}(),
- harvest2=>Vector{Float64}())
- height = Dict(janfirst=>Vector{Float64}(), sow=>Vector{Float64}(),
- marchfirst=>Vector{Float64}(), harvest1=>Vector{Float64}(),
- harvest2=>Vector{Float64}())
- for e in data
- append!(GDD[e.growth_phase], e.GDD)
- append!(LAItotal[e.growth_phase], e.LAI_total)
- append!(LAIgreen[e.growth_phase], e.LAI_green)
- append!(height[e.growth_phase], e.height)
- end
- CropCurveParams(data[1].curve_id, data[1].nutrient_status=="high",
- GDD, LAItotal, LAIgreen, height)
-end
-
-"""
- readcropparameters(generalcropfile, cropgrowthfile)
-
-Parse a CSV file containing the required parameter values for each crop
-(as produced from the original ALMaSS file by `convert_almass_data.py`).
-"""
-function readcropparameters(generalcropfile::String, growthfile::String)
- @debug "Reading crop parameters"
- cropdata = CSV.File(generalcropfile, missingstring="NA", dateformat="d U",
- types=[String,Date,Date,Date,Date,Float64])
- growthdata = CSV.File(growthfile, missingstring="NA",
- types=[Int,String,String,GrowthPhase,String,
- Float64,Float64,Float64,Float64])
- croptypes = Dict{String,CropType}()
- for crop in cropdata
- cropgrowthdata = growthdata |> filter(x -> !ismissing(x.crop_name) &&
- x.crop_name == crop.name)
- highnuts = buildgrowthcurve(cropgrowthdata |>
- filter(x -> x.nutrient_status=="high"))
- lownuts = buildgrowthcurve(cropgrowthdata |>
- filter(x -> x.nutrient_status=="low"))
- croptypes[crop.name] = CropType(crop.name, crop.minsowdate, crop.maxsowdate,
- crop.minharvestdate, crop.maxharvestdate,
- crop.mingrowthtemp, highnuts, lownuts)
- end
- croptypes
-end
diff --git a/src/crop/farmplot.jl b/src/crop/farmplot.jl
index 359c0c00d53ac50f5614f8cd71f8f3a7faf5ffad..c03d2c2d97f65265d9cfbcbc2f1fb5ad0bf45f9e 100644
--- a/src/crop/farmplot.jl
+++ b/src/crop/farmplot.jl
@@ -3,155 +3,48 @@
### This file contains code for the fields that farmers manage.
###
-#XXX not sure whether it makes sense to have this as an agent type,
-# or perhaps better a grid property?
-"""
- FarmPlot
-
-This represents one field, i.e. a collection of pixels with the same management.
-This is the spatial unit with which the crop growth model and the farm model work.
-"""
-mutable struct FarmPlot <: ModelAgent
- #TODO add Unitful
+mutable struct FarmPlot{T} <: ModelAgent
const id::Int64
pixels::Vector{Tuple{Int64, Int64}}
- croptype::CropType
- phase::GrowthPhase
- growingdegreedays::Float64
- height::Float64
- LAItotal::Float64
- LAIgreen::Float64
- #biomass::Float64 #XXX I need to figure out how to calculate this
- events::Vector{Management}
+ cropstate :: T
end
-"""
- initfields!(model)
-
-Initialise the model with its farm plots.
-"""
-function initfields!(model::SimulationModel)
- n = 0
- convertid = Dict{Int64,Int64}()
- width, height = size(model.landscape)
- for x in 1:width
- for y in 1:height
- # for each pixel, we need to extract the field ID given by the map input
- # file, and convert it into the internal object ID used by Agents.jl,
- # creating a new agent object if necessary
- rawid = model.landscape[x,y].fieldid
- (ismissing(rawid)) && continue
- if rawid in keys(convertid)
- objectid = convertid[rawid]
- model.landscape[x,y].fieldid = objectid
- push!(model.farmplots[objectid].pixels, (x,y))
- else
- #XXX does this phase calculation work out?
- month(model.date) < 3 ? phase = janfirst : phase = marchfirst
- fp = FarmPlot(length(model.farmplots)+1, [(x,y)],
- model.crops["natural grass"], phase,
- 0.0, 0.0, 0.0, 0.0, Vector{Management}())
- push!(model.farmplots, fp)
- model.landscape[x,y].fieldid = fp.id
- convertid[rawid] = fp.id
- n += 1
- end
- end
- end
- @info "Initialised $n farm plots."
-end
+croptype(f::FarmPlot{T}) where {T} = croptype(f.cropstate)
+cropname(f::FarmPlot{T}) where {T} = cropname(croptype(f))
+cropheight(f::FarmPlot{T}) where {T} = cropheight(f.cropstate)
+cropcover(f::FarmPlot{T}) where {T} = cropcover(f.cropstate)
+cropyield(f::FarmPlot{T}) where {T} = cropyield(f.cropstate)
"""
stepagent!(farmplot, model)
Update a farm plot by one day.
"""
-function stepagent!(farmplot::FarmPlot, model::SimulationModel)
- # update growing degree days
- # if no crop-specific base temperature is given, default to 5°C
- # (https://www.eea.europa.eu/publications/europes-changing-climate-hazards-1/heat-and-cold/heat-and-cold-2014-mean)
- basetemp = farmplot.croptype.mingrowthtemp
- ismissing(basetemp) && (basetemp = 5.0)
- gdd = (maxtemp(model)+mintemp(model))/2 - basetemp
- gdd > 0 && (farmplot.growingdegreedays += gdd)
- # update the phase on key dates
- monthday(model.date) == (1,1) && (farmplot.phase = janfirst)
- monthday(model.date) == (3,1) && (farmplot.phase = marchfirst)
- # update crop growth
- growcrop!(farmplot, model)
+function stepagent!(farmplot::FarmPlot{T}, model::SimulationModel) where T
+ stepagent!(farmplot.cropstate, model)
end
-
-## CROP MANAGEMENT AND GROWTH FUNCTIONS
-
"""
- sow!(cropname, farmplot, model)
+ sow!(farmplot, model, cropname)
-Sow the specified crop on this farmplot.
+Sow the specified crop on the farmplot.
"""
-function sow!(cropname::String, farmplot::FarmPlot, model::SimulationModel)
- createevent!(model, farmplot.pixels, sowing)
- farmplot.croptype = model.crops[cropname]
- farmplot.phase = sow
+function sow!(farmplot::FarmPlot, model::SimulationModel, cropname::String)
#XXX test if the crop is sowable?
+ createevent!(model, farmplot.pixels, sowing)
+ sow!(farmplot.cropstate, model, cropname)
end
"""
harvest!(farmplot, model)
-Harvest the crop on this farmplot.
+Harvest the crop of this farmplot.
"""
-function harvest!(farmplot::FarmPlot, model::SimulationModel)
+function harvest!(farmplot::FarmPlot{T}, model::SimulationModel) where T
createevent!(model, farmplot.pixels, harvesting)
- farmplot.phase in [harvest1, harvest2] ?
- farmplot.phase = harvest2 :
- farmplot.phase = harvest1
- # height & LAI will be automatically adjusted by the growth function
- #TODO calculate and return yield
+ harvest!(farmplot.cropstate, model) # TODO: multiply with area to return units of `g`
end
-#TODO fertilise!()
-#TODO spray!()
-#TODO till!()
-
-"""
- growcrop!(farmplot, model)
-
-Apply the relevant crop growth model to update the plants on this farm plot.
-Currently only supports the ALMaSS crop growth model by Topping et al.
-"""
-function growcrop!(farmplot::FarmPlot, model::SimulationModel)
- fertiliser in farmplot.events ?
- curve = farmplot.croptype.lownutrientgrowth :
- curve = farmplot.croptype.highnutrientgrowth
- points = curve.GDD[farmplot.phase]
- for p in 1:length(points)
- if points[p] == 99999
- return # the marker that there is no further growth this phase
- elseif points[p] == -1 # the marker to set all variables to specified values
- farmplot.height = curve.height[farmplot.phase][p]
- farmplot.LAItotal = curve.LAItotal[farmplot.phase][p]
- farmplot.LAIgreen = curve.LAIgreen[farmplot.phase][p]
- return
- else
- gdd = farmplot.growingdegreedays
- # figure out which is the correct slope value to use for growth
- if p == length(points) || gdd < points[p+1]
- farmplot.height += curve.height[farmplot.phase][p]
- farmplot.LAItotal += curve.LAItotal[farmplot.phase][p]
- farmplot.LAIgreen += curve.LAIgreen[farmplot.phase][p]
- return
- end
- #XXX To be precise, we ought to figure out if one or more inflection
- # points have been passed between yesterday and today, and calculate the
- # growth exactly up to the inflection point with the old slope, and onward
- # with the new slope. Not doing so will introduce a small amount of error,
- # although I think this is acceptable.
- end
- end
-end
-
-
## UTILITY FUNCTIONS
"""
@@ -175,7 +68,7 @@ Return the crop at this position, or nothing if there is no crop here (utility w
"""
function croptype(pos::Tuple{Int64,Int64}, model::SimulationModel)
ismissing(model.landscape[pos...].fieldid) ? nothing :
- model.farmplots[model.landscape[pos...].fieldid].croptype
+ croptype(model.farmplots[model.landscape[pos...].fieldid])
end
"""
@@ -186,29 +79,27 @@ Return the name of the crop at this position, or an empty string if there is no
"""
function cropname(pos::Tuple{Int64,Int64}, model::SimulationModel)
field = model.landscape[pos...].fieldid
- ismissing(field) ? "" : model.farmplots[field].croptype.name
+ ismissing(field) ? "" : cropname(model.farmplots[field])
end
"""
cropheight(model, position)
-Return the height of the crop at this position, or 0 if there is no crop here
+Return the height of the crop at this position, or nothing if there is no crop here
(utility wrapper).
"""
function cropheight(pos::Tuple{Int64,Int64}, model::SimulationModel)
- ismissing(model.landscape[pos...].fieldid) ? 0cm : #FIXME can I return something better than 0?
- model.farmplots[model.landscape[pos...].fieldid].height*1cm #FIXME units
+ ismissing(model.landscape[pos...].fieldid) ? 0.0cm : # TODO: 0.0cm correct here?
+ cropheight(model.farmplots[model.landscape[pos...].fieldid])
end
"""
cropcover(model, position)
-Return the percentage ground cover of the crop at this position, or 0 if there is no crop
-here (utility wrapper).
+Return the crop cover of the crop at this position, or nothing if
+there is no crop here (utility wrapper).
"""
function cropcover(pos::Tuple{Int64,Int64}, model::SimulationModel)
- #FIXME LAItotal != ground cover?
- ismissing(model.landscape[pos...].fieldid) ? 0 : #FIXME can I return something better than 0?
- model.farmplots[model.landscape[pos...].fieldid].LAItotal
+ ismissing(model.landscape[pos...].fieldid) ? 0.0 : # TODO: 0.0 correct here?
+ cropcover(model.farmplots[model.landscape[pos...].fieldid])
end
-
diff --git a/src/crop/simplecrop.jl b/src/crop/simplecrop.jl
new file mode 100644
index 0000000000000000000000000000000000000000..6bfe858364cf472af1a49b87cf0a8f455dfd5b8c
--- /dev/null
+++ b/src/crop/simplecrop.jl
@@ -0,0 +1,70 @@
+module SimpleCrop
+
+using Persefone:
+ FarmPlot,
+ Length,
+ cm,
+ SimulationModel
+
+import Persefone:
+ stepagent!,
+ croptype,
+ cropname,
+ cropheight,
+ cropcover,
+ cropyield
+
+struct CropType
+ name::String
+end
+
+cropname(ct::CropType) = ct.name
+
+mutable struct CropState
+ croptype::CropType
+ height::Length{Float64}
+end
+
+croptype(cs::CropState) = cs.croptype
+cropname(cs::CropState) = cropname(croptype(cs))
+cropheight(cs::CropState) = cs.height
+cropcover(cs::CropState) = 0.0
+cropyield(cs::CropState) = 0.0 # TODO: units?
+
+
+"""
+ stepagent!(farmplot, model)
+
+Update a crop state by one day.
+"""
+function stepagent!(cs::CropState, model::SimulationModel)
+ # height undergoes random diffusion, bounded by 0
+ delta_height = (2 * rand() - 1) * cm
+ cs.height = max(0.0cm, cs.height + delta_height)
+end
+
+"""
+ sow!(cropstate, model, cropname)
+
+Change the cropstate to sow the specified crop.
+"""
+function sow!(cs::CropState, model::SimulationModel, cropname::String)
+ cs.croptype = model.crops[cropname]
+ cs.height = 0cm
+end
+
+"""
+ harvest!(cropstate, model)
+
+Harvest the crop of this cropstate.
+"""
+function harvest!(cs::CropState, model::SimulationModel)
+ # TODO: set cs.croptype ?
+ # TODO: this 1.0g/cm/m^2 height_to_yield factor should be a param
+ # for every crop type
+ yield = cs.height * 1.0g/cm/m^2
+ cs.height = 0cm
+ return yield
+end
+
+end # module SimpleCrop
diff --git a/src/farm/farm.jl b/src/farm/farm.jl
index 391b825997502a7fc290c91e4b5b036b753a0adc..6d7b1588f33e54cb67e14a5fe86f4153e4379c15 100644
--- a/src/farm/farm.jl
+++ b/src/farm/farm.jl
@@ -12,8 +12,9 @@ mutable struct Farmer <: ModelAgent
#XXX make this into an abstract type and create subtypes for different
# farm submodels? (#69)
const id::Int64
- fields::Vector{FarmPlot}
- croprotation::Vector{CropType}
+ # TODO: hardcoded ALMaSS crop model
+ fields::Vector{FarmPlot{ALMaSS.CropState}}
+ croprotation::Vector{ALMaSS.CropType}
#TODO add AES
end
diff --git a/src/parameters.toml b/src/parameters.toml
index 2f2215b842291d7d0a5b65f6c7dc223cc7ee3299..659e8870ed158976062c91433ad7e06b69045427 100644
--- a/src/parameters.toml
+++ b/src/parameters.toml
@@ -40,7 +40,7 @@ indoutfreq = "end" # output frequency individual-level data, daily/monthly/yearl
insectmodel = ["season", "habitat", "pesticides", "weather"] # factors affecting insect growth
[crop]
-cropmodel = "almass" # crop growth model to use, "almass" or "aquacrop"
+cropmodel = "almass" # crop growth model to use, "almass", "aquacrop", or "simple"
cropfile = "data/crops/almass/crop_data_general.csv" # file with general crop parameters
growthfile = "data/crops/almass/almass_crop_growth_curves.csv" # file with crop growth parameters
diff --git a/test/crop_tests.jl b/test/crop_tests.jl
index 2ab0affb0eaad04c1a71e141c06c49e497d21185..91469250561ef844d1fd99969b6442701c8133e2 100644
--- a/test/crop_tests.jl
+++ b/test/crop_tests.jl
@@ -1,6 +1,45 @@
### Persefone.jl - a model of agricultural landscapes and ecosystems in Europe.
###
-### These are the tests for the crop growth model.
+### These are the tests for the crop growth models.
###
-#TODO
+const TESTPARAM_ALMASS = joinpath(pkgdir(Persefone), "test", "test_parameters_almass.toml")
+const TESTPARAM_SIMPLECROP = joinpath(pkgdir(Persefone), "test", "test_parameters_simplecrop.toml")
+
+@testset for paramfile in (TESTPARAM_ALMASS, TESTPARAM_SIMPLECROP)
+ @testset "Model initialisation" begin
+ model = initialise(paramfile)
+ @test model isa AgricultureModel
+ end
+ @testset "Time step" begin
+ model = initialise(paramfile)
+ stepsimulation!(model)
+ @test model isa AgricultureModel
+ end
+end
+
+@testset "Submodule ALMaSS" begin
+ ct = Ps.ALMaSS.CropType("olive tree", missing, missing, missing, missing, missing,
+ missing, missing)
+ fp = FarmPlot(0, [(0,0)],
+ Ps.ALMaSS.CropState(ct, Ps.ALMaSS.janfirst, 0.0, 0.0m, 0.0, 0.0, Ps.Management[]))
+ @test fp isa FarmPlot
+ @test fp isa FarmPlot{Ps.ALMaSS.CropState}
+ @test croptype(fp) isa Ps.ALMaSS.CropType
+ @test cropname(fp) isa String
+ @test cropheight(fp) isa Length{Float64}
+ @test cropcover(fp) isa Float64
+ @test cropyield(fp) isa Float64
+end
+
+@testset "Submodule SimpleCrop" begin
+ ct = Ps.SimpleCrop.CropType("olive tree")
+ fp = FarmPlot(0, [(0,0)], Ps.SimpleCrop.CropState(ct, 0.0cm))
+ @test fp isa FarmPlot
+ @test fp isa FarmPlot{Ps.SimpleCrop.CropState}
+ @test croptype(fp) isa Ps.SimpleCrop.CropType
+ @test cropname(fp) isa String
+ @test cropheight(fp) isa Length{Float64}
+ @test cropcover(fp) isa Float64
+ @test cropyield(fp) isa Float64
+end
diff --git a/test/nature_tests.jl b/test/nature_tests.jl
index afeab5ea92b0168c134ade932923350d7ec47bef..a5ef3f91df4d2c9dea80527dc129fda8615a9f09 100644
--- a/test/nature_tests.jl
+++ b/test/nature_tests.jl
@@ -47,17 +47,22 @@ end) # end eval
@testset "Habitat macros" begin
# set up the testing landscape
model = inittestmodel()
- model.landscape[6,6] = Pixel(Ps.agriculture, 1)
- push!(model.farmplots,
- FarmPlot(1, [(6,6)], model.crops["winter wheat"], Ps.janfirst,
- 0.0, 0.0, 0.0, 0.0, Vector{Ps.Management}()))
+ model.landscape[6,6] = Pixel(Ps.agriculture, 1, [], [], [])
+ fp = Ps.FarmPlot(
+ 1, [(6,6)],
+ Ps.ALMaSS.CropState(
+ model.crops["winter wheat"], Ps.ALMaSS.janfirst,
+ 0.0, 0.0m, 0.0, 0.0, Vector{Ps.Management}()
+ )
+ )
+ push!(model.farmplots, fp)
push!(model.animals,
Ps.Mermaid(1, Ps.male, (-1,-1), (3,3), Ps.life),
Ps.Mermaid(2, Ps.female, (-1,-1), (4,4), Ps.life))
# create a set of habitat descriptors
h1 = @habitat(@landcover() == Ps.water)
h2 = @habitat(@cropname() == "winter wheat" &&
- @cropheight() < 2) #XXX randomly chosen value
+ @cropheight() < 2m) #XXX randomly chosen value
h3 = @habitat(@distanceto(Ps.water) > 20m &&
@distancetoedge() <= 20m)
#FIXME nested macros don't work properly, counting seems wrong
diff --git a/test/paramscan.toml b/test/paramscan.toml
index d2db42fb9add8f360c2247adb5cc3e65f8223054..565e6855f876dbc6474137cd22b0353c384485f8 100644
--- a/test/paramscan.toml
+++ b/test/paramscan.toml
@@ -28,5 +28,5 @@ popoutfreq = "daily" # output frequency population-level data, daily/monthly/yea
indoutfreq = "end" # output frequency individual-level data, daily/monthly/yearly/end/never
[crop]
-cropmodel = "linear" # crop growth model to use, "linear" or "aquacrop" (not yet implemented)
+cropmodel = "simple" # crop growth model to use, "simple" or "almass"
diff --git a/test/runtests.jl b/test/runtests.jl
index b36a792bd36b7795f1efbe3a8603a2f5c41c63e0..4aaf78eb80ee893ba409eed6b396547f056721ab 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -43,22 +43,27 @@ function inittestmodel(smallmap=true)
TESTSETTINGS["world.weatherfile"]),
TESTSETTINGS["core.startdate"],
TESTSETTINGS["core.enddate"])
- crops = Ps.readcropparameters(TESTSETTINGS["crop.cropfile"],
- TESTSETTINGS["crop.growthfile"])
- model = AgricultureModel(TESTSETTINGS,
- StableRNG(TESTSETTINGS["core.seed"]),
- global_logger(),
- Vector{DataOutput}(),
- Dict{String, DataFrame}(),
- TESTSETTINGS["core.startdate"],
- landscape,
- weather,
- crops,
- Vector{Farmer}(),
- Vector{FarmPlot}(),
- Vector{Union{Animal,Nothing}}(),
- Vector{Pair{Animal, Date}}(),
- Vector{FarmEvent}())
+ # TODO: support other crop models besides ALMaSS
+ crops = Ps.ALMaSS.readcropparameters(TESTSETTINGS["crop.cropfile"],
+ TESTSETTINGS["crop.growthfile"])
+ farmers = Vector{Farmer}()
+ farmplots = Vector{Ps.FarmPlot{Ps.ALMaSS.CropState}}()
+ model = AgricultureModel{Ps.ALMaSS.CropType,Ps.ALMaSS.CropState}(
+ TESTSETTINGS,
+ StableRNG(TESTSETTINGS["core.seed"]),
+ global_logger(),
+ Vector{DataOutput}(),
+ Dict{String, DataFrame}(),
+ TESTSETTINGS["core.startdate"],
+ landscape,
+ weather,
+ crops,
+ farmers,
+ farmplots,
+ Vector{Union{Animal,Nothing}}(),
+ Vector{Pair{Animal, Date}}(),
+ Vector{FarmEvent}()
+ )
model
end
diff --git a/test/test_parameters.toml b/test/test_parameters.toml
index 4d88eb3bf0deeee889c106ce4ff90329e1706a89..94d6a9a8b25be61910d3a7d690079d2405256295 100644
--- a/test/test_parameters.toml
+++ b/test/test_parameters.toml
@@ -38,7 +38,7 @@ 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 = "almass" # crop growth model to use, "almass" or "aquacrop"
+cropmodel = "almass" # crop growth model to use, "almass", "aquacrop", or "simple"
cropfile = "crop_data_general.csv" # file with general crop parameters
growthfile = "almass_crop_growth_curves.csv" # file with crop growth parameters
diff --git a/test/test_parameters_almass.toml b/test/test_parameters_almass.toml
new file mode 100644
index 0000000000000000000000000000000000000000..1308029ee9b7b387064752b0c0e3a70174798c77
--- /dev/null
+++ b/test/test_parameters_almass.toml
@@ -0,0 +1,44 @@
+### Persefone.jl - a model of agricultural landscapes and ecosystems in Europe.
+###
+### This configuration file is used for the test suite.
+###
+
+#XXX remember that changes here may break tests!
+
+[core]
+configfile = "test_parameters.toml" # location of the configuration file
+outdir = "results_testsuite" # location and name of the output folder
+logoutput = "both"
+overwrite = true # overwrite the output directory? (true/false/"ask")
+csvoutput = true # save collected data in CSV files
+visualise = true # generate result graphs
+storedata = true # keep collected data in memory
+loglevel = "warn" # verbosity level: "debug", "info", "warn"
+processors = 6 # number of processors to use on parallel runs
+seed = 1 # seed value for the RNG (0 -> random value)
+# dates to start and end the simulation
+startdate = 2022-02-01
+enddate = 2022-03-31
+
+[world]
+mapdirectory = "." # the directory in which all geographic data are stored
+mapresolution = 10 # map resolution in meters
+landcovermap = "landcover_jena.tif" # location of the landcover map
+farmfieldsmap = "fields_jena.tif" # location of the field geometry map
+weatherfile = "weather_jena.csv" # location of the weather data file
+
+[farm]
+farmmodel = "FieldManager" # which version of the farm model to use (not yet implemented)
+
+[nature]
+targetspecies = ["Wolpertinger", "Wyvern"] # list of target species to simulate - example species
+#targetspecies = ["Skylark"] # list of target species to simulate
+popoutfreq = "daily" # output frequency population-level data, daily/monthly/yearly/end/never
+indoutfreq = "daily" # output frequency individual-level data, daily/monthly/yearly/end/never
+insectmodel = ["season", "habitat", "pesticides"] # which factors affect insect growth ("weather" is not yet implemented)
+
+[crop]
+cropmodel = "almass" # crop growth model to use, "almass", "aquacrop", or "simple"
+cropfile = "crop_data_general.csv" # file with general crop parameters
+growthfile = "almass_crop_growth_curves.csv" # file with crop growth parameters
+
diff --git a/test/test_parameters_simplecrop.toml b/test/test_parameters_simplecrop.toml
new file mode 100644
index 0000000000000000000000000000000000000000..fc199a8c8664b0d900e131d9b5b159dc006ee502
--- /dev/null
+++ b/test/test_parameters_simplecrop.toml
@@ -0,0 +1,45 @@
+### Persefone.jl - a model of agricultural landscapes and ecosystems in Europe.
+###
+### This configuration file is used for the test suite.
+###
+
+#XXX remember that changes here may break tests!
+
+[core]
+configfile = "test_parameters.toml" # location of the configuration file
+outdir = "results_testsuite" # location and name of the output folder
+logoutput = "both"
+overwrite = true # overwrite the output directory? (true/false/"ask")
+csvoutput = true # save collected data in CSV files
+visualise = true # generate result graphs
+storedata = true # keep collected data in memory
+loglevel = "warn" # verbosity level: "debug", "info", "warn"
+processors = 6 # number of processors to use on parallel runs
+seed = 1 # seed value for the RNG (0 -> random value)
+# dates to start and end the simulation
+startdate = 2022-02-01
+enddate = 2022-03-31
+
+[world]
+mapdirectory = "." # the directory in which all geographic data are stored
+mapresolution = 10 # map resolution in meters
+landcovermap = "landcover_jena.tif" # location of the landcover map
+farmfieldsmap = "fields_jena.tif" # location of the field geometry map
+weatherfile = "weather_jena.csv" # location of the weather data file
+
+[farm]
+farmmodel = "FieldManager" # which version of the farm model to use (not yet implemented)
+
+[nature]
+targetspecies = ["Wolpertinger", "Wyvern"] # list of target species to simulate - example species
+#targetspecies = ["Skylark"] # list of target species to simulate
+popoutfreq = "daily" # output frequency population-level data, daily/monthly/yearly/end/never
+indoutfreq = "daily" # output frequency individual-level data, daily/monthly/yearly/end/never
+insectmodel = ["season", "habitat", "pesticides"] # which factors affect insect growth ("weather" is not yet implemented)
+
+[crop]
+cropmodel = "simple" # crop growth model to use, "almass", "aquacrop", or "simple"
+cropfile = "crop_data_general.csv" # file with general crop parameters
+# TODO: this parameter should be unnecessary, but removing it
+# currently gives an error
+growthfile = "almass_crop_growth_curves.csv" # file with crop growth parameters