Implemented ALMaSS crop growth function

This is a major commit, probably rather too large, but anyway. There
are no tests yet, so no guarantee for correctness. But it's good
progress :-)

data/crop_data_general.csv

 name,minsowdate,maxsowdate,minharvestdate,maxharvestdate,mingrowthtemp
-"winter barley","15 September","30 September",NA,NA,NA
-"spring barley","1 March","10 April",NA,NA,NA
-"peas/beans",NA,NA,NA,NA,NA
+"winter barley","15 September","30 September",NA,NA,0
+"spring barley","1 March","10 April",NA,NA,0
+"peas/beans",NA,NA,NA,NA,5
 "spring rape",NA,NA,NA,NA,NA
 "winter rape",NA,NA,NA,NA,NA
 "winter rye",NA,NA,NA,NA,NA
-"winter wheat","15 October","31 October",NA,NA,NA
+"winter wheat","15 October","31 October",NA,NA,0
 "beet",NA,NA,NA,NA,NA
-"maize",NA,NA,NA,NA,NA
+"maize",NA,NA,NA,NA,8
 "permanent grassland (grazed)",NA,NA,NA,NA,NA
 "permanent grassland (seeded)",NA,NA,NA,NA,NA
 "fodder/clover",NA,NA,NA,NA,NA
 "natural grass",NA,NA,NA,NA,NA
-"potatoes",NA,NA,NA,NA,NA
-"undersown spring barley",NA,NA,NA,NA,NA
+"potatoes",NA,NA,NA,NA,4
+"undersown spring barley",NA,NA,NA,NA,0
 "carrots",NA,NA,NA,NA,NA
 "oats",NA,NA,NA,NA,NA
 "permanent set-aside",NA,NA,NA,NA,NA

src/core/simulation.jl

@@ -69,12 +69,15 @@ function initmodel(settings::Dict{String, Any})
         weather = initweather(settings["world.weatherfile"],
                               settings["core.startdate"],
                               settings["core.enddate"])
+        crops = readcropparameters(settings["crop.cropfile"],
+                                   settings["crop.growthfile"])
         space = GridSpace(size(landscape), periodic=false)
         properties = Dict{Symbol,Any}(:settings=>settings,
                                       :logger=>logger,
                                       :date=>settings["core.startdate"],
                                       :landscape=>landscape,
                                       :weather=>weather,
+                                      :crops=>crops,
                                       :dataoutputs=>dataoutputs,
                                       :events=>events)
         model = AgentBasedModel(Union{Farmer,Animal,FarmPlot}, space, properties=properties,

src/crop/crops.jl

@@ -4,7 +4,6 @@
 ###
 
 #TODO write tests for input functions
-#TODO write actual growth function
 
 """
     GrowthPhase
@@ -96,14 +95,14 @@ end
 Parse a CSV file containing the required parameter values for each crop
 (as produced from the original ALMaSS file by `convert_almass_data.py`).
 """
-function readcropparameters(generalcropfile::String, cropgrowthfile::String)
+function readcropparameters(generalcropfile::String, growthfile::String)
     @debug "Reading crop parameters"
     cropdata = CSV.File(generalcropfile, missingstring="NA", dateformat="d U",
                         types=[String,Date,Date,Date,Date,Float64])
-    growthdata = CSV.File(cropgrowthfile, missingstring="NA",
+    growthdata = CSV.File(growthfile, missingstring="NA",
                           types=[Int,String,String,GrowthPhase,String,
                                  Float64,Float64,Float64,Float64])
-    croptypes = Vector{CropType}()
+    croptypes = Dict{String,CropType}()
     for crop in cropdata
         cropgrowthdata = growthdata |> filter(x -> !ismissing(x.crop_name) &&
                                               x.crop_name == crop.name)
@@ -111,9 +110,9 @@ function readcropparameters(generalcropfile::String, cropgrowthfile::String)
                                     filter(x -> x.nutrient_status=="high"))
         lownuts = buildgrowthcurve(cropgrowthdata |>
                                    filter(x -> x.nutrient_status=="low"))
-        append!(croptypes, [CropType(crop.name, crop.minsowdate, crop.maxsowdate,
-                                     crop.minharvestdate, crop.maxharvestdate,
-                                     crop.mingrowthtemp, highnuts, lownuts)])
+        croptypes[crop.name] = CropType(crop.name, crop.minsowdate, crop.maxsowdate,
+                                        crop.minharvestdate, crop.maxharvestdate,
+                                        crop.mingrowthtemp, highnuts, lownuts)
     end
     croptypes
 end

src/crop/farmplot.jl

@@ -14,25 +14,13 @@ This is the spatial unit with which the crop growth model and the farm model wor
 @agent FarmPlot GridAgent{2} begin
     pixels::Vector{Tuple{Int64, Int64}}
     croptype::CropType
+    phase::GrowthPhase
     growingdegreedays::Float64
     height::Float64
-    biomass::Float64
+    LAItotal::Float64
+    LAIgreen::Float64
+    #biomass::Float64 #XXX I need to figure out how to calculate this
     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
 
 """
@@ -56,7 +44,11 @@ function initfields!(model::AgentBasedModel)
                 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)
+                #XXX does this phase calculation work out?
+                month(model.date) < 3 ? phase = janfirst : phase = marchfirst
+                fp = add_agent!((x,y), FarmPlot, model, [(x,y)],
+                                model.crops["natural grass"], phase, false,
+                                0.0, 0.0, 0.0, 0.0, Vector{EventType}())
                 model.landscape[x,y].fieldid = fp.id
                 convertid[rawid] = fp.id
                 n += 1
@@ -66,6 +58,99 @@ function initfields!(model::AgentBasedModel)
     @info "Initialised $n farm plots."
 end
 
+"""
+    stepagent!(farmplot, model)
+
+Update a farm plot by one day.
+"""
+function stepagent!(farmplot::FarmPlot, model::AgentBasedModel)
+    # 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
+    ismissing(basetemp) && (basetemp = 5.0)        
+    gdd = (maxtemp(model)+mintemp(model))/2 - basetemp
+    gdd > 0 && (farmplot.growingdegreedays += gdd)
+    # update the phase on key dates
+    monthday(model.date) == (1,1) && (farmplot.phase = janfirst)
+    monthday(model.date) == (3,1) && (farmplot.phase = marchfirst)
+    # update crop growth
+    growcrop!(farmplot, model)
+end
+
+
+## CROP MANAGEMENT AND GROWTH FUNCTIONS
+
+"""
+    sow!(cropname, farmplot, model)
+
+Sow the specified crop on this farmplot.
+"""
+function sow!(cropname::String, farmplot::FarmPlot, model::AgentBasedModel)
+    createevent!(model, farmplot.pixels, sowing)
+    farmplot.croptype = model.crops[cropname]
+    farmplot.phase = sow
+    #XXX test if the crop is sowable?
+end
+
+"""
+    harvest!(farmplot, model)
+
+Harvest the crop on this farmplot.
+"""
+function harvest!(farmplot::FarmPlot, model::AgentBasedModel)
+    createevent!(model, farmplot.pixels, harvesting)
+    farmplot.phase in [harvest1, harvest2] ?
+        farmplot.phase = harvest2 :
+        farmplot.phase = harvest1
+    # height & LAI will be automatically adjusted by the growth function
+    #TODO calculate and return yield
+end
+
+#TODO fertilise!()
+#TODO spray!()
+#TODO till!()
+
+"""
+    growcrop!(farmplot, model)
+
+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::AgentBasedModel)
+    fertiliser in farmplot.events ?
+        curve = farmplot.croptype.lownutrientgrowth :
+        curve = farmplot.croptype.highnutrientgrowth
+    points = curve.GDD[farmplot.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]
+            return
+        else
+            gdd = farmplot.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]
+                return
+            end
+            #XXX To be precise, we ought to figure out if one or more inflection
+            # points have been passed between yesterday and today, and calculate the
+            # growth exactly up to the inflection point with the old slope, and onward
+            # with the new slope. Not doing so will introduce a small amount of error,
+            # although I think this is acceptable.
+        end
+    end
+end
+                 
+
+## UTILITY FUNCTIONS
+
 """
     averagefieldsize(model)
 
@@ -101,14 +186,3 @@ 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/parameters.toml

@@ -35,5 +35,5 @@ insectmodel = ["season", "habitat", "pesticides", "weather"] # factors affecting
 [crop]
 cropmodel = "almass" # crop growth model to use, "almass" or "aquacrop"
 cropfile = "data/crop_data_general.csv" # file with general crop parameters
-cropgrowthfile = "data/almass_crop_growth_curves.csv" # file with crop growth parameters	
+growthfile = "data/almass_crop_growth_curves.csv" # file with crop growth parameters	
 

src/world/landscape.jl

@@ -8,7 +8,7 @@
 ## Do not change the order of this enum, or initlandscape() will break!
 
 "The types of landscape event that can be simulated"
-@enum EventType tillage sowing fertiliser pesticide harvest
+@enum EventType tillage sowing fertiliser pesticide harvesting
 
 """
     Pixel

test/test_parameters.toml

@@ -31,5 +31,7 @@ indoutfreq = "end" # output frequency individual-level data, daily/monthly/yearl
 insectmodel = ["season", "habitat", "pesticides"] # which factors affect insect growth ("weather" is not yet implemented)
 	
 [crop]
-cropmodel = "linear" # crop growth model to use, "linear" or "aquacrop" (not yet implemented)
+cropmodel = "almass" # crop growth model to use, "almass" or "aquacrop"
+cropfile = "data/crop_data_general.csv" # file with general crop parameters
+growthfile = "data/almass_crop_growth_curves.csv" # file with crop growth parameters