Added units with Unitful.jl

closes #85

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
 

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>`."

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)

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}}

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

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?)

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

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

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

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

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