farmplot.jl

### Persefone.jl - a model of agricultural landscapes and ecosystems in Europe.
###
### 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.
"""
mutable struct FarmPlot <: ModelAgent
    #TODO add Unitful
    const id::Int64
    pixels::Vector{Tuple{Int64, Int64}}
    croptype::CropType
    phase::GrowthPhase
    growingdegreedays::Float64
    height::Float64
    LAItotal::Float64
    LAIgreen::Float64
    #biomass::Float64 #XXX I need to figure out how to calculate this
    events::Vector{Management}
end

"""
    initfields!(model)

Initialise the model with its farm plots.
"""
function initfields!(model::SimulationModel)
    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.farmplots[objectid].pixels, (x,y))
            else
                #XXX does this phase calculation work out?
                month(model.date) < 3 ? phase = janfirst : phase = marchfirst
                fp = FarmPlot(length(model.farmplots)+1, [(x,y)],
                                model.crops["natural grass"], phase,
                              0.0, 0.0, 0.0, 0.0, Vector{Management}())
                push!(model.farmplots, fp)
                model.landscape[x,y].fieldid = fp.id
                convertid[rawid] = fp.id
                n += 1
            end
        end
    end
    @info "Initialised $n farm plots."
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 = 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::SimulationModel)
    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::SimulationModel)
    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::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) ? 0cm : #FIXME should not return 0
              model.farmplots[model.landscape[pos...].fieldid].height
end

"""
    cropcover(model, position)

Return the percentage ground cover of the crop at this position, or nothing if there is no crop
here (utility wrapper).
"""
function cropcover(pos::Tuple{Int64,Int64}, model::SimulationModel)
    #FIXME LAItotal != ground cover?
    ismissing(model.landscape[pos...].fieldid) ? 0 : #FIXME should not return 0
              model.farmplots[model.landscape[pos...].fieldid].LAItotal
end