From 580b6ae73e9905b4dddbdc2aa63f1ddbc06a9241 Mon Sep 17 00:00:00 2001
From: Marco Matthies <71844+marcom@users.noreply.github.com>
Date: Mon, 15 Jul 2024 10:10:38 +0200
Subject: [PATCH] Fix bugs from moving ALMaSS crop model into submodule
---
src/Persefone.jl | 4 +-
src/core/simulation.jl | 6 +-
src/crop/crops.jl | 203 +++++++++++++++++++++++++++++++++++++++-
src/farm/farm.jl | 2 +-
src/nature/macros.jl | 4 +-
test/landscape_tests.jl | 4 +-
test/runtests.jl | 2 +-
7 files changed, 214 insertions(+), 11 deletions(-)
diff --git a/src/Persefone.jl b/src/Persefone.jl
index ec269bf..8b55870 100644
--- a/src/Persefone.jl
+++ b/src/Persefone.jl
@@ -46,7 +46,7 @@ export
Weather,
FarmEvent,
ModelAgent,
- FarmPlot,
+ # FarmPlot,
Animal,
Farmer,
DataOutput,
@@ -127,6 +127,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")
diff --git a/src/core/simulation.jl b/src/core/simulation.jl
index b94b521..f745b24 100644
--- a/src/core/simulation.jl
+++ b/src/core/simulation.jl
@@ -23,7 +23,7 @@ mutable struct AgricultureModel <: SimulationModel
weather::Dict{Date,Weather}
crops::Dict{String,ALMaSS.CropType}
farmers::Vector{Farmer}
- farmplots::Vector{FarmPlot}
+ farmplots::Vector{ALMaSS.FarmPlot}
animals::Vector{Union{Animal,Nothing}}
migrants::Vector{Pair{Animal,Date}}
events::Vector{FarmEvent}
@@ -129,12 +129,12 @@ function initmodel(settings::Dict{String, Any})
weather,
crops,
Vector{Farmer}(),
- Vector{FarmPlot}(),
+ Vector{ALMaSS.FarmPlot}(),
Vector{Union{Animal,Nothing}}(),
Vector{Pair{Animal, Date}}(),
Vector{FarmEvent}())
saveinputfiles(model)
- initfields!(model)
+ ALMaSS.initfields!(model)
initfarms!(model)
initnature!(model)
model
diff --git a/src/crop/crops.jl b/src/crop/crops.jl
index 3b4ce0c..e66ce77 100644
--- a/src/crop/crops.jl
+++ b/src/crop/crops.jl
@@ -6,7 +6,10 @@
#TODO write tests for input functions
module ALMaSS
-using Dates: Date
+
+using Persefone: ModelAgent, EventType, SimulationModel, maxtemp, mintemp, fertiliser
+import Persefone: stepagent!
+using Dates: Date, month, monthday
using CSV: CSV
"""
@@ -124,4 +127,202 @@ function readcropparameters(generalcropfile::String, growthfile::String)
croptypes
end
+
+#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
+ const id::Int64
+ pixels::Vector{Tuple{Int64, Int64}}
+ croptype::ALMaSS.CropType
+ phase::ALMaSS.GrowthPhase
+ growingdegreedays::Float64
+ height::Float64
+ LAItotal::Float64
+ LAIgreen::Float64
+ #biomass::Float64 #XXX I need to figure out how to calculate this
+ events::Vector{EventType}
+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 = ALMaSS.janfirst : phase = ALMaSS.marchfirst
+ fp = FarmPlot(length(model.farmplots)+1, [(x,y)],
+ model.crops["natural grass"], phase,
+ 0.0, 0.0, 0.0, 0.0, Vector{EventType}())
+ 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
+
+"""
+ 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 = ALMaSS.janfirst)
+ monthday(model.date) == (3,1) && (farmplot.phase = ALMaSS.marchfirst)
+ # update crop growth
+ growcrop!(farmplot, model)
+end
+
+
+## CROP MANAGEMENT AND GROWTH FUNCTIONS
+
+"""
+ sow!(cropname, farmplot, model)
+
+Sow the specified crop on this farmplot.
+"""
+function sow!(cropname::String, farmplot::FarmPlot, model::SimulationModel)
+ createevent!(model, farmplot.pixels, sowing)
+ farmplot.croptype = model.crops[cropname]
+ farmplot.phase = ALMaSS.sow
+ #XXX test if the crop is sowable?
+end
+
+"""
+ harvest!(farmplot, model)
+
+Harvest the crop on this farmplot.
+"""
+function harvest!(farmplot::FarmPlot, model::SimulationModel)
+ createevent!(model, farmplot.pixels, harvesting)
+ farmplot.phase in [ALMaSS.harvest1, ALMaSS.harvest2] ?
+ farmplot.phase = ALMaSS.harvest2 :
+ farmplot.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!(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
+
+"""
+ averagefieldsize(model)
+
+Calculate the average field size in hectares for the model landscape.
+"""
+function averagefieldsize(model::SimulationModel)
+ conversionfactor = 100 #our pixels are currently 10x10m, so 100 pixels per hectare
+ sizes::Vector{Float64} = []
+ for fp in model.farmplots
+ push!(sizes, size(fp.pixels)[1]/conversionfactor)
+ end
+ round(sum(sizes)/size(sizes)[1], digits=2)
+end
+
+"""
+ croptype(model, position)
+
+Return the crop at this position, or nothing if there is no crop here (utility wrapper).
+"""
+function croptype(pos::Tuple{Int64,Int64}, model::SimulationModel)
+ ismissing(model.landscape[pos...].fieldid) ? nothing :
+ model.farmplots[model.landscape[pos...].fieldid].croptype
+end
+
+"""
+ cropname(model, position)
+
+Return the name of the crop at this position, or an empty string if there is no crop here
+(utility wrapper).
+"""
+function cropname(pos::Tuple{Int64,Int64}, model::SimulationModel)
+ field = model.landscape[pos...].fieldid
+ ismissing(field) ? "" : model.farmplots[field].croptype.name
+end
+
+"""
+ cropheight(model, position)
+
+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) ? nothing :
+ model.farmplots[model.landscape[pos...].fieldid].height
+end
+
end # module ALMaSS
diff --git a/src/farm/farm.jl b/src/farm/farm.jl
index ef3d34c..e9e70c3 100644
--- a/src/farm/farm.jl
+++ b/src/farm/farm.jl
@@ -12,7 +12,7 @@ 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}
+ fields::Vector{ALMaSS.FarmPlot}
croprotation::Vector{ALMaSS.CropType}
#TODO add AES
end
diff --git a/src/nature/macros.jl b/src/nature/macros.jl
index add747c..b9345b8 100644
--- a/src/nature/macros.jl
+++ b/src/nature/macros.jl
@@ -342,7 +342,7 @@ This is a utility wrapper that can only be used nested within [`@phase`](@ref)
or [`@habitat`](@ref).
"""
macro cropname()
- :(cropname($(esc(:pos)), $(esc(:model))))
+ :(ALMaSS.cropname($(esc(:pos)), $(esc(:model))))
end
"""
@@ -353,7 +353,7 @@ This is a utility wrapper that can only be used nested within [`@phase`](@ref)
or [`@habitat`](@ref).
"""
macro cropheight()
- :(cropheight($(esc(:pos)), $(esc(:model))))
+ :(ALMaSS.cropheight($(esc(:pos)), $(esc(:model))))
end
"""
diff --git a/test/landscape_tests.jl b/test/landscape_tests.jl
index b825765..9fcc4ec 100644
--- a/test/landscape_tests.jl
+++ b/test/landscape_tests.jl
@@ -6,7 +6,7 @@
@testset "Landscape initialisation" begin
model = inittestmodel(false)
# these tests are specific to the Jena maps
- @test_logs (:info, "Initialised 2092 farm plots.") match_mode=:any Ps.initfields!(model)
+ @test_logs (:info, "Initialised 2092 farm plots.") match_mode=:any Ps.ALMaSS.initfields!(model)
@test size(model.landscape) == (1754, 1602)
@test Ps.landcover((100,100), model) == Ps.forest
@test Ps.landcover((300,1), model) == Ps.soil
@@ -14,7 +14,7 @@
@test Ps.landcover((400,400), model) == Ps.grass
@test Ps.landcover((800,800), model) == Ps.agriculture
@test Ps.landcover((1100,1100), model) == Ps.builtup
- @test Ps.averagefieldsize(model) == 5.37
+ @test Ps.ALMaSS.averagefieldsize(model) == 5.37
@test count(f -> ismissing(f.fieldid), model.landscape) == 1685573
@test length(Ps.farmplot((800,800), model).pixels) == 4049
end
diff --git a/test/runtests.jl b/test/runtests.jl
index 39aef77..b318010 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -55,7 +55,7 @@ function inittestmodel(smallmap=true)
weather,
crops,
Vector{Farmer}(),
- Vector{FarmPlot}(),
+ Vector{Ps.ALMaSS.FarmPlot}(),
Vector{Union{Animal,Nothing}}(),
Vector{Pair{Animal, Date}}(),
Vector{FarmEvent}())
--
GitLab