diff --git a/src/Persefone.jl b/src/Persefone.jl index 7e2fdd5c24efa789fa4259f2c203066f1852d85b..af2b38fdabe09670a594e93efe37bab614e1c2dc 100644 --- a/src/Persefone.jl +++ b/src/Persefone.jl @@ -29,14 +29,14 @@ using Random, Serialization, StableRNGs, - TOML + TOML, + Unitful ## Packages that may be useful later on: # MacroTools, http://fluxml.ai/MacroTools.jl/stable/utilities/ # Debugger, https://github.com/JuliaDebug/Debugger.jl # PackageCompiler, https://julialang.github.io/PackageCompiler.jl/stable/ # SpatialEcology, https://github.com/EcoJulia/SpatialEcology.jl -# OmniScape, https://circuitscape.org/ # PlantSimEngine, https://virtualplantlab.github.io/PlantSimEngine.jl/stable/ # Overseer, https://juliapackages.com/p/overseer -> ECS @@ -99,6 +99,14 @@ export savemodelobject, loadmodelobject +## Import and define units and dimensions +import Unitful: cm, m, km, ha, Length, Area +const m² = m^2 +const km² = km^2 +import Base./ # enable division with different length/area unit types +/(x::S,y::T) where {S<:Length, T<:Length} = (upreferred(x)/m) / (upreferred(y)/m) +/(x::S,y::T) where {S<:Area, T<:Area} = (upreferred(x)/m²) / (upreferred(y)/m²) + """ SimulationModel diff --git a/src/core/input.jl b/src/core/input.jl index d7011fccb6bfac37adaec1f9d635d47ffd02743c..5c22b60bf3887d43d9aaaff63e2956d104efb32c 100644 --- a/src/core/input.jl +++ b/src/core/input.jl @@ -87,6 +87,7 @@ function preprocessparameters(settings::Dict{String,Any}, defaultoutdir::String) if settings["core.startdate"] > settings["core.enddate"] Base.error("Enddate is earlier than startdate.") #TODO replace with exception end + settings["world.mapresolution"] = settings["world.mapresolution"] * 1m #FIXME enable parallelisation # if !isempty(settings["internal.scanparams"]) && (nprocs() < 2) # @warn "To parallelise multiple simulation runs, use `julia -p <n>`." diff --git a/src/core/output.jl b/src/core/output.jl index 9842feff3606def29cb43ec229c9596699a73537..5028ddee623ea5006461bd0ef309b9c5ad404313 100644 --- a/src/core/output.jl +++ b/src/core/output.jl @@ -94,7 +94,7 @@ function saveinputfiles(model::SimulationModel) println(f, "# WARNING: Your repository contains uncommitted changes. This may") println(f, "# compromise the reproducibility of this simulation run.\n") end - TOML.print(f, prepareTOML(model.settings)) + TOML.print(f, prepareTOML(deepcopy(model.settings))) end # Copy the map files to the output folder lcmap = joinpath(@param(world.mapdirectory), @param(world.landcovermap)) @@ -118,6 +118,7 @@ function prepareTOML(settings) settings["core.loglevel"] == Logging.Debug ? settings["core.loglevel"] = "debug" : settings["core.loglevel"] == Logging.Warn ? settings["core.loglevel"] = "warn" : settings["core.loglevel"] = "info" + settings["world.mapresolution"] = settings["world.mapresolution"] / m # convert dict structure fulldict = Dict{String, Dict{String, Any}}() for parameter in keys(settings) diff --git a/src/crop/crops.jl b/src/crop/crops.jl index 927b2373cee95af9275301ad430a1b8b7df9d35e..996c6dabd0187189e6dd701e62669598003f4f63 100644 --- a/src/crop/crops.jl +++ b/src/crop/crops.jl @@ -19,6 +19,7 @@ seasonal dates or agricultural events. The values in this struct define one crop growth curve. """ struct CropCurveParams + #TODO add Unitful curveID::Int highnutrients::Bool GDD::Dict{GrowthPhase,Vector{Float64}} diff --git a/src/crop/farmplot.jl b/src/crop/farmplot.jl index cd18822843a68bbaf3e4c8511e2c6cb2362a5282..3f037e09a49bea0519410d778bb8f39197c5df5d 100644 --- a/src/crop/farmplot.jl +++ b/src/crop/farmplot.jl @@ -12,6 +12,7 @@ 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 diff --git a/src/nature/populations.jl b/src/nature/populations.jl index 7494dbc3ebab3bdea94f16f941fb27ddad7f643a..4292ba142163172b8882b0542b18fb1e6e19f57c 100644 --- a/src/nature/populations.jl +++ b/src/nature/populations.jl @@ -22,11 +22,11 @@ using [`@populate`](@ref). - `popsize` determines the number of individuals that will be created, dispersed over the suitable locations in the landscape. If this is zero or negative, one individual will be created in every suitable location. If it is greater than the number of suitable - locations, multiple individuals will be created per location. Alternately, use `popdensity`. + locations, multiple individuals will be created per location. Alternately, use `indarea`. -- `popdensity`: if this is greater than zero, the chance of creating an individual (or - pair of individuals) at a suitable location is 1/popdensity. Use this as an alternative - to `popsize`. +- `indarea`: if this is greater than zero, it determines the habitat area allocated to each + individual or pair. To be precise, the chance of creating an individual (or pair of + individuals) at a suitable location is 1/indarea. Use this as an alternative to `popsize`. - If `pairs` is true, a male and a female individual will be created in each selected location, otherwise, only one individual will be created at a time. (default: false) @@ -40,7 +40,7 @@ using [`@populate`](@ref). birthphase::Function habitat::Function = @habitat(true) popsize::Int64 = -1 - popdensity::Int64 = -1 #XXX this is counterintuitive + indarea::Area = -1m² pairs::Bool = false asexual::Bool = false end @@ -74,21 +74,24 @@ This is an internal function called by initpopulation!(), and was split off from better testing. """ function initpopulation!(species::Type, p::PopInitParams, model::SimulationModel) - (p.popsize <= 0 && p.popdensity <= 0) && # can be legit if a habitat descriptor is provided - @warn("initpopulation!() called with popsize and popdensity both <= 0") - (p.popsize > 0 && p.popdensity > 0) && #XXX not sure what this would do - @warn("initpopulation!() called with popsize and popdensity both > 0") + (p.popsize <= 0 && p.indarea <= 0m²) && # can be legit if a habitat descriptor is provided + @warn("initpopulation!() called with popsize and indarea both <= 0") + (p.popsize > 0 && p.indarea > 0m²) && #XXX not sure what this would do + @warn("initpopulation!() called with popsize and indarea both > 0") # create as many individuals as necessary in the landscape n = 0 lastn = 0 width, height = size(model.landscape) + if p.indarea > 0m² + pixelsperind = Int(round(p.indarea / @param(world.mapresolution)^2)) + end while n == 0 || n < p.popsize for x in @shuffle!(Vector(1:width)) for y in @shuffle!(Vector(1:height)) if p.habitat((x,y), model) && - (p.popdensity <= 0 || n == 0 || @chance(1/p.popdensity)) #XXX what if pd==0? + (p.indarea <= 0m² || n == 0 || @chance(1/pixelsperind)) #XXX what if ppi==0? #XXX `n==0` above guarantees that at least one individual is created, even - # in a landscape that is otherwise too small for the specified popdensity - + # in a landscape that is otherwise too small for the specified indarea - # do we want this? if p.pairs a1 = species(length(model.animals)+1, male, (-1, -1), (x,y), p.initphase) @@ -202,9 +205,10 @@ end Return a list of IDs of the animals within a given radius of the position. """ -function nearby_ids(pos::Tuple{Int64,Int64}, model::SimulationModel, radius::Int64) +function nearby_ids(pos::Tuple{Int64,Int64}, model::SimulationModel, radius::Length) ids = [] msize = size(model.landscape) + radius = Int(floor(radius / @param(world.mapresolution))) for x in (pos[1]-radius):(pos[1]+radius) (x < 1 || x > msize[1]) && continue for y in (pos[2]-radius):(pos[2]+radius) @@ -222,7 +226,7 @@ end Return a list of animals in the given radius around this position, optionally filtering by species. """ function nearby_animals(pos::Tuple{Int64,Int64}, model::SimulationModel; - radius::Int64=0, species="") #XXX add type for species + radius::Length=0m, species="") #XXX add type for species neighbours = nearby_ids(pos, model, radius) isempty(neighbours) && return neighbours if isempty(species) @@ -239,7 +243,7 @@ Return the number of animals in the given radius around this position, optionall by species. """ function countanimals(pos::Tuple{Int64,Int64}, model::SimulationModel; - radius::Int64=0, species="") #XXX add type for species + radius::Length=0m, species="") #XXX add type for species length(nearby_animals(pos, model, radius=radius, species=species)) end @@ -249,7 +253,8 @@ end Return a list of animals in the given radius around this animal, excluding itself. By default, only return conspecific animals. """ -function neighbours(animal::Animal, model::SimulationModel, radius::Int64=0, conspecifics::Bool=true) +function neighbours(animal::Animal, model::SimulationModel, + radius::Length=0m, conspecifics::Bool=true) filter(a -> a.id != animal.id, nearby_animals(animal.pos, model, radius = radius, species = conspecifics ? speciesof(animal) : "")) @@ -272,7 +277,8 @@ Calculate the distance from the given position to the animal. """ function distanceto(pos::Tuple{Int64,Int64}, model::SimulationModel, animal::Animal) # have to use a coordinate as first argument rather than an animal because of @distanceto - maximum(abs.(animal.pos .- pos)) + #XXX this is very imprecise because diagonal distances are not calculated trigonometrically + maximum(abs.(animal.pos .- pos)) * @param(world.mapresolution) end """ @@ -281,9 +287,10 @@ end Move the follower animal to a location near the leading animal. """ function followanimal!(follower::Animal, leader::Animal, model::SimulationModel, - distance::Int64=0) + distance::Length=0m) #TODO test function - spacing = Tuple(@rand(-distance:distance, 2)) + dist = Int(floor(distance / @param(world.mapresolution))) + spacing = Tuple(@rand(-dist:dist, 2)) targetposition = safebounds(spacing .+ leader.pos, model) move!(follower, model, targetposition) end @@ -303,12 +310,13 @@ function move!(animal::Animal, model::SimulationModel, position::Tuple{Int64,Int end """ - walk!(animal, model, direction, steps=1) + walk!(animal, model, direction, distance=1pixel) -Let the animal move a given number of step in the given direction ("north", "northeast", +Let the animal move a given number of steps in the given direction ("north", "northeast", "east", "southeast", "south", "southwest", "west", "northwest", "random"). """ -function walk!(animal::Animal, model::SimulationModel, direction::String, steps=1) +function walk!(animal::Animal, model::SimulationModel, direction::String, distance::Length=-1m) + distance < 0m ? steps = 1 : steps = Int(floor(distance/@param(world.mapresolution))) if direction == "north" shift = (0,-steps) elseif direction == "northeast" @@ -326,30 +334,33 @@ function walk!(animal::Animal, model::SimulationModel, direction::String, steps= elseif direction == "northwest" shift = (-steps,-steps) elseif direction == "random" - shift = Tuple(@rand([-steps,steps], 2)) + shift = Tuple(@rand([-steps,0,steps], 2)) else @error "Invalid direction in @walk: "*direction end - move!(animal, model, animal.pos .+ shift) + walk!(animal, model, shift) end """ - walk!(animal, model, direction, speed=-1) + walk!(animal, model, direction, distance=-1) Let the animal move in the given direction, where the direction is defined by an (x, y) tuple to specify the shift in coordinates. -If speed >= 0, walk no more than the given number of cells. -""" -function walk!(animal::Animal, model::SimulationModel, direction::Tuple{Int64,Int64}, speed::Int64=-1) - #TODO add a habitat descriptor? - if speed >= 0 - direction[1] > speed && (direction[1] = speed) - direction[2] > speed && (direction[2] = speed) - direction[1] < -speed && (direction[1] = -speed) - direction[2] < -speed && (direction[2] = -speed) +If maxdist >= 0, move no further than the specified distance. +""" +function walk!(animal::Animal, model::SimulationModel, direction::Tuple{Int64,Int64}, + maxdist::Length=-1m) + #TODO test + distance = Int(floor(maxdist/@param(world.mapresolution))) + if distance >= 0 + direction[1] > distance && (direction[1] = distance) + direction[2] > distance && (direction[2] = distance) + direction[1] < -distance && (direction[1] = -distance) + direction[2] < -distance && (direction[2] = -distance) end newpos = animal.pos .+ direction move!(animal, model, newpos) end +#TODO add a walk function with a habitat descriptor ##TODO add walktoward or similar function (incl. pathfinding?) diff --git a/src/nature/species/skylark.jl b/src/nature/species/skylark.jl index 72e410247373dc431ca1402741ece889642219a4..b639630b22f5a3b891308e1f04b69a7286a5e21d 100644 --- a/src/nature/species/skylark.jl +++ b/src/nature/species/skylark.jl @@ -7,7 +7,7 @@ skylarkhabitat = @habitat((@landcover() == grass || # settle on grass or arable land (but not maize) (@landcover() == agriculture && @cropname() != "maize")) && - @distancetoedge() > 5) # at least 50m from other habitats + @distancetoedge() > 50m) # at least 50m from other habitats #XXX this ought to check for distance to forest and builtup, # but that's very expensive (see below) # @distanceto(forest) > 5 && # at least 50m from forest edges @@ -31,8 +31,6 @@ At the moment, this implementation is still in development. ISBN 3-89104-019-9 """ @species Skylark begin - #XXX use Unitful.jl - # species parameters const eggtime::Int64 = 11 # 11 days from laying to hatching const nestlingtime::UnitRange{Int64} = 7:11 # 7-11 days from hatching to leaving nest @@ -57,11 +55,11 @@ At the moment, this implementation is still in development. # individual variables daystonextphase::Int64 = 0 # days remaining until fledging or maturity migrationdates::Tuple = () # is defined by each individual in @create(Skylark) + territory::Vector = [] # pixels that this skylark claims as its territory mate::Int64 = -1 # the agent ID of the mate (-1 if none) nest::Tuple = () # coordinates of current nest nestcompletion::Int64 = 0 # days left until the nest is built clutch::Vector{Int64} = Vector{Int64}() # IDs of offspring in current clutch - end """ @@ -101,7 +99,7 @@ check mortality. @phase Skylark fledgling begin #TODO add feeding & growth @kill(self.fledglingpredationmortality, "predation") - @walk("random") #TODO add movement following the parents + @walk("random", 10m) #TODO add movement following the parents # if self.age == self.fledglingtime+self.eggtime # @kill(self.firstyearmortality, "first year mortality") #XXX mechanistic? # @setphase(nonbreeding) @@ -120,15 +118,15 @@ As a non-breeding adult, move around with other individuals and check for migrat @phase Skylark nonbreeding begin # flocking behaviour - follow a random neighbour or move randomly #TODO add feeding and mortality, respect habitat when moving - neighbours = @neighbours(10) #XXX magic number + neighbours = @neighbours(100m) #XXX magic number isempty(neighbours) ? - @walk("random", 5) : - @follow(@rand(neighbours), 2) + @walk("random", 50m) : + @follow(@rand(neighbours), 20m) # check if the bird migrates leave, arrive = self.migrationdates - m, d = monthday(model.date) - migrate = (((m < arrive[1]) || (m == arrive[1] && d < arrive[2])) || - ((m > leave[1]) || (m == leave[1] && d >= leave[2]))) + month, day = monthday(model.date) + migrate = (((month < arrive[1]) || (month == arrive[1] && day < arrive[2])) || + ((month > leave[1]) || (month == leave[1] && day >= leave[2]))) if migrate #FIXME not all migrate? @kill(self.migrationmortality, "migration") returndate = Date(year(model.date)+1, arrive[1], arrive[2]) @@ -149,14 +147,14 @@ Move around until a mate is found. self.mate = -1 return end - m, d = monthday(model.date) - nest = ((m == self.nestingbegin[1] && d >= self.nestingbegin[2] - && @chance(0.05)) || (m > self.nestingbegin[1])) #XXX why the chance? + month, day = monthday(model.date) + nest = ((month == self.nestingbegin[1] && day >= self.nestingbegin[2] + && @chance(0.05)) || (month > self.nestingbegin[1])) #XXX why the chance? nest && @setphase(nestbuilding) return end # look for a mate among the neighbouring birds, or move randomly - for n in @neighbours(50) #XXX magic number + for n in @neighbours(500m) #XXX magic number if n.sex != self.sex && n.phase == mating && n.mate == -1 self.mate = n.id n.mate = self.id @@ -165,7 +163,7 @@ Move around until a mate is found. end end #@debug("$(animalid(self)) didn't find a mate.") - @walk("random", 10) #XXX magic number + @walk("random", 100m) #XXX magic number end """ @@ -181,9 +179,9 @@ Females select a location and build a nest. Males do nothing. (Sound familiar?) if self.sex == female if isempty(self.nest) # try to find a nest in the neighbourhood, or move on - nestlocation = @randompixel(10, self.habitats) #XXX magic number + nestlocation = @randompixel(100m, self.habitats) #XXX magic number if isnothing(nestlocation) - @walk("random", 20) #XXX magic number + @walk("random", 200m) #XXX magic number else # if we've found a location, start the clock on the building time # (building time doubles for the first nest of the year) @@ -205,7 +203,7 @@ Females select a location and build a nest. Males do nothing. (Sound familiar?) else # males stay near the female mate = @animal(self.mate) - @follow(mate, 5) + @follow(mate, 50m) mate.phase == breeding && @setphase(breeding) end end @@ -256,9 +254,9 @@ should currently be on migration. Also sets other individual-specific variables. # calculate migration dates for this individual self.migrationdates = migrationdates(self, model) leave, arrive = self.migrationdates - m, d = monthday(model.date) - migrate = (((m < arrive[1]) || (m == arrive[1] && d < arrive[2])) || - ((m > leave[1]) || (m == leave[1] && d >= leave[2]))) + month, day = monthday(model.date) + migrate = (((month < arrive[1]) || (month == arrive[1] && day < arrive[2])) || + ((month > leave[1]) || (month == leave[1] && day >= leave[2]))) if migrate returndate = Date(year(model.date), arrive[1], arrive[2]) model.date != @param(core.startdate) && (returndate += Year(1)) @@ -271,6 +269,6 @@ end habitat = skylarkhabitat initphase = mating birthphase = egg - popdensity = 300 #XXX use Unitful.jl + indarea = 3ha pairs = true end diff --git a/src/nature/species/wolpertinger.jl b/src/nature/species/wolpertinger.jl index 83177fcebb049657ca6a87a9d30c7ab707b7737a..31fa70281328b14d7dfa7123737204055364db7c 100644 --- a/src/nature/species/wolpertinger.jl +++ b/src/nature/species/wolpertinger.jl @@ -27,7 +27,7 @@ and occasionally reproduce by spontaneous parthenogenesis... #walk!(animal, direction, model; ifempty=false) end - if @rand() < self.fecundity && length(@neighbours(10)) < self.crowding + if @rand() < self.fecundity && length(@neighbours(100m)) < self.crowding @reproduce() end @@ -49,5 +49,5 @@ Population densities of the endangered Wolpertinger are down to 1 animal per km asexual = true initphase = lifephase birthphase = lifephase - popdensity = 10000 #XXX use Unitful.jl for conversion? + indarea = 1km² end diff --git a/src/nature/species/wyvern.jl b/src/nature/species/wyvern.jl index d12c3952f335dca726bfb20e0228853d8633013b..da00307520e1f91dfc9f365b4f2e668bff31d2ce 100644 --- a/src/nature/species/wyvern.jl +++ b/src/nature/species/wyvern.jl @@ -12,8 +12,8 @@ legs, but that doesn't make it any less dangerous... @species Wyvern begin fecundity = 0.02 maxage = 365 - speed = 20 - vision = 50 + speed = 20m + vision = 150m aggression = 0.2 huntsuccess = 0.8 end diff --git a/src/parameters.toml b/src/parameters.toml index f9da32f0abdeda332bb44caadb8abfc2126d9409..2f2215b842291d7d0a5b65f6c7dc223cc7ee3299 100644 --- a/src/parameters.toml +++ b/src/parameters.toml @@ -10,7 +10,7 @@ configfile = "src/parameters.toml" # location of the configuration file outdir = "results" # location and name of the output folder overwrite = "ask" # overwrite the output directory? (true/false/"ask") -logoutput = "both" # log output to screen/file/none/both #XXX default "both" +logoutput = "both" # log output to screen/file/none/both csvoutput = true # save collected data in CSV files visualise = true # generate result graphs storedata = true # keep collected data in memory @@ -24,6 +24,7 @@ enddate = 2022-12-31 # last day of the simulation [world] mapdirectory = "data/regions/jena-small" # the directory in which all geographic data are stored +mapresolution = 10 # map resolution in meters landcovermap = "landcover.tif" # name of the landcover map in the map directory farmfieldsmap = "fields.tif" # name of the field geometry map in the map directory weatherfile = "weather.csv" # name of the weather data file in the map directory diff --git a/src/world/landscape.jl b/src/world/landscape.jl index 36a47febe98f6f201b09434e2f6cd62a37bcbf3d..d78bd7084fc6479a7064c7035ac9b401c7f0dbd6 100644 --- a/src/world/landscape.jl +++ b/src/world/landscape.jl @@ -180,9 +180,10 @@ Calculate the distance from the given location to the closest location matching 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)) + return maximum(abs.(target))*@param(world.mapresolution) end """ @@ -213,8 +214,9 @@ end 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::Int64=1, +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 @@ -225,11 +227,12 @@ function randompixel(pos::Tuple{Int64,Int64}, model::SimulationModel, range::Int end """ - randomdirection(model, range=1) + randomdirection(model, distance) -Get a random direction coordinate tuple within the specified range. +Get a random direction coordinate tuple within the specified distance. """ -function randomdirection(model::SimulationModel, range::Int64=1) +function randomdirection(model::SimulationModel, distance::Length) + range = Int(floor(distance / @param(world.mapresolution))) Tuple(@rand(-range:range, 2)) end