diff --git a/src/analysis/makieplots.jl b/src/analysis/makieplots.jl
index 30cacefedaa00499f0999144024df95f5aad39da..70ba8ce93fbeab50d95d77cec4c3be01d910ebe5 100644
--- a/src/analysis/makieplots.jl
+++ b/src/analysis/makieplots.jl
@@ -59,7 +59,7 @@ function populationtrends(model::SimulationModel)
update_theme!(palette=(color=cgrad(:seaborn_bright, ncolors),), cycle=[:color])
f = Figure()
dates = @param(core.startdate):@param(core.enddate)
- axlimits = (1, length(dates), 0, maximum(pops[!,:Abundance]))
+ axlimits = (1, length(dates), 0, maximum(pops[!,:Abundance])*1.1)
ax = Axis(f[1,1], xlabel="Date", ylabel="Population size",
limits=axlimits, xticks=datetickmarks(dates))
for s in @param(nature.targetspecies)
@@ -133,10 +133,10 @@ Returns a Makie figure object.
"""
function croptrends(model::SimulationModel)
cropdata = @data("fields")
- f = Figure(size=(1000,1000))
+ f = Figure(size=(1200,1000))
dates = @param(core.startdate):@param(core.enddate)
# plot cropped area over time
- axlimitsarea = (1, length(dates), 0, maximum(cropdata[!,:Area])*1.2)
+ axlimitsarea = (1, length(dates), 0, maximum(cropdata[!,:Area])*1.1)
ax1 = Axis(f[1,1], xlabel="Date", ylabel="Cropped area (ha)",
limits=axlimitsarea, xticks=datetickmarks(dates))
for c in unique(cropdata.Crop)
@@ -144,16 +144,17 @@ function croptrends(model::SimulationModel)
(iszero(size(points)[1]) || iszero(sum(points.Area))) && continue
lines!(f[1,1], Vector{Float32}(points.Area), linewidth=3, label=c)
end
- axislegend("Crop"; position=:lt)
+ axislegend("Crop"; position=:rt)
# plot average crop height over time
- axlimitsheight = (1, length(dates), 0, maximum(cropdata[!,:Height])*1.2)
+ axlimitsheight = (1, length(dates), 0, maximum(cropdata[!,:Height])*1.1)
ax2 = Axis(f[2,1], xlabel="Date", ylabel="Average plant height (cm)",
limits=axlimitsheight, xticks=datetickmarks(dates))
for c in unique(cropdata.Crop)
points = @select!(@subset(cropdata, :Crop .== c), :Height)
- (iszero(size(points)[1]) || iszero(sum(points.Height))) && continue
+ (iszero(size(points)[1]) || iszero(sum(points.Height))) && c != "no growth" && continue
lines!(f[2,1], Vector{Float32}(points.Height), linewidth=3, label=c)
end
+ axislegend("Crop"; position=:rt)
f
end
diff --git a/src/crop/almass.jl b/src/crop/almass.jl
index c6fc177748a9bd30ea914f6aaf171d0dcb98750c..891cbcf828276423c4f55a14dbefa9b4dbd5590b 100644
--- a/src/crop/almass.jl
+++ b/src/crop/almass.jl
@@ -3,7 +3,11 @@
### This file is responsible for managing the crop growth modules.
###
-#TODO write tests for input functions
+#FIXME Looking at fields.pdf in the Persefone output, crop maturity is not recognised
+# properly, so most crops (except maize) are never planted. At the same time, maize never
+# seems to start growing.
+
+#TODO write tests to compare our output to the output of the original ALMaSS algorithm
module ALMaSS
@@ -26,7 +30,8 @@ import Persefone:
cropyield,
sow!,
harvest!,
- isharvestable
+ isharvestable,
+ bounds
using Dates: Date, month, monthday
using CSV: CSV
@@ -92,7 +97,7 @@ mutable struct CropState
LAIgreen::Float64
#biomass::Float64 #XXX I need to figure out how to calculate this
mature::Bool #TODO how do we determine this?
- events::Vector{Management}
+ events::Vector{Management} #FIXME does this do anything?
end
croptype(cs::CropState) = cs.croptype
@@ -180,6 +185,14 @@ end
Update a farm plot by one day.
"""
function stepagent!(cs::CropState, model::SimulationModel)
+ # update the phase on key dates
+ if monthday(model.date) == (1, 1)
+ cs.phase = ALMaSS.janfirst
+ cs.growingdegreedays = 0.0
+ elseif monthday(model.date) == (3, 1)
+ cs.phase = ALMaSS.marchfirst
+ cs.growingdegreedays = 0.0
+ end
# 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)
@@ -187,11 +200,8 @@ function stepagent!(cs::CropState, model::SimulationModel)
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)
+ growcrop!(cs, gdd, model)
end
@@ -203,9 +213,16 @@ end
Change the cropstate to sow the specified crop.
"""
function sow!(cs::CropState, model::SimulationModel, cropname::String)
- #XXX test if the crop is sowable?
+ !ismissing(cs.croptype.minsowdate) && model.date < cs.croptype.minsowdate &&
+ @warn "$(model.date) is earlier than the minimum sowing date for $(cropname)."
cs.croptype = model.crops[cropname]
cs.phase = ALMaSS.sow
+ cs.growingdegreedays = 0.0
+ cs.height = 0.0cm
+ cs.LAItotal = 0.0
+ cs.LAIgreen = 0.0
+ cs.mature = false
+ cs.events = Vector{Management}()
end
"""
@@ -217,6 +234,8 @@ function harvest!(cs::CropState, model::SimulationModel)
cs.phase in [ALMaSS.harvest1, ALMaSS.harvest2] ?
cs.phase = ALMaSS.harvest2 :
cs.phase = ALMaSS.harvest1
+ cs.growingdegreedays = 0.0
+ cs.mature = false
# height & LAI will be automatically adjusted by the growth function
#TODO calculate and return yield
end
@@ -226,17 +245,18 @@ end
#TODO till!()
"""
- growcrop!(cropstate, model)
+ growcrop!(cropstate, gdd, 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.
+et al. (see `ALMaSS/Landscape/plant.cpp:PlantGrowthData::FindDiff()` and
+`ALMaSS/Landscape/elements.cpp:VegElement::DoDevelopment()`).
"""
-function growcrop!(cs::CropState, model::SimulationModel)
+function growcrop!(cs::CropState, gdd::Float64, model::SimulationModel)
fertiliser in cs.events ?
curve = cs.croptype.lownutrientgrowth :
curve = cs.croptype.highnutrientgrowth
- points = curve.GDD[cs.phase]
+ points = curve.GDD[cs.phase] #FIXME what if the curve is empty?
for p in 1:length(points)
if points[p] == 99999
!(cs.phase in (janfirst, sow)) && (cs.mature = true) #FIXME only in the last phase?
@@ -250,9 +270,14 @@ function growcrop!(cs::CropState, model::SimulationModel)
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]
+ #FIXME it appears to me from the ALMaSS source code that the curve value
+ # for the day is multiplied with that day's growingdegreedays to give the
+ # diff, which is then added to the previous values. However, this results
+ # in values that are staggeringly too large. I'm not quite sure what the
+ # issue is here...
+ cs.height += bounds(curve.height[cs.phase][p]*gdd)
+ cs.LAItotal += bounds(curve.LAItotal[cs.phase][p]*gdd)
+ cs.LAIgreen += bounds(curve.LAIgreen[cs.phase][p]*gdd)
return
end
#XXX To be precise, we ought to figure out if one or more inflection
diff --git a/src/crop/cropmodels.jl b/src/crop/cropmodels.jl
index aef409f0f531326deaa706dcb58b8303bf713c71..22d6b1470849aa46b922d82ec930e97018eadfec 100644
--- a/src/crop/cropmodels.jl
+++ b/src/crop/cropmodels.jl
@@ -66,7 +66,7 @@ function makecropstate(model::SimulationModel)
cs = ALMaSS.CropState(
model.crops["natural grass"],
phase,
- 0.0, 0.0m, 0.0, 0.0, false, Vector{Management}()
+ 0.0, 0.0cm, 0.0, 0.0, false, Vector{Management}()
)
elseif @param(crop.cropmodel) == "simple"
cs = SimpleCrop.CropState(
diff --git a/src/parameters.toml b/src/parameters.toml
index 9a5920f8b3c84e0e4ae9c7d25e7aea14b18062f4..484c1a9b02d4245f7e48b40ccfdf9de1f5221895 100644
--- a/src/parameters.toml
+++ b/src/parameters.toml
@@ -15,12 +15,12 @@ csvoutput = true # save collected data in CSV files
visualise = true # generate result graphs
storedata = true # keep collected data in memory
figureformat = "pdf" # file format to use for graphical output
-loglevel = "info" # verbosity level: "debug", "info", "warn"
+loglevel = "debug" # verbosity level: "debug", "info", "warn"
processors = 2 # number of processors to use on parallel runs
seed = 2 # seed value for the RNG (0 -> random value)
-startdate = 2022-01-01 # first day of the simulation
-enddate = 2022-12-31 # last day of the simulation
-#enddate = 2022-03-31 # last day of the simulation (test value)
+startdate = 2021-01-01 # first day of the simulation
+enddate = 2023-12-31 # last day of the simulation
+#enddate = 2022-01-02 # last day of the simulation (test value)
[world]
mapdirectory = "data/regions/jena-small" # the directory in which all geographic data are stored
@@ -30,14 +30,18 @@ farmfieldsmap = "fields.tif" # name of the field geometry map in the map directo
weatherfile = "weather.csv" # name of the weather data file in the map directory
[farm]
-farmmodel = "FieldManager" # which version of the farm model to use (not yet implemented)
+farmmodel = "BasicFarmer" # which version of the farm model to use
setaside = 0.04 # proportion of farm area set aside as fallow
+fieldoutfreq = "daily" # output frequency for crop/field data, daily/monthly/yearly/end/never
[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 = "end" # output frequency individual-level data, daily/monthly/yearly/end/never
+targetspecies = [] # XXX disable all species for farm model testing
+#targetspecies = ["Skylark"] # list of target species to simulate
+popoutfreq = "never" # output frequency population-level data, daily/monthly/yearly/end/never
+indoutfreq = "never" # output frequency individual-level data, daily/monthly/yearly/end/never
+#popoutfreq = "daily" # output frequency population-level data, daily/monthly/yearly/end/never
+#indoutfreq = "end" # output frequency individual-level data, daily/monthly/yearly/end/never
insectmodel = ["season", "habitat", "pesticides", "weather"] # factors affecting insect growth
[crop]