Rework crop models, introduce AbstractCropType and AbstractCropState

src/Persefone.jl

@@ -100,7 +100,10 @@ export
     populationtrends,
     visualiseoutput,
     savemodelobject,
-    loadmodelobject
+    loadmodelobject,
+    croptype,
+    cropname,
+    cropheight
 
 ## Import and define units and dimensions
 import Unitful: cm, m, km, ha, mg, g, kg, Length, Area, Mass
@@ -138,9 +141,11 @@ include("analysis/makieplots.jl")
 include("world/landscape.jl")
 include("world/weather.jl")
 
+include("crop/croptypes.jl")
 include("crop/farmplot.jl")
 include("crop/almass.jl")
 include("crop/simplecrop.jl")
+
 include("farm/farm.jl")
 
 include("nature/insects.jl")

src/core/simulation.jl

@@ -21,7 +21,7 @@ mutable struct AgricultureModel <: SimulationModel
     date::Date
     landscape::Matrix{Pixel}
     weather::Dict{Date,Weather}
-    crops::Dict{String,ALMaSS.CropType}
+    crops::Dict{String,AbstractCropType}
     farmers::Vector{Farmer}
     farmplots::Vector{AbstractFarmPlot}
     animals::Vector{Union{Animal,Nothing}}

src/crop/almass.jl

@@ -8,14 +8,20 @@
 module ALMaSS
 
 using Persefone:
-    AbstractFarmPlot,
+    AbstractCropType,
+    AbstractCropState,
     EventType,
     SimulationModel,
     fertiliser,
     initfields_fill_with!,
     maxtemp,
     mintemp
-import Persefone: stepagent!
+import Persefone:
+    Persefone,
+    stepagent!,
+    croptype,
+    cropname,
+    cropheight
 using Dates: Date, month, monthday
 using CSV: CSV
 
@@ -48,7 +54,7 @@ end
 The type struct for all crops. Currently follows the crop growth model as
 implemented in ALMaSS.
 """
-struct CropType
+struct CropType <: AbstractCropType
     #TODO make this into an abstract type and create subtypes for different
     # crop submodels (#70)
     name::String
@@ -62,6 +68,32 @@ struct CropType
     #issowable::Union{Function,Bool}
 end
 
+cropname(ct::CropType) = ct.name
+
+"""
+    CropState
+
+The state data for an ALMaSS vegetation point calculation as used in
+FarmPlot.
+"""
+mutable struct CropState <: AbstractCropState
+    #TODO add Unitful
+    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{EventType}
+end
+
+croptype(cs::CropState) = cs.croptype
+cropname(cs::CropState) = cropname(croptype(cs))
+cropheight(cs::CropState) = cs.height
+
+const FarmPlot = Persefone.FarmPlot{CropState}
+
 """
     Base.tryparse(type, str)
 
@@ -134,29 +166,6 @@ 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 <: AbstractFarmPlot
-    #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)
 
@@ -165,11 +174,11 @@ Initialise the model with its farm plots.
 function initfields!(model::SimulationModel)
     initfields_fill_with!(model) do model, x, y
         month(model.date) < 3 ? phase = ALMaSS.janfirst : phase = ALMaSS.marchfirst
-        FarmPlot(length(model.farmplots) + 1,
-                 [(x,y)],
-                 model.crops["natural grass"],
-                 phase,
-                 0.0, 0.0, 0.0, 0.0, Vector{EventType}())
+        FarmPlot(length(model.farmplots) + 1, [(x,y)],
+                 CropState(model.crops["natural grass"],
+                           phase,
+                           0.0, 0.0, 0.0, 0.0, Vector{EventType}())
+                 )
     end
 end
 
@@ -182,13 +191,14 @@ 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
+    cs = farmplot.crop_state
+    basetemp = cs.croptype.mingrowthtemp
     ismissing(basetemp) && (basetemp = 5.0)
     gdd = (maxtemp(model)+mintemp(model))/2 - basetemp
-    gdd > 0 && (farmplot.growingdegreedays += gdd)
+    gdd > 0 && (cs.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)
+    monthday(model.date) == (1,1) && (cs.phase = ALMaSS.janfirst)
+    monthday(model.date) == (3,1) && (cs.phase = ALMaSS.marchfirst)
     # update crop growth
     growcrop!(farmplot, model)
 end
@@ -202,9 +212,10 @@ end
 Sow the specified crop on this farmplot.
 """
 function sow!(cropname::String, farmplot::FarmPlot, model::SimulationModel)
+    cs = farmplot.crop_state
     createevent!(model, farmplot.pixels, sowing)
-    farmplot.croptype = model.crops[cropname]
-    farmplot.phase = ALMaSS.sow
+    cs.croptype = model.crops[cropname]
+    cs.phase = ALMaSS.sow
     #XXX test if the crop is sowable?
 end
 
@@ -214,10 +225,11 @@ end
 Harvest the crop on this farmplot.
 """
 function harvest!(farmplot::FarmPlot, model::SimulationModel)
+    cs = farmplot.crop_state
     createevent!(model, farmplot.pixels, harvesting)
-    farmplot.phase in [ALMaSS.harvest1, ALMaSS.harvest2] ?
-        farmplot.phase = ALMaSS.harvest2 :
-        farmplot.phase = ALMaSS.harvest1
+    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
@@ -233,25 +245,26 @@ 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]
+    cs = farmplot.crop_state
+    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
-            farmplot.height = curve.height[farmplot.phase][p]
-            farmplot.LAItotal = curve.LAItotal[farmplot.phase][p]
-            farmplot.LAIgreen = curve.LAIgreen[farmplot.phase][p]
+            cs.height = curve.height[cs.phase][p]
+            cs.LAItotal = curve.LAItotal[cs.phase][p]
+            cs.LAIgreen = curve.LAIgreen[cs.phase][p]
             return
         else
-            gdd = farmplot.growingdegreedays
+            gdd = cs.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]
+                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
@@ -263,53 +276,4 @@ function growcrop!(farmplot::FarmPlot, model::SimulationModel)
     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

src/crop/croptypes.jl

+### Persefone.jl - a model of agricultural landscapes and ecosystems in Europe.
+###
+### This file contains abstract types that crop models must derive from.
+###
+
+abstract type AbstractCropType end
+cropname(t::AbstractCropType) = t.name
+
+abstract type AbstractCropState end
+croptype(s::AbstractCropState) = s.croptype
+cropname(s::AbstractCropState) = cropname(cropytype(s))
+# cropheight(s::AbstractCropState) = s.cropheight
+# cropcover(s::AbstractCropState) = s.cropcover

src/crop/farmplot.jl

@@ -214,6 +214,16 @@ end
 
 
 
+mutable struct FarmPlot{T <: AbstractCropState} <: AbstractFarmPlot
+    const id::Int64
+    pixels::Vector{Tuple{Int64, Int64}}
+    crop_state :: T
+end
+
+croptype(f::FarmPlot{T}) where {T} = croptype(f.crop_state)
+cropname(f::FarmPlot{T}) where {T} = cropname(croptype(f))
+cropheight(f::FarmPlot{T}) where {T} = cropheight(f.crop_state)
+
 """
     initfields_fill_with!(make_farmplot_fn, model)
 
@@ -246,3 +256,52 @@ function initfields_fill_with!(make_farmplot_fn::Function, model::SimulationMode
     end
     @info "Initialised $n farm plots."
 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 :
+              croptype(model.farmplots[model.landscape[pos...].fieldid])
+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) ? "" : cropname(model.farmplots[field])
+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 :
+              cropheight(model.farmplots[model.landscape[pos...].fieldid])
+end

src/crop/simplecrop.jl

 module SimpleCrop
 
-using Persefone: AbstractFarmPlot, SimulationModel, initfields_fill_with!
-import Persefone: stepagent!
+using Persefone:
+    AbstractCropType,
+    AbstractCropState,
+    SimulationModel,
+    initfields_fill_with!
+import Persefone:
+    Persefone,
+    stepagent!,
+    croptype,
+    cropname,
+    cropheight
 
-mutable struct FarmPlot <: AbstractFarmPlot
-    const id::Int64
-    pixels::Vector{Tuple{Int64, Int64}}
+struct CropType <: AbstractCropType
+    name::String
 end
 
+cropname(ct::CropType) = ct.name
+
+mutable struct CropState <: AbstractCropState
+    croptype::CropType
+    cropheight::Float64
+end
+
+croptype(cs::CropState) = cs.croptype
+cropname(cs::CropState) = cropname(croptype(cs))
+cropheight(cs::CropState) = cs.height
+
+const FarmPlot = Persefone.FarmPlot{CropState}
+
 """
     stepagent!(farmplot, model)
 
@@ -24,7 +45,7 @@ Initialise the model with its farm plots.
 """
 function initfields!(model::SimulationModel)
     initfields_fill_with!(model) do model, x, y
-        FarmPlot(length(model.farmplots) + 1, [(x,y)])
+        FarmPlot(length(model.farmplots) + 1, [(x,y)], CropState(CropType("nothing"), 0.0))
     end
 end
 

test/landscape_tests.jl

@@ -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.ALMaSS.averagefieldsize(model) == 5.37
+    @test Ps.averagefieldsize(model) == 5.37
     @test count(f -> ismissing(f.fieldid), model.landscape) == 1685573
     @test length(Ps.farmplot((800,800), model).pixels) == 4049               
 end

test/nature_tests.jl

@@ -47,10 +47,15 @@ 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.ALMaSS.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.0, 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))