Started working on crop growth model

data/convert_almass_data.py

 #!/usr/bin/python3
 ###
-### ALMaSS uses a pretty crazy file format to specify its crop growth parameters.
-### This script converts the data to a saner long-table CSV format.
+### ALMaSS uses a pretty convoluted file format to specify its crop growth parameters.
+### This script converts the data to a more easily understandable long-table CSV format.
 ###
 ### Daniel Vedder, 02/08/2023
 ###
 
-import time
-
 input_file = "almass_crop_growth_curves.pre"
 output_file = "almass_crop_growth_curves.csv"
 
@@ -82,13 +80,13 @@ def interpret_crop_number(i):
     elif i in [4, 104]: crop_type = "winter wheat"
     elif i in [5, 105]: crop_type = "winter rye"
     elif i in [6, 106]: crop_type = "oats"
-    #elif i in [7, 107]: crop_type = "triticale" #FIXME 107 not defined, but probable
-    elif i == 7: crop_type = "triticale"
+    elif i in [7, 107]: crop_type = "triticale" #XXX 107 not defined, but probable
+    #elif i == 7: crop_type = "triticale"
     elif i in [8, 108]: crop_type = "maize"
     elif i in [13, 113]: crop_type = "undersown spring barley"
     elif i in [21, 121]: crop_type = "spring rape"
-    elif i == 22: crop_type = "winter rape"
-    #elif i in [22, 122]: crop_type = "winter rape" #FIXME 122 not defined, but probable
+    #elif i == 22: crop_type = "winter rape"
+    elif i in [22, 122]: crop_type = "winter rape" #XXX 122 not defined, but probable
     elif i == 25: crop_type = "permanent grassland, low yield"
     elif i == 26: crop_type = "permanent grassland, grazed"
     elif i == 27: crop_type = "permanent grassland, seeded"
@@ -97,8 +95,8 @@ def interpret_crop_number(i):
     elif i == 30: crop_type = "peas/beans"
     elif i in [41, 141]: crop_type = "carrots"
     elif i in [50, 150]: crop_type = "potatoes"
-    #elif i in [60, 160]: crop_type = "beet" #FIXME 160 not defined, but probable
-    elif i == 60: crop_type = "beet"
+    elif i in [60, 160]: crop_type = "beet" #XXX 160 not defined, but probable
+    #elif i == 60: crop_type = "beet"
     #elif i in [65, 165]: crop_type = "NA" #FIXME unknown numbers
     elif i == 70: crop_type = "lucerne"
     elif i == 80: crop_type = "tulips"
@@ -108,7 +106,6 @@ def interpret_crop_number(i):
     elif i == 94: crop_type = "lawn"
     #elif i == 99: crop_type = "NA" #FIXME unknown number
     elif i == 112: crop_type = "heath"
-    #elif i == 121: crop_type = "NA" #FIXME unknown number
     else:
         crop_type = "NA"
         print("WARNING: unknown crop type "+str(i))

src/Persefone.jl

@@ -80,8 +80,9 @@ include("core/output.jl")
 include("world/landscape.jl")
 include("world/weather.jl")
 
-include("farm/farm.jl")
 include("crop/crops.jl")
+include("crop/farmplot.jl")
+include("farm/farm.jl")
 
 include("nature/insects.jl")
 include("nature/energy.jl")

src/core/simulation.jl

@@ -80,8 +80,8 @@ function initmodel(settings::Dict{String, Any})
         model = AgentBasedModel(Union{Farmer,Animal,FarmPlot}, space, properties=properties,
                                 rng=StableRNG(settings["core.seed"]), warn=false)
         saveinputfiles(model)
-        initfarms!(model)
         initfields!(model)
+        initfarms!(model)
         initnature!(model)
         model
     end

src/crop/crops.jl

@@ -3,97 +3,85 @@
 ### This file is responsible for managing the crop growth modules.
 ###
 
-"The crop types simulated by the model"
-@enum CropType fallow wheat barley maize rapeseed
+@enum GrowthPhase janfirst sow marchfirst harvest1 harvest2 vegphase
 
-#XXX not sure whether it makes sense to have this as an agent type,
-# or perhaps better a grid property?
 """
-    FarmPlot
+    CropType
 
-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.
+The type struct for all crops. Currently follows the crop
+growth model as implemented in ALMaSS.
 """
-@agent FarmPlot GridAgent{2} begin
-    pixels::Vector{Tuple{Int64, Int64}}
-    croptype::CropType
-    height::Float64
-    biomass::Float64
-    #TODO
+@kwdef struct CropType
+    name::String
+    sowdate::Date
+    harvestdate::Date
+    mingrowthtemp::Float64
+    inflectionpoints::Dict{GrowthPhase,Vector{Float64}}
+    slopes::Dict{GrowthPhase,Vector{Float64}}
+    #issowable::Union{Function,Bool}=true
 end
 
-"""
-    stepagent!(farmplot, model)
+let croptypes = Dict{String,CropType}()
+    global function registercrop(ct::CropType)
+        croptypes[ct.name] = ct
+    end
 
-Update a farm plot by one day.
-"""
-function stepagent!(farmplot::FarmPlot, model::AgentBasedModel)
-    #TODO
+    global function getcrop(name::String)
+        croptypes[name]
+    end
 end
 
 """
-    initfields!(model)
+    @dday(datestring)
 
-Initialise the model with its farm plots.
+A utility macro to make it quicker to create dates in the
+format "21 August".
 """
-function initfields!(model::AgentBasedModel)
-    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[objectid].pixels, (x,y))
-            else
-                fp = add_agent!((x,y), FarmPlot, model, [(x,y)], fallow, 0.0, 0.0)
-                model.landscape[x,y].fieldid = fp.id
-                convertid[rawid] = fp.id
-                n += 1
-            end
-        end
-    end
-    @info "Initialised $n farm plots."
+macro dday(date::String)
+    :(Date($(date), dateformat"d U"))
 end
 
-"""
-    averagefieldsize(model)
+let croplist = Dict{String,CropType}()
+    
+    global function newcrop(;kwargs...)
+        crop = CropType(;kwargs...)
+        croplist[crop.name] = crop
+    end
 
-Calculate the average field size in hectares for the model landscape.
-"""
-function averagefieldsize(model::AgentBasedModel)
-    conversionfactor = 100 #our pixels are currently 10x10m, so 100 pixels per hectare
-    sizes::Vector{Float64} = []
-    for a in allagents(model)
-        (typeof(a) != FarmPlot) && continue
-        push!(sizes, size(a.pixels)[1]/conversionfactor)
+    global function cropdefinition(name)
+        croplist[name]
     end
-    round(sum(sizes)/size(sizes)[1], digits=2)
 end
 
-"""
-    croptype(model, position)
+#
+# XXX It would be nicer to be able to define crops with a macro
+# (like @species), but since macros are so tricky to get right,
+# I'll leave that for later.
+#
+# """
+#     @crop
+#
+# A macro to create and register a new crop type. Pass a list
+# of keyword arguments to instantiate each field of the
+# [`CropType`](@ref) struct.
+# """
+# macro crop(;kwargs...) # macros don't take kwargs
+#     quote
+#         registercrop(CropType($(esc(kwargs))))
+#     end
+# end
+#
+# @crop "winter wheat" begin
+#     minsowdate = @dday("1 November")
+#     maxsowdate = @dday("31 November")
+#     growthparameter = 1.0
+#     issowable = true
+# end
 
-Return the crop at this position, or nothing if there is no crop here (utility wrapper).
-"""
-function croptype(pos::Tuple{Int64,Int64}, model::AgentBasedModel)
-    ismissing(model.landscape[pos...].fieldid) ? nothing :
-              model[model.landscape[pos...].fieldid].croptype
-end
+newcrop(name="winter wheat",
+        sowdate = @dday("1 November"),
+        harvestdate = @dday("1 July"),
+        mingrowthtemp = 5.0,
+        growthparameter = 1.0,
+        issowable = true)
 
-"""
-    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::AgentBasedModel)
-    ismissing(model.landscape[pos...].fieldid) ? nothing :
-              model[model.landscape[pos...].fieldid].height
-end

src/crop/farmplot.jl

+### Persefone - a socio-economic-ecological model of European agricultural landscapes.
+###
+### This file contains code for the fields that farmers manage.
+###
+
+#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.
+"""
+@agent FarmPlot GridAgent{2} begin
+    pixels::Vector{Tuple{Int64, Int64}}
+    croptype::CropType
+    growingdegreedays::Float64
+    height::Float64
+    biomass::Float64
+    events::Vector{EventType}
+    #TODO
+end
+
+"""
+    stepagent!(farmplot, model)
+
+Update a farm plot by one day.
+"""
+function stepagent!(farmplot::FarmPlot, model::AgentBasedModel)
+    # update growing degree days
+    gdd = meantemp(model) - farmplot.croptype.mingrowthtemp
+    gdd > 0 && (farmplot.growingdegreedays += gdd)
+    # update crop growth
+    growcrop!(farmplot, model)
+    #TODO expand?
+end
+
+"""
+    initfields!(model)
+
+Initialise the model with its farm plots.
+"""
+function initfields!(model::AgentBasedModel)
+    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[objectid].pixels, (x,y))
+            else
+                fp = add_agent!((x,y), FarmPlot, model, [(x,y)], fallow, 0.0, 0.0)
+                model.landscape[x,y].fieldid = fp.id
+                convertid[rawid] = fp.id
+                n += 1
+            end
+        end
+    end
+    @info "Initialised $n farm plots."
+end
+
+"""
+    averagefieldsize(model)
+
+Calculate the average field size in hectares for the model landscape.
+"""
+function averagefieldsize(model::AgentBasedModel)
+    conversionfactor = 100 #our pixels are currently 10x10m, so 100 pixels per hectare
+    sizes::Vector{Float64} = []
+    for a in allagents(model)
+        (typeof(a) != FarmPlot) && continue
+        push!(sizes, size(a.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::AgentBasedModel)
+    ismissing(model.landscape[pos...].fieldid) ? nothing :
+              model[model.landscape[pos...].fieldid].croptype
+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::AgentBasedModel)
+    ismissing(model.landscape[pos...].fieldid) ? nothing :
+              model[model.landscape[pos...].fieldid].height
+end
+
+
+"""
+    growcrop!(farmplot, model)
+
+Apply the relevant crop growth model to update the plant height on this farm plot.
+"""
+function growcrop!(farmplot::FarmPlot, model::AgentBasedModel)
+    crop = farmplot.croptype
+    #TODO
+end

src/farm/farm.jl

@@ -9,7 +9,9 @@
 This is the agent type for the farm ABM. (Not yet implemented.)
 """
 @agent Farmer GridAgent{2} begin
-    #TODO
+    fields::Vector{FarmPlot}
+    croprotation::Vector{CropType}
+    #TODO add AES
 end
 
 """