almass.jl

### Persefone.jl - a model of agricultural landscapes and ecosystems in Europe.
###
### This file is responsible for managing the crop growth modules.
###

#FIXME Looking at fields.pdf in the Persefone output, crop maturity is not recognised
# properly, so most crops (except maize) are never planted. At the same time, maize never
# seems to start growing.

#TODO write tests to compare our output to the output of the original ALMaSS algorithm

module ALMaSS

using Persefone:
    AnnualDate,
    Management,
    Length,
    cm,
    SimulationModel,
    fertiliser,
    maxtemp,
    mintemp

import Persefone:
    stepagent!,
    croptype,
    cropname,
    cropheight,
    cropcover,
    cropyield,
    sow!,
    harvest!,
    isharvestable,
    bounds

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{Length{Float64}}}
end

function Base.show(io::IO, ::MIME"text/plain", curve::CropCurveParams)
    println(io, "CropCurveParams")
    println(io, "    curveID       = $(curve.curveID)")
    println(io, "    highnutrients = $(curve.highnutrients)")
    for (dict, name) in [
        (curve.GDD, "GDD"),
        (curve.LAItotal, "LAItotal"),
        (curve.LAIgreen, "LAIgreen"),
        ]
        println("    $name:")
        for (phase, arr) in dict
            println("        $phase = $arr")
        end
    end
    # Manually print height array for each phase so that we avoid the
    # lengthy type information before printing the array.
    println("    height:")
    for (phase, arr) in curve.height
        arrstr = "[" * join(split(string(arr), "[")[2:end], "[")
        println("        $phase = $arrstr")
    end
end

"""
    CropType

The type struct for all crops. Currently follows the crop growth model as
implemented in ALMaSS.
"""
struct CropType
    #FIXME this needs thinking about. The sowing and harvest dates belong in the farm model,
    # not here. Also, we need to harmonise crops across the crop growth models.
    name::String
    group::String
    minsowdate::Union{Missing,AnnualDate}
    maxsowdate::Union{Missing,AnnualDate}
    minharvestdate::Union{Missing,AnnualDate}
    maxharvestdate::Union{Missing,AnnualDate}
    mingrowthtemp::Union{Missing,Float64}
    highnutrientgrowth::Union{Missing,CropCurveParams}
    lownutrientgrowth::Union{Missing,CropCurveParams}
end

cropname(ct::CropType) = ct.name

"""
    CropState

The state data for an ALMaSS vegetation point calculation.  Usually
part of a `FarmPlot`.
"""
mutable struct CropState
    #TODO add Unitful
    croptype::CropType
    phase::GrowthPhase
    growingdegreedays::Float64
    height::Length{Float64}
    LAItotal::Float64
    LAIgreen::Float64
    #biomass::Float64 #XXX I need to figure out how to calculate this
    mature::Bool #TODO how do we determine this?
    events::Vector{Management} #FIXME does this do anything?
end

croptype(cs::CropState) = cs.croptype
cropname(cs::CropState) = cropname(croptype(cs))
cropheight(cs::CropState) = cs.height
cropcover(cs::CropState) = 0.0  # TODO: related to LAItotal, LAIgreen?
cropyield(cs::CropState) = 0.0  # TODO: units? needs biomass?
isharvestable(cs::CropState) = cs.mature

"Temperature to solar conversion factor for C3 plants."
const temperature_to_solar_conversion_c3::Vector{Float64} = [
    0,   0,    0,    0,    0,    0,    0,    0,    0,    0,     # -30°C to -21°C
    0,   0,    0,    0,    0,    0,    0,    0,    0,    0,     # -20°C to -11°C
    0,   0,    0,    0,    0,    0,    0,    0,    0,    0,     # -10°C to  -1°C
    0,   0,    0,    0,    0,    0.28, 0.56, 0.84, 1.12, 1.4,   #   0°C to   9°C
    1.4, 1.4,  1.4,  1.4,  1.4,  1.4,  1.4,  1.4,  1.4,  1.4,   #  10°C to  19°C
    1.4, 1.4,  1.4,  1.4,  1.4,  1.4,  1.26, 1.12, 0.98, 0.84,  #  20°C to  29°C
    0.7, 0.56, 0.42, 0.28, 0.14, 0,    0,    0,    0,    0,     #  30°C to  39°C
    0,   0,    0,    0,    0,    0,    0,    0,    0,    0,     #  40°C to  49°C
    0                                                           #  50°C
]

"Temperature to solar conversion factor for C4 plants."
const temperature_to_solar_conversion_c4::Vector{Float64} = [
    0,        0,        0,        0,        0,    0,   0,    0,    0,        0,         # -30°C to -21°C
    0,        0,        0,        0,        0,    0,   0,    0,    0,        0,         # -20°C to -11°C
    0,        0,        0,        0,        0,    0,   0,    0,    0,        0,         #   0°C to   9°C
    0,        0,        0,        0,        0,    0,   0,    0,    0.242857, 0.485714,  #  10°C to  19°C
    0.728571, 0.971429, 1.214286, 1.457143, 1.7,  1.7, 1.7,  1.7,  1.7,      1.7,       #  20°C to  29°C
    1.7,      1.7,      1.7,      1.7,      1.7,  1.7, 1.7,  1.7,  1.7,      1.7,       #  30°C to  39°C
    1.7,      1.7,      1.7,      1.7,      1.7,  1.7, 1.53, 1.36, 1.19,     1.02,      #  40°C to  49°C
    0.85,     0.68,     0.51,     0.34,     0.17, 0,   0,    0,    0,        0,         #  50°C
    0
]

"""
    solar_conversion_c3(temperature)

Solar conversion factor (no units) for C3 plants.
"""
function solar_conversion_c3(temperature)
    if temperature < -30 || temperature > 50
        return 0.0
    end
    idx = Int(floor(0.5 + temperature)) + 30 + 1
    return temperature_to_solar_conversion_c3[idx]
end

"""
    solar_conversion_c3(temperature)

Solar conversion factor (no units) for C4 plants.
"""
function solar_conversion_c4(temperature)
    if temperature < -30 || temperature > 50
        return 0.0
    end
    idx = Int(floor(0.5 + temperature)) + 30 + 1
    return temperature_to_solar_conversion_c4[idx]
end

"""
    Base.tryparse(type, str)

Extend `tryparse` to allow parsing GrowthPhase values.
(Needed to read in the CSV parameter file.)
"""
function Base.tryparse(::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{Length{Float64}}(), sow=>Vector{Length{Float64}}(),
                  marchfirst=>Vector{Length{Float64}}(), harvest1=>Vector{Length{Float64}}(),
                  harvest2=>Vector{Length{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 * cm)  # assumes `height` is in units of `cm`
    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",
                        types=[String,AnnualDate,AnnualDate,AnnualDate,AnnualDate,Float64,String])
    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"))
        # TODO: set crop group temporarily until there is a column in
        # the csv file
        crop_group = "CROP_GROUP_NOT_SET"
        croptypes[crop.name] = CropType(crop.name, crop_group, crop.minsowdate, crop.maxsowdate,
                                        crop.minharvestdate, crop.maxharvestdate,
                                        crop.mingrowthtemp, highnuts, lownuts)
    end
    croptypes
end

function setphase!(cs::CropState, phase::GrowthPhase)
    cs.phase = phase
    # TODO: should cs.growingdegreedays be set here?
    # cs.growingdegreedays = 0.0

    fertiliser in cs.events ?
        curve = cs.croptype.lownutrientgrowth :
        curve = cs.croptype.highnutrientgrowth
    if ismissing(curve)
        error("Missing curve = $curve")
    end
    GDD = curve.GDD[cs.phase]
    height = curve.height[cs.phase]
    LAItotal = curve.LAItotal[cs.phase]
    LAIgreen = curve.LAIgreen[cs.phase]

    if (length(GDD) == 0
        || length(height) == 0
        || length(LAItotal) == 0
        || length(LAIgreen) == 0)
        error("One of these has length zero:\n  GDD = $GDD\n  height = $height\n  LAItotal = $LAItotal\n  LAIgreen = $LAIgreen")
    end
    if GDD[1] == -1.0
        # set height, LAItotal, LAIgreen
        # TODO: what about when
        #           GDD[1] == curve.GDD[cs.phase][1] == 99999.0
        #       ???
        # TODO: cropname(cs) == "no growth" has curve.GDD[harvest1] == 99999.0
        cs.height = height[1]
        cs.LAItotal = LAItotal[1]
        cs.LAIgreen = LAIgreen[1]
    end
end

"""
    stepagent!(cropstate, model)

Update a farm plot by one day.
"""
function stepagent!(cs::CropState, model::SimulationModel)
    # update the phase on key dates
    if monthday(model.date) == (1, 1)
        setphase!(cs, ALMaSS.janfirst)
        cs.growingdegreedays = 0.0
    elseif monthday(model.date) == (3, 1)
        setphase!(cs, ALMaSS.marchfirst)
        cs.growingdegreedays = 0.0
    end
    # 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 = cs.croptype.mingrowthtemp
    ismissing(basetemp) && (basetemp = 5.0)
    delta_gdd = bounds((maxtemp(model) + mintemp(model)) / 2 - basetemp)
    cs.growingdegreedays += delta_gdd
    # update crop growth
    growcrop!(cs, delta_gdd, model)
end


## CROP MANAGEMENT AND GROWTH FUNCTIONS

"""
    sow!(cropstate, model, cropname)

Change the cropstate to sow the specified crop.
"""
function sow!(cs::CropState, model::SimulationModel, cropname::String)
    !ismissing(cs.croptype.minsowdate) && model.date < cs.croptype.minsowdate &&
        @warn "$(model.date) is earlier than the minimum sowing date for $(cropname)."
    cs.croptype = model.crops[cropname]

    setphase!(cs, ALMaSS.sow)
    cs.growingdegreedays = 0.0
    # cs.phase = ALMaSS.sow
    # cs.growingdegreedays = 0.0
    # # TODO: set from crop curve params
    # #       cs.height = cs.croptype.lownutrientgrowth.height[sow]  # or highnutrientgrowth
    # cs.height = 0.0cm
    # cs.LAItotal = 0.0
    # cs.LAIgreen = 0.0

    cs.mature = false
    cs.events = Vector{Management}()
end

"""
    harvest!(cropstate, model)

Harvest the crop of this cropstate.
"""
function harvest!(cs::CropState, model::SimulationModel)
    if cs.phase in (ALMaSS.harvest1, ALMaSS.harvest2)
        setphase!(cs, ALMaSS.harvest2)
    else
        setphase!(cs, ALMaSS.harvest1)
    end
    cs.mature = false

    #TODO calculate and return yield
end

#TODO fertilise!()
#TODO spray!()
#TODO till!()

"""
    growcrop!(cropstate, delta_gdd, model)

Update the crop state on this farm plot for a given daily change in
growing-degree days `delta_gdd` (units of 'Kelvin Days').

Implements the ALMaSS crop growth model by Topping et al. (see
`ALMaSS/Landscape/plant.cpp:PlantGrowthData::FindDiff()` and
`ALMaSS/Landscape/elements.cpp:VegElement::DoDevelopment()`).
"""
function growcrop!(cs::CropState, delta_gdd::Real, model::SimulationModel)
    total_gdd = cs.growingdegreedays
    # TODO: logic the wrong way around for fertiliser?
    fertiliser in cs.events ?
        curve = cs.croptype.lownutrientgrowth :
        curve = cs.croptype.highnutrientgrowth
    if ismissing(curve)
         #FIXME what to do if the curve is empty?
        error("Missing curve")
    end
    points = curve.GDD[cs.phase]

    if length(points) == 1 && points[1] == 99999.0
        !(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
    end

    # TODO: calc oldidx, catch cases where the delta_gdd crosses more
    #       than one inflection points
    # TODO: check logic for idx
    idx = findfirst(x -> x > total_gdd, points)
    if isnothing(idx)
        idx = lastindex(points)
    elseif idx != firstindex(points)
        idx -= 1  # we need the last index that is smaller than total_gdd
    end
    cs.height += curve.height[cs.phase][idx] * delta_gdd
    cs.LAItotal += curve.LAItotal[cs.phase][idx] * delta_gdd
    cs.LAIgreen += curve.LAIgreen[cs.phase][idx] * delta_gdd

    function lastphase(cs::CropState)
        if  cs.croptype.name == "maize"
            return sow
        else
            return marchfirst
        end
    end

    # set mature if it's the last growth point of the last phase
    if idx == lastindex(points) && cs.phase == lastphase(cs)
        cs.mature = true
    end
end

end # module ALMaSS