landscape.jl

### Persefone.jl - a model of agricultural landscapes and ecosystems in Europe.
###
### This file manages the landscape maps that underlie the model.
###

## IMPORTANT: do not change the order of this enum, or initlandscape() will break!
"The land cover classes encoded in the Mundialis Sentinel data."
@enum LandCover nodata forest grass water builtup soil agriculture

"The types of management event that can be simulated"
@enum Management tillage sowing fertiliser pesticide harvesting

"""
    Pixel

A pixel is a simple data structure to combine land use and ownership information
in a single object. The model landscape consists of a matrix of pixels.
(Note: further landscape information may be added here in future.)
"""
mutable struct Pixel
    landcover::LandCover          # land cover class at this position
    fieldid::Union{Missing,Int64} # ID of the farmplot (if any) at this position
    events::Vector{Management}    # management events that have been applied to this pixel
    animals::Vector{Int64}        # IDs of animals currently at this position
    territories::Vector{Int64}    # IDs of animals that claim this pixel as part of their territory
end

Pixel(landcover::LandCover, fieldid::Union{Missing,Int64}) =
    Pixel(landcover, fieldid, Vector{Management}(), Vector{Int64}(), Vector{Int64}())
Pixel(landcover::LandCover) = Pixel(landcover, missing)

"""
    FarmEvent

A data structure to define a landscape event, giving its type,
spatial extent, and duration.
"""
mutable struct FarmEvent
    management::Management
    pixels::Vector{Tuple{Int64,Int64}}
    duration::Int64
end

"""
    initlandscape(directory, landcovermap, farmfieldsmap)

Initialise the model landscape based on the map files specified in the
configuration. Returns a matrix of pixels.
"""
function initlandscape(directory::String, landcovermap::String, farmfieldsmap::String)
    @debug "Initialising landscape"
    landcovermap = joinpath(directory, landcovermap)
    farmfieldsmap = joinpath(directory, farmfieldsmap)
    #TODO replace errors with exception
    !(isfile(landcovermap)) && Base.error("Landcover map $(landcovermap) doesn't exist.")
    !(isfile(farmfieldsmap)) && Base.error("Field map $(farmfieldsmap) doesn't exist.")
    #XXX can I read the files in using FileIO, to avoid the GeoArrays dependency?
    # read in TIFFs, replacing missing values with 0
    landcover = coalesce(GeoArrays.read(landcovermap), 0)
    farmfields = coalesce(GeoArrays.read(farmfieldsmap), 0)
    (size(landcover) != size(farmfields)) && Base.error("Input map sizes don't match.")
    width, height = size(landcover)[1:2]
    landscape = Matrix{Pixel}(undef, width, height)
    for x in 1:width
        for y in 1:height
            # convert the numeric landcover value to LandCover, then create the Pixel object
            lcv = LandCover(Int(landcover[x,y][1]/10))
            ff = Int64(farmfields[x,y][1])
            (iszero(ff)) && (ff = missing)
            landscape[x,y] = Pixel(lcv, ff)
        end
    end
    return landscape
end

"""
    updateevents!(model)

Cycle through the list of events, removing those that have expired.
"""
function updateevents!(model::SimulationModel)
    expiredevents = []
    for e in 1:length(model.events)
        event = model.events[e]
        event.duration -= 1
        if event.duration <= 0
            push!(expiredevents, e)
            for p in event.pixels
                i = findnext(x -> x == event.management, model.landscape[p...].events, 1)
                deleteat!(model.landscape[p...].events, i)
            end
        end
    end
    deleteat!(model.events, expiredevents)
end

"""
    createevent!(model, pixels, name, duration=1)

Add a farm event to the specified pixels (a vector of position tuples) for a given duration.
"""
function createevent!(model::SimulationModel, pixels::Vector{Tuple{Int64,Int64}},
                     name::Management, duration::Int64=1)
    push!(model.events, FarmEvent(name, pixels, duration))
    for p in pixels
        push!(model.landscape[p...].events, name)
    end
end

"""
    landcover(position, model)

Return the land cover class at this position (utility wrapper).
"""
function landcover(pos::Tuple{Int64,Int64}, model::SimulationModel)
    model.landscape[pos...].landcover
end

"""
    farmplot(position, model)

Return the farm plot at this position, or nothing if there is none (utility wrapper).
"""
function farmplot(pos::Tuple{Int64,Int64}, model::SimulationModel)
    id = model.landscape[pos...].fieldid
    ismissing(id) ? nothing : model.farmplots[id]
end

"""
    directionto(pos, model, habitatdescriptor)

Calculate the direction from the given location to the closest location matching the
habitat descriptor function. Returns a coordinate tuple (target - position), or nothing
if no matching habitat is found. Caution: can be computationally expensive!
"""
function directionto(pos::Tuple{Int64,Int64}, model::SimulationModel, habitatdescriptor::Function)
    #XXX split this up into a findnearest() and a directionto() function?
    # search in expanding rectangles around the origin
    (habitatdescriptor(pos, model)) && (return (0,0))
    dist = 1
    width, height = size(model.landscape)
    while dist <= width || dist <= height #XXX is there a quicker method?
        # check the upper and lower bounds of the enclosing rectangle
        y1 = pos[2] - dist
        y2 = pos[2] + dist
        minx = maximum((pos[1]-dist, 1))
        maxx = minimum((pos[1]+dist, width))
        #XXX alternately, we could check individually whether y1 and y2 are in bounds, and loop
        # over x for each separately (not sure if that is quicker, but it may be)
        for x in minx:maxx
            (y1 > 0 && habitatdescriptor((x,y1), model)) && (return (x, y1).-pos)
            (y2 <= height && habitatdescriptor((x,y2), model)) && (return (x, y2).-pos)
        end
        # check the left and right bounds of the enclosing rectangle
        x1 = pos[1] - dist
        x2 = pos[1] + dist
        miny = maximum((pos[2]-dist+1, 1))
        maxy = minimum((pos[2]+dist-1, height))
        for y in miny:maxy
            (x1 > 0 && habitatdescriptor((x1,y), model)) && (return (x1,y).-pos)
            (x2 <= width && habitatdescriptor((x2,y), model)) && (return (x2,y).-pos)
        end
        dist += 1
    end
    return nothing
end

"""
    directionto(pos, model, habitattype)

Calculate the direction from the given location to the closest habitat of the specified type.
Returns a coordinate tuple (target - position), or nothing if no matching habitat is found.
Caution: can be computationally expensive!
"""
function directionto(pos::Tuple{Int64,Int64}, model::SimulationModel, habitattype::LandCover)
    # can't use @habitat here because macros.jl is loaded later than this file
    directionto(pos, model, function(p,m) landcover(p,m) == habitattype end)
end

"""
    distanceto(pos, model, habitatdescriptor)

Calculate the distance from the given location to the closest location matching the
habitat descriptor function. Caution: can be computationally expensive!
"""
function distanceto(pos::Tuple{Int64,Int64}, model::SimulationModel, habitatdescriptor::Function)
    #XXX this is very imprecise, as diagonal distances are not calculated trigonometrically
    target = directionto(pos, model, habitatdescriptor)
    isnothing(target) && return Inf
    return maximum(abs.(target))*@param(world.mapresolution)
end

"""
    distanceto(pos, model, habitattype)

Calculate the distance from the given location to the closest habitat of the specified type.
Caution: can be computationally expensive!
"""
function distanceto(pos::Tuple{Int64,Int64}, model::SimulationModel, habitattype::LandCover)
    # can't use @habitat here because macros.jl is loaded later than this file
    distanceto(pos, model, function(p,m) landcover(p,m) == habitattype end)
end

"""
    distancetoedge(pos, model)

Calculate the distance from the given location to the closest neighbouring habitat.
Caution: can be computationally expensive!
"""
function distancetoedge(pos::Tuple{Int64,Int64}, model::SimulationModel)
    # can't use @habitat here because macros.jl is loaded later than this file
    distanceto(pos, model, function(p,m) landcover(p,m) != landcover(pos, model) end)
end

"""
    randompixel(position, model, range, habitatdescriptor)

Find a random pixel within a given `range` of the `position` that matches the
habitatdescriptor (create this using [`@habitat`](@ref)).
"""
function randompixel(pos::Tuple{Int64,Int64}, model::SimulationModel, range::Length,
                     habitatdescriptor::Function=(pos,model)->nothing)
    range = Int(floor(range / @param(world.mapresolution)))
    for x in @shuffle!(collect((pos[1]-range):(pos[1]+range)))
        for y in @shuffle!(collect((pos[2]-range):(pos[2]+range)))
            !inbounds((x,y), model) && continue
            habitatdescriptor((x,y), model) && return (x,y)
        end
    end
    nothing
end

"""
    randomdirection(model, distance)

Get a random direction coordinate tuple within the specified distance.
"""
function randomdirection(model::SimulationModel, distance::Length)
    range = Int(floor(distance / @param(world.mapresolution)))
    Tuple(@rand(-range:range, 2))
end

"""
    bounds(x; max=Inf, min=0)

A utility function to make sure that a number is within a given set of bounds.
Returns `max`/`min` if `x` is greater/less than this.
"""
function bounds(x::Number; max::Number=Inf, min::Number=0)
    x > max ? max :
        x < min ? min :
        x
end

"""
    inbounds(pos, model)

Is the given position within the bounds of the model landscape?
"""
function inbounds(pos, model)
    dims = size(model.landscape)
    pos[1] > 0 && pos[1] <= dims[1] && pos[2] > 0 && pos[2] <= dims[2]    
end

"""
    safebounds(pos, model)

Make sure that a given position is within the bounds of the model landscape.
"""
function safebounds(pos::Tuple{Int64,Int64}, model::SimulationModel)
    width, height = size(model.landscape)
    (bounds(pos[1], max=width, min=1), bounds(pos[2], max=height, min=1))
end