aquacrop.jl

# Persefone.jl wrapper for AquaCrop.jl

module AquaCropWrapper

const CROPFILE = "crop_data.csv"

import AquaCrop
import CSV
using Dates: Date
using DataFrames: DataFrame

using Persefone:
    AbstractCropState,
    AbstractCropType,
    AnnualDate,
    cm,
    daynumber,
    SoilType,
    SimulationModel

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

using Unitful: @u_str

# We can't use Length{Float64} as that is a Union type
const Tlength = typeof(1.0cm)

# TODO: read crop names directly from AquaCrop.jl
# names extracted by running this command in the AquaCrop.jl src dir:
#     grep -nir 'crop_type ==' src/ | cut -d '=' -f 3 | sed 's/$/,/' | sort | uniq
const AQUACROP_CROPNAMES = [
    "alfalfaGDD",
    "barley",
    "barleyGDD",
    "cotton",
    "cottonGDD",
    "drybean",
    "drybeanGDD",
    "maize",
    "maizeGDD",
    "oat",
    "paddyrice",
    "paddyriceGDD",
    "potato",
    "potatoGDD",
    "quinoa",
    "rapeseed",
    "sorghum",
    "sorghumGDD",
    "soybean",
    "soybeanGDD",
    "sugarbeet",
    "sugarbeetGDD",
    "sugarcane",
    "sunflower",
    "sunflowerGDD",
    "tef",
    "tomato",
    "tomatoGDD",
    "wheat",
    "wheatGDD",
]

@kwdef struct CropType <: AbstractCropType
    name::String
    aquacrop_name::String
    group::String
    minsowdate::Union{Missing,AnnualDate}
    maxsowdate::Union{Missing,AnnualDate}
end

cropname(ct::CropType) = ct.name

"""
    readcropparameters(cropdirectory)

Parse a CSV file containing some parameters required to map AquaCrop
crop names to Persefone crop names as well as additional crop data
needed for Persefone (cropgroup, minsowdate, maxsowdate).
"""
function readcropparameters(cropdirectory::String)
    @debug "Reading extra crop parameters for AquaCrop crop model"
    croptypes = Dict{String, AbstractCropType}()
    df = CSV.read(joinpath(cropdirectory, CROPFILE), DataFrame;
                  missingstring="NA", types=[String, String, String, AnnualDate, AnnualDate])
    for row in eachrow(df)
        croptypes[row.persefone_name] =
            CropType(name = row.persefone_name,
                     aquacrop_name = row.aquacrop_name,
                     group = row.crop_group,
                     minsowdate = row.min_sow_date,
                     maxsowdate = row.max_sow_date)
    end
    return croptypes
end

function make_aquacrop_cropfield(croptype::CropType,
                                 model::SimulationModel,
                                 soiltype::SoilType)
    # TODO: get this mapping for soil types from a CSV file?
    soiltype_to_aquacrop = Dict(soiltype => replace(string(soiltype), "_" => " ")
                                for soiltype in instances(SoilType))
    aquacrop_soiltype = soiltype_to_aquacrop[soiltype]

    start_date = model.date            # TODO: maybe `today(model)`
    end_date = model.weather.lastdate  # TODO: maybe `lastdate(model)`
    start_daynum = daynumber(model.weather, start_date)
    end_daynum = daynumber(model.weather, end_date)
    dayrange = start_daynum:end_daynum
    input_Tmin = @view model.weather.mintemp[dayrange]
    input_Tmax = @view model.weather.maxtemp[dayrange]
    input_ETo = @view model.weather.evapotranspiration[dayrange]
    input_Rain = @view model.weather.precipitation[dayrange]

    # Generate the keyword object for the AquaCrop simulation
    kwargs = (
        ## Necessary keywords
        # runtype
        runtype = AquaCrop.NoFileRun(),
        # project input
        Simulation_DayNr1 = start_date,
        Simulation_DayNrN = end_date,
        Crop_Day1 = start_date,    # TODO: originally start_date + Week(1)
        Crop_DayN = end_date,
        # soil
        soil_type = aquacrop_soiltype,
        # crop
        crop_type = croptype.aquacrop_name,
        # Climate
        InitialClimDate = start_date,
        ## Optional keyworkds
        # Climate
        Tmin = input_Tmin,
        Tmax = input_Tmax,
        ETo = input_ETo,
        Rain = input_Rain,
        # change soil properties
        soil_layers = Dict("Thickness" => 5.0)
    )

    aquacrop_cropfield, all_ok = AquaCrop.start_cropfield(; kwargs...)
    if ! all_ok.logi
        @error "Failed calling AquaCrop.start_cropfield()\nAquaCrop error: $(all_ok.msg)"
    end

    return aquacrop_cropfield
end

mutable struct CropState <: AbstractCropState
    croptype::CropType
    height::Tlength  # TODO: remove height field, supply from cropstate
    soiltype::SoilType
    cropstate::AquaCrop.AquaCropField

    function CropState(croptype::CropType, soiltype::SoilType,
                       model::SimulationModel, height::Tlength=0.0cm)
        aquacrop_cropfield = make_aquacrop_cropfield(croptype, model, soiltype)
        return new(croptype, height, soiltype, aquacrop_cropfield)
    end
end

croptype(cs::CropState) = cs.croptype
cropname(cs::CropState) = cropname(croptype(cs))
cropheight(cs::CropState) = cs.height  # TODO: calculate from AquaCrop state info
cropcover(cs::CropState) = AquaCrop.canopycover(cs.cropstate)
cropyield(cs::CropState) = AquaCrop.dryyield(cs.cropstate)  # TODO: there is also freshyield
function isharvestable(cs::CropState)
    stage = get_aquacrop_stage(cs)
    return !ismissing(stage) && stage == 0
end
makecropstate(croptype::CropType, model::SimulationModel, soiltype::SoilType) =
    CropState(croptype, soiltype, model)
function setsoiltype!(cs::CropState, soiltype::SoilType)
    # TODO: this does not affect the actual soil type of the state of
    # the aquacrop simulation
    cs.soiltype = soiltype
    return nothing
end

function get_aquacrop_stage(cs::CropState)
    stages = cs.cropstate.dayout.Stage
    length(stages) > 0 ? last(stages) : missing
end

"""
    stepagent!(cropstate, model)

Update a crop state by one day.
"""
function stepagent!(cs::CropState, model::SimulationModel)
    AquaCrop.dailyupdate!(cs.cropstate)
end

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

Change the cropstate to sow the specified crop.
"""
function sow!(cs::CropState, model::SimulationModel, cropname::String)
    if cropname ∉ keys(model.crops)
        @error "cropname \"$cropname\" not found"
    end
    new_croptype = model.crops[cropname]
    if typeof(new_croptype) !== CropType
        error("Cannot change crop model of AquaCrop crop state.")
    end

    cs.croptype = new_croptype
    cs.height = 0.0cm
    # cs.soiltype stays the way it is
    cs.cropstate = make_aquacrop_cropfield(cs.croptype, model, soiltype)
end

"""
    harvest!(cropstate, model)

Harvest the crop of this cropstate.
"""
function harvest!(cs::CropState, model::SimulationModel)
    # TODO: implement this
    return 0.0u"g/m^2"
end

end # module