### Persefone.jl - a model of agricultural landscapes and ecosystems in Europe. ### ### This file is responsible for managing the crop growth modules. ### #TODO write tests for input functions module ALMaSS using Persefone: AbstractFarmPlot, EventType, SimulationModel, fertiliser, initfields_fill_with!, maxtemp, mintemp import Persefone: stepagent! using Dates: Date, month, monthday using CSV: CSV """ GrowthPhase ALMaSS crop growth curves are split into five phases, triggered by seasonal dates or agricultural events. """ @enum GrowthPhase janfirst sow marchfirst harvest1 harvest2 """ CropCurveParams The values in this struct define one crop growth curve. """ struct CropCurveParams #TODO add Unitful curveID::Int highnutrients::Bool GDD::Dict{GrowthPhase,Vector{Float64}} LAItotal::Dict{GrowthPhase,Vector{Float64}} LAIgreen::Dict{GrowthPhase,Vector{Float64}} height::Dict{GrowthPhase,Vector{Float64}} end """ CropType The type struct for all crops. Currently follows the crop growth model as implemented in ALMaSS. """ struct CropType #TODO make this into an abstract type and create subtypes for different # crop submodels (#70) name::String minsowdate::Union{Missing,Date} maxsowdate::Union{Missing,Date} minharvestdate::Union{Missing,Date} maxharvestdate::Union{Missing,Date} mingrowthtemp::Union{Missing,Float64} highnutrientgrowth::Union{Missing,CropCurveParams} lownutrientgrowth::Union{Missing,CropCurveParams} #issowable::Union{Function,Bool} end """ Base.tryparse(type, str) Extend `tryparse` to allow parsing GrowthPhase values. (Needed to read in the CSV parameter file.) """ function Base.tryparse(type::Type{GrowthPhase}, str::String) str == "janfirst" ? janfirst : str == "sow" ? sow : str == "marchfirst" ? marchfirst : str == "harvest1" ? harvest1 : str == "harvest2" ? harvest2 : nothing end """ buildgrowthcurve(data) Convert a list of rows from the crop growth data into a CropCurveParams object. """ function buildgrowthcurve(data::Vector{CSV.Row}) isempty(data) && return missing GDD = Dict(janfirst=>Vector{Float64}(), sow=>Vector{Float64}(), marchfirst=>Vector{Float64}(), harvest1=>Vector{Float64}(), harvest2=>Vector{Float64}()) LAItotal = Dict(janfirst=>Vector{Float64}(), sow=>Vector{Float64}(), marchfirst=>Vector{Float64}(), harvest1=>Vector{Float64}(), harvest2=>Vector{Float64}()) LAIgreen = Dict(janfirst=>Vector{Float64}(), sow=>Vector{Float64}(), marchfirst=>Vector{Float64}(), harvest1=>Vector{Float64}(), harvest2=>Vector{Float64}()) height = Dict(janfirst=>Vector{Float64}(), sow=>Vector{Float64}(), marchfirst=>Vector{Float64}(), harvest1=>Vector{Float64}(), harvest2=>Vector{Float64}()) for e in data append!(GDD[e.growth_phase], e.GDD) append!(LAItotal[e.growth_phase], e.LAI_total) append!(LAIgreen[e.growth_phase], e.LAI_green) append!(height[e.growth_phase], e.height) end CropCurveParams(data[1].curve_id, data[1].nutrient_status=="high", GDD, LAItotal, LAIgreen, height) end """ readcropparameters(generalcropfile, cropgrowthfile) 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, 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(growthfile, missingstring="NA", types=[Int,String,String,GrowthPhase,String, Float64,Float64,Float64,Float64]) croptypes = Dict{String,CropType}() for crop in cropdata cropgrowthdata = growthdata |> filter(x -> !ismissing(x.crop_name) && x.crop_name == crop.name) highnuts = buildgrowthcurve(cropgrowthdata |> filter(x -> x.nutrient_status=="high")) lownuts = buildgrowthcurve(cropgrowthdata |> filter(x -> x.nutrient_status=="low")) croptypes[crop.name] = CropType(crop.name, crop.minsowdate, crop.maxsowdate, crop.minharvestdate, crop.maxharvestdate, crop.mingrowthtemp, highnuts, lownuts) end 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) 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}()) end end """ stepagent!(farmplot, model) Update a farm plot by one day. """ 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 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 = ALMaSS.janfirst) monthday(model.date) == (3,1) && (farmplot.phase = ALMaSS.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::SimulationModel) createevent!(model, farmplot.pixels, sowing) farmplot.croptype = model.crops[cropname] farmplot.phase = ALMaSS.sow #XXX test if the crop is sowable? end """ harvest!(farmplot, model) Harvest the crop on this farmplot. """ function harvest!(farmplot::FarmPlot, model::SimulationModel) createevent!(model, farmplot.pixels, harvesting) farmplot.phase in [ALMaSS.harvest1, ALMaSS.harvest2] ? farmplot.phase = ALMaSS.harvest2 : farmplot.phase = ALMaSS.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::SimulationModel) 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) 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