From 8e3bed2df1c9c10bdce0d389f62e9f212fe3326b Mon Sep 17 00:00:00 2001 From: Daniel Vedder <daniel.vedder@idiv.de> Date: Wed, 7 Aug 2024 14:37:01 +0200 Subject: [PATCH] Created skylark stats plot --- src/Persefone.jl | 9 ++++++++- src/analysis/makieplots.jl | 37 ++++++++++++++++++++++++++++++++--- src/crop/farmplot.jl | 2 +- src/nature/ecologicaldata.jl | 32 +++--------------------------- src/nature/individuals.jl | 2 +- src/nature/macros.jl | 10 ++++++++++ src/nature/populations.jl | 13 ++++++++++++ src/nature/species/skylark.jl | 5 ++++- src/parameters.toml | 2 +- 9 files changed, 75 insertions(+), 37 deletions(-) diff --git a/src/Persefone.jl b/src/Persefone.jl index 3ff4da7..37ed428 100644 --- a/src/Persefone.jl +++ b/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 diff --git a/src/analysis/makieplots.jl b/src/analysis/makieplots.jl index cc10094..a0b5d37 100644 --- a/src/analysis/makieplots.jl +++ b/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 diff --git a/src/crop/farmplot.jl b/src/crop/farmplot.jl index c03d2c2..ebe7c84 100644 --- a/src/crop/farmplot.jl +++ b/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 """ diff --git a/src/nature/ecologicaldata.jl b/src/nature/ecologicaldata.jl index 8f0b7e8..3a1e538 100644 --- a/src/nature/ecologicaldata.jl +++ b/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 diff --git a/src/nature/individuals.jl b/src/nature/individuals.jl index 7cb9da5..126aa6a 100644 --- a/src/nature/individuals.jl +++ b/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) diff --git a/src/nature/macros.jl b/src/nature/macros.jl index 250d688..53da2a6 100644 --- a/src/nature/macros.jl +++ b/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) diff --git a/src/nature/populations.jl b/src/nature/populations.jl index a31a840..2f45c71 100644 --- a/src/nature/populations.jl +++ b/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) diff --git a/src/nature/species/skylark.jl b/src/nature/species/skylark.jl index 547f1cd..ac582dc 100644 --- a/src/nature/species/skylark.jl +++ b/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 diff --git a/src/parameters.toml b/src/parameters.toml index 659e887..ee391ab 100644 --- a/src/parameters.toml +++ b/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 -- GitLab