Created skylark stats plot

src/Persefone.jl

@@ -56,6 +56,11 @@ export
     @param,
     @rand,
     @shuffle!,
+    @thisyear,
+    @nextyear,
+    @lastyear,
+    @record,
+    @data,
     @species,
     @populate,
     @create,
@@ -88,6 +93,7 @@ export
     @move,
     @walk,
     @follow,
+    @destroynest,
     #functions
     simulate,
     simulate!,
@@ -130,7 +136,6 @@ function stepagent! end
 include("core/utils.jl")
 include("core/input.jl")
 include("core/output.jl")
-include("analysis/makieplots.jl")
 
 include("world/landscape.jl")
 include("world/weather.jl")
@@ -153,6 +158,8 @@ include("nature/species/skylark.jl")
 include("nature/species/wolpertinger.jl")
 include("nature/species/wyvern.jl")
 
+include("analysis/makieplots.jl")
+
 include("core/simulation.jl") #this must be last
 
 # precompile important functions - TODO use PrecompileTools.jl

src/analysis/makieplots.jl

@@ -25,7 +25,7 @@ function visualisemap(model::SimulationModel,date=nothing,landcover=nothing)
     image!(f[1,1], landcover)
     ax.aspect = DataAspect()
     # check if there are individuals and plot them
-    inds = @subset(model.dataoutputs["individuals"].datastore, :X .== -1) #remove migrants
+    inds = @subset(@data("individuals"), :X .== -1) #remove migrants
     if iszero(size(inds)[1])
         @debug "No individual data to map"
         return f
@@ -54,7 +54,7 @@ Plot a line graph of population sizes of each species over time.
 Returns a Makie figure object.
 """
 function populationtrends(model::SimulationModel)
-    pops = model.dataoutputs["populations"].datastore
+    pops = @data("populations")
     ncolors = max(2, length(@param(nature.targetspecies)))
     update_theme!(palette=(color=cgrad(:seaborn_bright, ncolors),), cycle=[:color])
     f = Figure()
@@ -78,7 +78,7 @@ Plot a line graph of total population size and individual demographics of skylar
 Returns a Makie figure object.
 """
 function skylarkpopulation(model::SimulationModel)
-    pops = model.dataoutputs["skylark_abundance"].datastore
+    pops = @data("skylark_abundance")
     f = Figure()
     dates = @param(core.startdate):@param(core.enddate)
     axlimits = (1, length(dates), 0, maximum(pops[!,:TotalAbundance]))
@@ -94,6 +94,35 @@ function skylarkpopulation(model::SimulationModel)
     f
 end
 
+"""
+    skylarkstats(model)
+
+Plot various statistics from the skylark model: nesting habitat, territory size, mortality.
+"""
+function skylarkstats(model::SimulationModel)
+    f = Figure()
+    nestingdata = @data("skylark_breeding")
+    mortalitydata = @subset(@data("mortality"), :Species .== "Skylark")
+    # nesting habitat
+    habitattype = unique(nestingdata.Landcover)
+    htfrequency = [count(x -> x==h, nestingdata.Landcover) for h in habitattype]
+    barplot(f[1,1], 1:length(habitattype), htfrequency,
+             axis = (xticks=(1:length(habitattype), habitattype), title="Nesting habitat"))
+    croptype = unique(nestingdata.Crop)
+    ctfrequency = [count(x -> x==c, nestingdata.Crop) for c in croptype]
+    barplot(f[2,1], 1:length(croptype), ctfrequency,
+             axis = (xticks=(1:length(croptype), croptype), title="Nesting crop type"))
+    # mortality
+    causes = unique(mortalitydata.Cause)
+    mtfrequency = [count(x -> x==c, mortalitydata.Cause) for c in causes]
+    barplot(f[1,2], 1:length(causes), mtfrequency,
+            axis = (xticks=(1:length(causes), causes), title="Causes of Mortality"))
+    # territory size
+    hist(f[2,2], nestingdata.Territory, bins=20,
+         axis=(ylabel="Frequency", title="Territory size (ha)"))
+    f
+end
+
 """
     gettickmarks(dates)
 
@@ -116,3 +145,5 @@ function gettickmarks(dates)
 end
 
 #XXX add animation with record()? https://docs.makie.org/stable/explanations/animation/
+
+#TODO causes of mortality, nesting habitats, territory sizes

src/crop/farmplot.jl

@@ -79,7 +79,7 @@ Return the name of the crop at this position, or an empty string if there is no
 """
 function cropname(pos::Tuple{Int64,Int64}, model::SimulationModel)
     field = model.landscape[pos...].fieldid
-    ismissing(field) ? "" : cropname(model.farmplots[field])
+    ismissing(field) ? "NA" : cropname(model.farmplots[field])
 end
 
 """

src/nature/ecologicaldata.jl

@@ -72,10 +72,9 @@ function initskylarkdata(model::SimulationModel)
                    "daily", skylarkabundance, skylarkpopulation)
     # newdataoutput!(model, "skylark_territories", ["Date", "ID", "X", "Y"],
     #                skylarkterritories, "monthly") #TODO add plotting function
-    newdataoutput!(model, "skylark_nests", ["Date", "ID", "X", "Y", "Landcover", "Crop"],
-                   "monthly", skylarknests) #TODO add plotting function
-    # newdataoutput!(model, "skylark_mortality", ["Date", "N", "Cause"],
-    #                skylarkmortality, "daily") #TODO add plotting function
+    newdataoutput!(model, "skylark_breeding",
+                   ["Date", "Female", "Male", "X", "Y", "Landcover", "Crop", "Territory"],
+                   @param(nature.popoutfreq), nothing, skylarkstats)
 end
 
 """
@@ -130,28 +129,3 @@ function skylarkterritories(model::SimulationModel)
     end
     data
 end
-
-"""
-    skylarknests(model)
-
-Return a list of coordinates of active skylark nests, and the habitat type they occupy.
-"""
-function skylarknests(model::SimulationModel)
-    data = []
-    for a in model.animals
-        (isnothing(a) || a.sex != female || speciesof(a) != "Skylark" || isempty(a.nest)) && continue
-        push!(data, [model.date, a.id, a.nest[1], a.nest[2],
-                     landcover(a.nest, model), cropname(a.nest, model)])
-    end
-    data 
-end
-
-"""
-    skylarkmortality(model)
-
-
-"""
-function skylarkmortality(model::SimulationModel)
-    #TODO
-    #XXX may require an additional feature in output.jl to account for a different output type
-end

src/nature/individuals.jl

@@ -82,7 +82,7 @@ Add the given location to the animal's territory. Returns `true` if successful
 """
 function occupy!(animal::Animal, model::SimulationModel, position::Tuple{Int64,Int64})
     if isoccupied(model, speciesof(animal), position) #XXX should this be an error?
-        @warn "Position $position is already occupied by a $(speciesof(animal))." animal
+        @warn "Position $position is already occupied by a $(speciesof(animal))." # animal
         return false
     else
         push!(animal.territory, position)

src/nature/macros.jl

@@ -531,6 +531,16 @@ macro record(args...)
     :(record!($(esc(:model)), $(map(esc, args)...)))
 end
 
+"""
+    @data(outputname)
+
+Return the data stored in the given output (assumes `core.storedata` is true).
+Only use in scopes where `model` is available.
+"""
+macro data(outputname)
+    :(data($(esc(:model)).dataoutputs[$(esc(outputname))]))
+end
+
 """
     @destroynest(reason)
 

src/nature/populations.jl

@@ -151,6 +151,19 @@ function isoccupied(model::SimulationModel, species::String, position::Tuple{Int
     return false
 end
 
+"""
+    territorysize(animal, model, stripunits=false)
+
+Calculate the size of this animal's territory in the given unit. If `stripunits` is true,
+return the size as a plain number.
+"""
+function territorysize(a::Union{Animal,Int64}, model::SimulationModel,
+                       units::Unitful.Units=ha, stripunits::Bool=false)
+    a isa Int && (a = @animal(a))
+    size = length(a.territory) * @param(world.mapresolution)^2 |> Float64
+    stripunits ? ustrip(units, size) : size |> units
+end
+
 """
    isalive(id, model)
 

src/nature/species/skylark.jl

@@ -211,6 +211,9 @@ Females that have found a partner build a nest and lay eggs in a suitable locati
                 #XXX all skylarks laying on the same day lay the same number of eggs? RNG?!
                 self.timer = @rand(self.nestbuildingtime)
                 (model.date == self.firstnest) && (self.timer *= 2) # the first nest takes twice as long to build
+                @record("skylark_breeding",
+                        [model.date, self.id, self.mate, pos[1], pos[2], string(@landcover()),
+                         @cropname(), territorysize(self.mate, model, ha, true)])
                 break
             end
         end
@@ -373,7 +376,7 @@ function destroynest!(self::Skylark, model::SimulationModel, reason::String)
     self.clutch = 0
     @debug("$(animalid(self)) had her nest destroyed by $reason.")
 end
-
+    
 
 ## INITIALISATION
 

src/parameters.toml

@@ -15,7 +15,7 @@ csvoutput = true # save collected data in CSV files
 visualise = true # generate result graphs
 storedata = true # keep collected data in memory
 figureformat = "pdf" # file format to use for graphical output
-loglevel = "debug" # verbosity level: "debug", "info", "warn"
+loglevel = "info" # verbosity level: "debug", "info", "warn"
 processors = 2 # number of processors to use on parallel runs
 seed = 2 # seed value for the RNG (0 -> random value)
 startdate = 2022-01-01 # first day of the simulation