Started working on the insect biomass submodel.

src/Persefone.jl

@@ -31,6 +31,7 @@ using
 # 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/
 
 ## define exported functions and variables
 export
@@ -82,6 +83,7 @@ include("farm/farm.jl")
 include("crop/crops.jl")
 
 include("nature/nature.jl")
+include("nature/insects.jl")
 include("nature/populations.jl")
 include("nature/ecologicaldata.jl")
 #TODO loop over files in nature/species directory

src/core/simulation.jl

@@ -118,6 +118,7 @@ Execute one update of the model.
 function stepsimulation!(model::AgentBasedModel)
     with_logger(model.logger) do
         @info "Simulating day $(model.date)."
+        updateinsectbiomass!(model)
         for a in Schedulers.ByType((Farmer,FarmPlot,Animal), true)(model)
             try #The animal may have been killed
                 stepagent!(model[a], model)

src/nature/insects.jl

+### Persefone - a socio-economic-ecological model of European agricultural landscapes.
+###
+### This file contains the functions for calculating insect biomass.
+###
+
+##FIXME Updating every pixel every day is really expensive! I should seriously
+## think about switching this to a "calculate-on-demand" design.
+## (The problem is that that works for all models except the pesticide one, which
+## is a Markov Chain model that cannot be calculated without knowing the previous
+## value at the location.)
+
+##TODO write tests for this
+
+"""
+    updateinsectbiomass!(model)
+
+Calculate a new insect biomass value for each pixel in the landscape.
+"""
+function updateinsectbiomass!(model::AgentBasedModel)
+    @debug "Updating insect biomass"
+    if @param(nature.insectmodel) == "seasonal"
+        updatefunction = seasonalmodel
+    elseif @param(nature.insectmodel) == "seasonalhabitat"
+        updatefunction = seasonalhabitatmodel
+    else
+        Base.error("Unknown insect model $(@param(nature.insectmodel)).")
+        #TODO replace with exception
+    end
+    for p in model.landscape
+        p.insectbiomass = updatefunction(p, model)
+    end
+end
+
+
+"""
+    seasonalmodel(pixel)
+
+Calculate an insect biomass value based purely on the calendar day, following
+a parabolic function as per Paquette et al. (2013). Assumes a maximum insect
+density of 300 mg/m² (cf. Odderskaer et al. 1997, Püttmanns et al. 2022).
+
+*[Note that this is a very approximate calculation!]*
+"""
+function seasonalmodel(pixel::Pixel, model::AgentBasedModel)
+    calendarday = dayofyear(model.date)
+    biomass = -0.085(calendarday-187)^2+300
+    biomass > 0 ? biomass : 0
+end
+
+"""
+    seasonalhabitatmodel(pixel)
+
+Calculate an insect biomass value based on the calendar day and habitat type,
+following a parabolic function.
+
+See Odderskaer et al. (1997), Grübler et al. (2008), Paquette et al. (2013),
+and Püttmanns et al. (2022).
+
+*[Note that this is a very approximate calculation!]*
+"""
+function seasonalhabitatmodel(pixel::Pixel, model::AgentBasedModel)
+    calendarday = dayofyear(model.date)
+    habitatfactor = 0.0
+    if pixel.landcover == soil || pixel.landcover == forest
+        habitatfactor = 2.0
+    elseif pixel.landcover == grass
+        habitatfactor = 1.5
+    elseif pixel.landcover == agriculture
+        habitatfactor = 1.0
+    end
+    biomass = -0.085(calendarday-187)^2+300*habitatfactor
+    biomass > 0 ? biomass : 0
+end
+
+#TODO add models that include weather and pesticide application

src/parameters.toml

@@ -29,6 +29,7 @@ farmmodel = "FieldManager" # which version of the farm model to use (not yet imp
 targetspecies = ["Wolpertinger", "Wyvern"] # list of target species to simulate
 popoutfreq = "daily" # output frequency population-level data, daily/monthly/yearly/end/never
 indoutfreq = "end" # output frequency individual-level data, daily/monthly/yearly/end/never
+insectmodel = "seasonalhabitat" # how to calculate insect biomass, seasonal/seasonalhabitat/full
 	
 [crop]
 cropmodel = "linear" # crop growth model to use, "linear" or "aquacrop" (not yet implemented)

src/world/landscape.jl

@@ -21,6 +21,7 @@ mutable struct Pixel
     landcover::LandCover
     fieldid::Union{Missing,Int64}
     events::Vector{EventType}
+    insectbiomass::Float64 # in mg/m²
 end
 
 """
@@ -57,7 +58,7 @@ function initlandscape(landcovermap::String, farmfieldsmap::String)
             lcv = LandCover(Int(lc/10))
             ff = farmfields[x,y][1]
             !(ismissing(ff)) && (ff = Int64(ff))
-            landscape[x,y] = Pixel(lcv, ff, Vector{Symbol}())
+            landscape[x,y] = Pixel(lcv, ff, Vector{Symbol}(), 0.0)
         end
     end
     return landscape

test/nature_tests.jl

@@ -6,7 +6,7 @@
 @testset "Habitat macros" begin
     # set up the testing landscape
     model = smalltestlandscape()
-    model.landscape[6,6] = Pixel(Ps.agriculture, 1, [])
+    model.landscape[6,6] = Pixel(Ps.agriculture, 1, [], 0.0)
     species::Dict{String,Any} = Dict("name"=>"test_animal")
     add_agent!((6,6), FarmPlot, model, [(6,6)], Ps.wheat, 1.2, 3.4)
     add_agent!((3,3), Animal, model, species, Ps.male, 1)

test/runtests.jl

@@ -37,10 +37,10 @@ function smalltestlandscape()
     for x in 1:6
         for y in 1:6
             (x in (1:2) || y in (1:4)) ? lc = Ps.forest : lc = Ps.grass
-            landscape[x,y] = Pixel(lc, missing, [])
+            landscape[x,y] = Pixel(lc, missing, [], 0.0)
         end
     end
-    landscape[6,4] = Pixel(Ps.water, 0, [])
+    landscape[6,4] = Pixel(Ps.water, 0, [], 0.0)
     space = GridSpace(size(landscape), periodic=false)
     properties = Dict{Symbol,Any}(:date=>TESTSETTINGS["core.startdate"],
                                   :landscape=>landscape,