diff --git a/src/Persefone.jl b/src/Persefone.jl
index 65c3b2937d39821ba84bfc42cd774349daa07474..509f728a8ca4304df5575f171becfa0dff3e77ee 100644
--- a/src/Persefone.jl
+++ b/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
diff --git a/src/core/simulation.jl b/src/core/simulation.jl
index 18c0dd03d55048fe21b1714185ce937be04c1af3..96a9f92d60af2e1b50259ecc3b3087d73772d398 100644
--- a/src/core/simulation.jl
+++ b/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)
diff --git a/src/nature/insects.jl b/src/nature/insects.jl
new file mode 100644
index 0000000000000000000000000000000000000000..783849279a5ff398cf1330be227e34a0bf477841
--- /dev/null
+++ b/src/nature/insects.jl
@@ -0,0 +1,75 @@
+### 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
diff --git a/src/parameters.toml b/src/parameters.toml
index 4e888a4a9fa28abb1999209f3d4f42c709a8fb51..0838a3fe94aa857d6ce074bc48ba1a8b6ca821e9 100644
--- a/src/parameters.toml
+++ b/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)
diff --git a/src/world/landscape.jl b/src/world/landscape.jl
index d6b16f608a5693dcd48027065b7d9658736f45f7..e3dcd150ed4173a2be266d133d6c806e1e964eae 100644
--- a/src/world/landscape.jl
+++ b/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
diff --git a/test/nature_tests.jl b/test/nature_tests.jl
index b357d7b53a775a21238666cc47c55a152788cad4..2f483f9cd098cba54a8f9ef4e25b4c484a0e5edb 100644
--- a/test/nature_tests.jl
+++ b/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)
diff --git a/test/runtests.jl b/test/runtests.jl
index a1de751da724b45a6d3678412639a97c65e48929..db3a518bea9151c39b8191e1910a1da527fd0ee3 100644
--- a/test/runtests.jl
+++ b/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,