From 77d355310ebf34281bbe1006fadc491b78f44e73 Mon Sep 17 00:00:00 2001 From: Daniel Vedder <daniel.vedder@idiv.de> Date: Wed, 2 Aug 2023 18:23:54 +0200 Subject: [PATCH] Started working on crop growth model --- data/convert_almass_data.py | 19 +++-- src/Persefone.jl | 3 +- src/core/simulation.jl | 2 +- src/crop/crops.jl | 136 ++++++++++++++++-------------------- src/crop/farmplot.jl | 114 ++++++++++++++++++++++++++++++ src/farm/farm.jl | 4 +- 6 files changed, 190 insertions(+), 88 deletions(-) create mode 100644 src/crop/farmplot.jl diff --git a/data/convert_almass_data.py b/data/convert_almass_data.py index 271519b..f071051 100755 --- a/data/convert_almass_data.py +++ b/data/convert_almass_data.py @@ -1,13 +1,11 @@ #!/usr/bin/python3 ### -### ALMaSS uses a pretty crazy file format to specify its crop growth parameters. -### This script converts the data to a saner long-table CSV format. +### ALMaSS uses a pretty convoluted file format to specify its crop growth parameters. +### This script converts the data to a more easily understandable long-table CSV format. ### ### Daniel Vedder, 02/08/2023 ### -import time - input_file = "almass_crop_growth_curves.pre" output_file = "almass_crop_growth_curves.csv" @@ -82,13 +80,13 @@ def interpret_crop_number(i): elif i in [4, 104]: crop_type = "winter wheat" elif i in [5, 105]: crop_type = "winter rye" elif i in [6, 106]: crop_type = "oats" - #elif i in [7, 107]: crop_type = "triticale" #FIXME 107 not defined, but probable - elif i == 7: crop_type = "triticale" + elif i in [7, 107]: crop_type = "triticale" #XXX 107 not defined, but probable + #elif i == 7: crop_type = "triticale" elif i in [8, 108]: crop_type = "maize" elif i in [13, 113]: crop_type = "undersown spring barley" elif i in [21, 121]: crop_type = "spring rape" - elif i == 22: crop_type = "winter rape" - #elif i in [22, 122]: crop_type = "winter rape" #FIXME 122 not defined, but probable + #elif i == 22: crop_type = "winter rape" + elif i in [22, 122]: crop_type = "winter rape" #XXX 122 not defined, but probable elif i == 25: crop_type = "permanent grassland, low yield" elif i == 26: crop_type = "permanent grassland, grazed" elif i == 27: crop_type = "permanent grassland, seeded" @@ -97,8 +95,8 @@ def interpret_crop_number(i): elif i == 30: crop_type = "peas/beans" elif i in [41, 141]: crop_type = "carrots" elif i in [50, 150]: crop_type = "potatoes" - #elif i in [60, 160]: crop_type = "beet" #FIXME 160 not defined, but probable - elif i == 60: crop_type = "beet" + elif i in [60, 160]: crop_type = "beet" #XXX 160 not defined, but probable + #elif i == 60: crop_type = "beet" #elif i in [65, 165]: crop_type = "NA" #FIXME unknown numbers elif i == 70: crop_type = "lucerne" elif i == 80: crop_type = "tulips" @@ -108,7 +106,6 @@ def interpret_crop_number(i): elif i == 94: crop_type = "lawn" #elif i == 99: crop_type = "NA" #FIXME unknown number elif i == 112: crop_type = "heath" - #elif i == 121: crop_type = "NA" #FIXME unknown number else: crop_type = "NA" print("WARNING: unknown crop type "+str(i)) diff --git a/src/Persefone.jl b/src/Persefone.jl index 1164f18..a9912f6 100644 --- a/src/Persefone.jl +++ b/src/Persefone.jl @@ -80,8 +80,9 @@ include("core/output.jl") include("world/landscape.jl") include("world/weather.jl") -include("farm/farm.jl") include("crop/crops.jl") +include("crop/farmplot.jl") +include("farm/farm.jl") include("nature/insects.jl") include("nature/energy.jl") diff --git a/src/core/simulation.jl b/src/core/simulation.jl index cff0692..ae42934 100644 --- a/src/core/simulation.jl +++ b/src/core/simulation.jl @@ -80,8 +80,8 @@ function initmodel(settings::Dict{String, Any}) model = AgentBasedModel(Union{Farmer,Animal,FarmPlot}, space, properties=properties, rng=StableRNG(settings["core.seed"]), warn=false) saveinputfiles(model) - initfarms!(model) initfields!(model) + initfarms!(model) initnature!(model) model end diff --git a/src/crop/crops.jl b/src/crop/crops.jl index 301b4d7..036cb1f 100644 --- a/src/crop/crops.jl +++ b/src/crop/crops.jl @@ -3,97 +3,85 @@ ### This file is responsible for managing the crop growth modules. ### -"The crop types simulated by the model" -@enum CropType fallow wheat barley maize rapeseed +@enum GrowthPhase janfirst sow marchfirst harvest1 harvest2 vegphase -#XXX not sure whether it makes sense to have this as an agent type, -# or perhaps better a grid property? """ - FarmPlot + CropType -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. +The type struct for all crops. Currently follows the crop +growth model as implemented in ALMaSS. """ -@agent FarmPlot GridAgent{2} begin - pixels::Vector{Tuple{Int64, Int64}} - croptype::CropType - height::Float64 - biomass::Float64 - #TODO +@kwdef struct CropType + name::String + sowdate::Date + harvestdate::Date + mingrowthtemp::Float64 + inflectionpoints::Dict{GrowthPhase,Vector{Float64}} + slopes::Dict{GrowthPhase,Vector{Float64}} + #issowable::Union{Function,Bool}=true end -""" - stepagent!(farmplot, model) +let croptypes = Dict{String,CropType}() + global function registercrop(ct::CropType) + croptypes[ct.name] = ct + end -Update a farm plot by one day. -""" -function stepagent!(farmplot::FarmPlot, model::AgentBasedModel) - #TODO + global function getcrop(name::String) + croptypes[name] + end end """ - initfields!(model) + @dday(datestring) -Initialise the model with its farm plots. +A utility macro to make it quicker to create dates in the +format "21 August". """ -function initfields!(model::AgentBasedModel) - n = 0 - convertid = Dict{Int64,Int64}() - width, height = size(model.landscape) - for x in 1:width - for y in 1:height - # for each pixel, we need to extract the field ID given by the map input file, - # and convert it into the internal object ID used by Agents.jl, creating a new - # agent object if necessary - rawid = model.landscape[x,y].fieldid - (ismissing(rawid)) && continue - if rawid in keys(convertid) - objectid = convertid[rawid] - model.landscape[x,y].fieldid = objectid - push!(model[objectid].pixels, (x,y)) - else - fp = add_agent!((x,y), FarmPlot, model, [(x,y)], fallow, 0.0, 0.0) - model.landscape[x,y].fieldid = fp.id - convertid[rawid] = fp.id - n += 1 - end - end - end - @info "Initialised $n farm plots." +macro dday(date::String) + :(Date($(date), dateformat"d U")) end -""" - averagefieldsize(model) +let croplist = Dict{String,CropType}() + + global function newcrop(;kwargs...) + crop = CropType(;kwargs...) + croplist[crop.name] = crop + end -Calculate the average field size in hectares for the model landscape. -""" -function averagefieldsize(model::AgentBasedModel) - conversionfactor = 100 #our pixels are currently 10x10m, so 100 pixels per hectare - sizes::Vector{Float64} = [] - for a in allagents(model) - (typeof(a) != FarmPlot) && continue - push!(sizes, size(a.pixels)[1]/conversionfactor) + global function cropdefinition(name) + croplist[name] end - round(sum(sizes)/size(sizes)[1], digits=2) end -""" - croptype(model, position) +# +# XXX It would be nicer to be able to define crops with a macro +# (like @species), but since macros are so tricky to get right, +# I'll leave that for later. +# +# """ +# @crop +# +# A macro to create and register a new crop type. Pass a list +# of keyword arguments to instantiate each field of the +# [`CropType`](@ref) struct. +# """ +# macro crop(;kwargs...) # macros don't take kwargs +# quote +# registercrop(CropType($(esc(kwargs)))) +# end +# end +# +# @crop "winter wheat" begin +# minsowdate = @dday("1 November") +# maxsowdate = @dday("31 November") +# growthparameter = 1.0 +# issowable = true +# end -Return the crop at this position, or nothing if there is no crop here (utility wrapper). -""" -function croptype(pos::Tuple{Int64,Int64}, model::AgentBasedModel) - ismissing(model.landscape[pos...].fieldid) ? nothing : - model[model.landscape[pos...].fieldid].croptype -end +newcrop(name="winter wheat", + sowdate = @dday("1 November"), + harvestdate = @dday("1 July"), + mingrowthtemp = 5.0, + growthparameter = 1.0, + issowable = true) -""" - cropheight(model, position) - -Return the height of the crop at this position, or nothing if there is no crop here -(utility wrapper). -""" -function cropheight(pos::Tuple{Int64,Int64}, model::AgentBasedModel) - ismissing(model.landscape[pos...].fieldid) ? nothing : - model[model.landscape[pos...].fieldid].height -end diff --git a/src/crop/farmplot.jl b/src/crop/farmplot.jl new file mode 100644 index 0000000..a77d0f8 --- /dev/null +++ b/src/crop/farmplot.jl @@ -0,0 +1,114 @@ +### Persefone - a socio-economic-ecological model of European agricultural landscapes. +### +### This file contains code for the fields that farmers manage. +### + +#XXX not sure whether it makes sense to have this as an agent type, +# or perhaps better a grid property? +""" + FarmPlot + +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. +""" +@agent FarmPlot GridAgent{2} begin + pixels::Vector{Tuple{Int64, Int64}} + croptype::CropType + growingdegreedays::Float64 + height::Float64 + biomass::Float64 + events::Vector{EventType} + #TODO +end + +""" + stepagent!(farmplot, model) + +Update a farm plot by one day. +""" +function stepagent!(farmplot::FarmPlot, model::AgentBasedModel) + # update growing degree days + gdd = meantemp(model) - farmplot.croptype.mingrowthtemp + gdd > 0 && (farmplot.growingdegreedays += gdd) + # update crop growth + growcrop!(farmplot, model) + #TODO expand? +end + +""" + initfields!(model) + +Initialise the model with its farm plots. +""" +function initfields!(model::AgentBasedModel) + n = 0 + convertid = Dict{Int64,Int64}() + width, height = size(model.landscape) + for x in 1:width + for y in 1:height + # for each pixel, we need to extract the field ID given by the map input file, + # and convert it into the internal object ID used by Agents.jl, creating a new + # agent object if necessary + rawid = model.landscape[x,y].fieldid + (ismissing(rawid)) && continue + if rawid in keys(convertid) + objectid = convertid[rawid] + model.landscape[x,y].fieldid = objectid + push!(model[objectid].pixels, (x,y)) + else + fp = add_agent!((x,y), FarmPlot, model, [(x,y)], fallow, 0.0, 0.0) + model.landscape[x,y].fieldid = fp.id + convertid[rawid] = fp.id + n += 1 + end + end + end + @info "Initialised $n farm plots." +end + +""" + averagefieldsize(model) + +Calculate the average field size in hectares for the model landscape. +""" +function averagefieldsize(model::AgentBasedModel) + conversionfactor = 100 #our pixels are currently 10x10m, so 100 pixels per hectare + sizes::Vector{Float64} = [] + for a in allagents(model) + (typeof(a) != FarmPlot) && continue + push!(sizes, size(a.pixels)[1]/conversionfactor) + end + round(sum(sizes)/size(sizes)[1], digits=2) +end + +""" + croptype(model, position) + +Return the crop at this position, or nothing if there is no crop here (utility wrapper). +""" +function croptype(pos::Tuple{Int64,Int64}, model::AgentBasedModel) + ismissing(model.landscape[pos...].fieldid) ? nothing : + model[model.landscape[pos...].fieldid].croptype +end + +""" + cropheight(model, position) + +Return the height of the crop at this position, or nothing if there is no crop here +(utility wrapper). +""" +function cropheight(pos::Tuple{Int64,Int64}, model::AgentBasedModel) + ismissing(model.landscape[pos...].fieldid) ? nothing : + model[model.landscape[pos...].fieldid].height +end + + +""" + growcrop!(farmplot, model) + +Apply the relevant crop growth model to update the plant height on this farm plot. +""" +function growcrop!(farmplot::FarmPlot, model::AgentBasedModel) + crop = farmplot.croptype + #TODO +end diff --git a/src/farm/farm.jl b/src/farm/farm.jl index 5a1089a..c42b18d 100644 --- a/src/farm/farm.jl +++ b/src/farm/farm.jl @@ -9,7 +9,9 @@ This is the agent type for the farm ABM. (Not yet implemented.) """ @agent Farmer GridAgent{2} begin - #TODO + fields::Vector{FarmPlot} + croprotation::Vector{CropType} + #TODO add AES end """ -- GitLab