From 74d2d078f626e5bbba16ca142d3903257d646c1f Mon Sep 17 00:00:00 2001
From: Daniel Vedder <daniel.vedder@idiv.de>
Date: Fri, 9 Aug 2024 12:01:44 +0200
Subject: [PATCH] Added data recording and visualisation for the farm model
---
src/Persefone.jl | 1 +
src/analysis/makieplots.jl | 42 ++++++++++++++++++++++++++++++++------
src/crop/almass.jl | 2 +-
src/farm/farm.jl | 6 +-----
src/farm/farmdata.jl | 41 +++++++++++++++++++++++++++++++++++++
5 files changed, 80 insertions(+), 12 deletions(-)
create mode 100644 src/farm/farmdata.jl
diff --git a/src/Persefone.jl b/src/Persefone.jl
index bff3e22..f5aa7b4 100644
--- a/src/Persefone.jl
+++ b/src/Persefone.jl
@@ -149,6 +149,7 @@ include("crop/almass.jl")
include("crop/simplecrop.jl")
include("farm/farm.jl")
+include("farm/farmdata.jl")
include("nature/insects.jl")
include("nature/energy.jl")
diff --git a/src/analysis/makieplots.jl b/src/analysis/makieplots.jl
index fd92184..30cacef 100644
--- a/src/analysis/makieplots.jl
+++ b/src/analysis/makieplots.jl
@@ -61,7 +61,7 @@ function populationtrends(model::SimulationModel)
dates = @param(core.startdate):@param(core.enddate)
axlimits = (1, length(dates), 0, maximum(pops[!,:Abundance]))
ax = Axis(f[1,1], xlabel="Date", ylabel="Population size",
- limits=axlimits, xticks=gettickmarks(dates))
+ limits=axlimits, xticks=datetickmarks(dates))
for s in @param(nature.targetspecies)
points = @select!(@subset(pops, :Species .== s), :Abundance)
iszero(size(points)[1]) && continue
@@ -84,7 +84,7 @@ function skylarkpopulation(model::SimulationModel)
dates = @param(core.startdate):@param(core.enddate)
axlimits = (1, length(dates), 0, maximum(pops[!,:TotalAbundance]))
ax = Axis(f[1,1], xlabel="Date", ylabel="Population size",
- limits=axlimits, xticks=gettickmarks(dates))
+ limits=axlimits, xticks=datetickmarks(dates))
lines!(f[1,1], Vector{Float32}(pops.TotalAbundance), linewidth=3, label="Total abundance")
lines!(f[1,1], Vector{Float32}(pops.Mating), linewidth=3, label="Mating")
lines!(f[1,1], Vector{Float32}(pops.Breeding), linewidth=3, label="Breeding")
@@ -126,12 +126,44 @@ function skylarkstats(model::SimulationModel)
end
"""
- gettickmarks(dates)
+ croptrends(model)
+
+Plot a dual line graph of cropped area and average plant height per crop over time.
+Returns a Makie figure object.
+"""
+function croptrends(model::SimulationModel)
+ cropdata = @data("fields")
+ f = Figure(size=(1000,1000))
+ dates = @param(core.startdate):@param(core.enddate)
+ # plot cropped area over time
+ axlimitsarea = (1, length(dates), 0, maximum(cropdata[!,:Area])*1.2)
+ ax1 = Axis(f[1,1], xlabel="Date", ylabel="Cropped area (ha)",
+ limits=axlimitsarea, xticks=datetickmarks(dates))
+ for c in unique(cropdata.Crop)
+ points = @select!(@subset(cropdata, :Crop .== c), :Area)
+ (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)
+ # plot average crop height over time
+ axlimitsheight = (1, length(dates), 0, maximum(cropdata[!,:Height])*1.2)
+ 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
+ lines!(f[2,1], Vector{Float32}(points.Height), linewidth=3, label=c)
+ end
+ f
+end
+
+"""
+ datetickmarks(dates)
Given a vector of dates, construct a selection to use as tick mark locations.
Helper function for `[populationtrends](@ref)`
"""
-function gettickmarks(dates)
+function datetickmarks(dates)
ticks = (Int[], String[])
for i in 1:length(dates)
if Day(dates[i]) == Day(1)
@@ -147,5 +179,3 @@ function gettickmarks(dates)
end
#XXX add animation with record()? https://docs.makie.org/stable/explanations/animation/
-
-#TODO causes of mortality, nesting habitats, territory sizes
diff --git a/src/crop/almass.jl b/src/crop/almass.jl
index 6c020cc..c6fc177 100644
--- a/src/crop/almass.jl
+++ b/src/crop/almass.jl
@@ -239,7 +239,7 @@ function growcrop!(cs::CropState, model::SimulationModel)
points = curve.GDD[cs.phase]
for p in 1:length(points)
if points[p] == 99999
- cs.mature = true #FIXME is this correct?
+ !(cs.phase in (janfirst, sow)) && (cs.mature = true) #FIXME only in the last phase?
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]
diff --git a/src/farm/farm.jl b/src/farm/farm.jl
index a3aec89..7a3ad80 100644
--- a/src/farm/farm.jl
+++ b/src/farm/farm.jl
@@ -3,11 +3,6 @@
### This file is responsible for managing the farm module(s).
###
-##TODO what data do we need to gather from the farm submodel?
-## - area covered by each crop over time
-## - average height of each crop over time
-## - total income per year
-
#XXX Initially, we're only working with a single simple crop rotation.
# Later on, we need to figure out how to integrate several.
const CROPROTATION = ["winter rape", "winter wheat", "maize", "winter barley"]
@@ -79,6 +74,7 @@ function initfarms!(model::SimulationModel)
farmer.sowdates[field.id] = @rand(nextcrop.minsowdate:nextcrop.maxsowdate)
end
end
+ initfarmdata(model)
end
"""
diff --git a/src/farm/farmdata.jl b/src/farm/farmdata.jl
new file mode 100644
index 0000000..945d9ef
--- /dev/null
+++ b/src/farm/farmdata.jl
@@ -0,0 +1,41 @@
+### Persefone.jl - a model of agricultural landscapes and ecosystems in Europe.
+###
+### This file handles the data output of the farm submodel.
+###
+
+"""
+ initfarmdata()
+
+Create output files for each data group collected by the farm model.
+"""
+function initfarmdata(model::SimulationModel)
+ newdataoutput!(model, "fields", ["Date", "Crop", "Area", "Height"],
+ @param(farm.fieldoutfreq), savefielddata, croptrends)
+ #XXX add later: income per year
+end
+
+"""
+ savefielddata(model)
+
+Return a data table (to be printed to `fields.csv`), giving the current
+date, and the area and average of each crop in the landscape. May be called
+never, daily, monthly, yearly, or at the end of a simulation, depending on
+the parameter `farm.fieldoutfreq`.
+"""
+function savefielddata(model::SimulationModel)
+ croparea = Dict(c => 0.0 for c in keys(model.crops))
+ cropheights = Dict(c => [0.0, 0.0] for c in keys(model.crops))
+ for f in model.farmplots
+ c = cropname(f)
+ croparea[c] += ustrip(ha, @areaof(length(f.pixels))) # total area with this crop
+ cropheights[c][1] += 1 # number of fields with this crop
+ cropheights[c][2] += ustrip(cm, cropheight(f)) # summed height of the plants of this crop
+ end
+ data = []
+ for p in keys(croparea)
+ averageheight = cropheights[p][2] == 0 ? 0 :
+ round(cropheights[p][2]/cropheights[p][1], digits=2)
+ push!(data, [model.date, p, croparea[p], averageheight])
+ end
+ data
+end
--
GitLab