Added data recording and visualisation for the farm model

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")

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

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]

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
 
 """

src/farm/farmdata.jl

+### 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