### 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 landscape event that can be simulated" @enum EventType 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 fieldid::Union{Missing,Int64} events::Vector{EventType} end """ FarmEvent A data structure to define a landscape event, giving its type, spatial extent, and duration. """ mutable struct FarmEvent type::EventType pixels::Vector{Tuple{Int64,Int64}} duration::Int64 end """ initlandscape(landcovermap, farmfieldsmap) Initialise the model landscape based on the map files specified in the configuration. Returns a matrix of pixels. """ function initlandscape(landcovermap::String, farmfieldsmap::String) @debug "Initialising landscape" #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, Vector{Symbol}()) end end return landscape end """ updateevents!(model) Cycle through the list of events, removing those that have expired. """ function updateevents!(model::AgentBasedModel) 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.type, 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::AgentBasedModel, pixels::Vector{Tuple{Int64,Int64}}, name::EventType, 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::AgentBasedModel) 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::AgentBasedModel) ismissing(model.landscape[pos...].fieldid) ? nothing : model[model.landscape[pos...].fieldid] 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::AgentBasedModel, habitatdescriptor::Function) (habitatdescriptor(pos, model)) && (return 0) dist = 1 width, height = size(model.landscape) while dist <= width || dist <= height # 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)) for x in minx:maxx (y1 > 0 && habitatdescriptor((x,y1), model)) && (return dist) (y2 <= height && habitatdescriptor((x,y2), model)) && (return dist) 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 dist) (x2 <= width && habitatdescriptor((x2,y), model)) && (return dist) end dist += 1 end return Inf 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::AgentBasedModel, habitattype::LandCover) #XXX allow testing for multiple habitat types? # can't use @habitat here because nature.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::AgentBasedModel) # can't use @habitat here because nature.jl is loaded later than this file distanceto(pos, model, function(p,m) landcover(p,m) != landcover(pos, model) end) end