From 290411d83801f8ca89d2120f41bbfe7b6116e906 Mon Sep 17 00:00:00 2001
From: Marco Matthies <71844+marcom@users.noreply.github.com>
Date: Mon, 24 Feb 2025 21:39:07 +0100
Subject: [PATCH] Make AquaCrop usable with weather and soil data and crop
types
---
data/crops/aquacrop/crop_data.csv | 9 ++
src/crop/almass.jl | 2 +-
src/crop/aquacrop.jl | 161 ++++++++++++++++++++++++-----
src/crop/cropmodels.jl | 37 ++++---
src/parameters.toml | 6 +-
test/crop_tests.jl | 13 ++-
test/test_parameters_aquacrop.toml | 3 +-
7 files changed, 181 insertions(+), 50 deletions(-)
create mode 100644 data/crops/aquacrop/crop_data.csv
diff --git a/data/crops/aquacrop/crop_data.csv b/data/crops/aquacrop/crop_data.csv
new file mode 100644
index 0000000..c9a871f
--- /dev/null
+++ b/data/crops/aquacrop/crop_data.csv
@@ -0,0 +1,9 @@
+persefone_name,aquacrop_name,crop_group,min_sow_date,max_sow_date
+natural grass,alfalfaGDD,semi-natural,NA,NA
+permanent grassland (grazed),alfalfaGDD,grass,NA,NA
+permanent set-aside,alfalfaGDD,grass,NA,NA
+no growth,alfalfaGDD,semi-natural,NA,NA
+winter rape,rapeseed,grain,20 August,25 August
+winter wheat,wheatGDD,grain,15 October,31 October
+maize,maize,grain,15 April,30 April
+winter barley,barleyGDD,grain,15 September,30 September
\ No newline at end of file
diff --git a/src/crop/almass.jl b/src/crop/almass.jl
index 59b94b9..8edc318 100644
--- a/src/crop/almass.jl
+++ b/src/crop/almass.jl
@@ -398,7 +398,7 @@ Parse a CSV file containing the required parameter values for each crop
(as produced from the original ALMaSS file by `convert_almass_data.py`).
"""
function readcropparameters(cropdirectory::String)
- @debug "Reading crop parameters"
+ @debug "Reading crop parameters for ALMaSS crop model"
# TODO: the second-last column (sowingdensity) uses `String` as
# type. This is because the entry for "peas/beans" has a value of
# "35/80" for the sowingdensity.
diff --git a/src/crop/aquacrop.jl b/src/crop/aquacrop.jl
index 7e09b2c..205078d 100644
--- a/src/crop/aquacrop.jl
+++ b/src/crop/aquacrop.jl
@@ -1,13 +1,19 @@
-module WrapAquaCrop
-
# Persefone.jl wrapper for AquaCrop.jl
-# Use different name from `AquaCrop` to avoid name clash
+
+module AquaCropWrapper
+
+const CROPFILE = "crop_data.csv"
import AquaCrop
+import CSV
+using Dates: Date
+using DataFrames: DataFrame
using Persefone:
AnnualDate,
cm,
+ daynumber,
+ SoilType,
SimulationModel
import Persefone:
@@ -26,48 +32,146 @@ using Unitful: @u_str
# We can't use Length{Float64} as that is a Union type
const Tlength = typeof(1.0cm)
-struct AquaCropType
+# TODO: read crop names directly from AquaCrop.jl
+# names extracted by running this command in the AquaCrop.jl src dir:
+# grep -nir 'crop_type ==' src/ | cut -d '=' -f 3 | sed 's/$/,/' | sort | uniq
+const AQUACROP_CROPNAMES = [
+ "alfalfaGDD",
+ "barley",
+ "barleyGDD",
+ "cotton",
+ "cottonGDD",
+ "drybean",
+ "drybeanGDD",
+ "maize",
+ "maizeGDD",
+ "oat",
+ "paddyrice",
+ "paddyriceGDD",
+ "potato",
+ "potatoGDD",
+ "quinoa",
+ "rapeseed",
+ "sorghum",
+ "sorghumGDD",
+ "soybean",
+ "soybeanGDD",
+ "sugarbeet",
+ "sugarbeetGDD",
+ "sugarcane",
+ "sunflower",
+ "sunflowerGDD",
+ "tef",
+ "tomato",
+ "tomatoGDD",
+ "wheat",
+ "wheatGDD",
+]
+
+@kwdef struct CropType
name::String
+ aquacrop_name::String
group::String
minsowdate::Union{Missing,AnnualDate}
maxsowdate::Union{Missing,AnnualDate}
end
-cropname(ct::AquaCropType) = ct.name
+cropname(ct::CropType) = ct.name
-mutable struct AquaCropState
- croptype::AquaCropType
+"""
+ readcropparameters(cropdirectory)
+
+Parse a CSV file containing some parameters required to map AquaCrop
+crop names to Persefone crop names as well as additional crop data
+needed for Persefone (cropgroup, minsowdate, maxsowdate).
+"""
+function readcropparameters(cropdirectory::String)
+ @debug "Reading extra crop parameters for AquaCrop crop model"
+ croptypes = Dict{String,CropType}()
+ df = CSV.read(joinpath(cropdirectory, CROPFILE), DataFrame;
+ missingstring="NA", types=[String, String, String, AnnualDate, AnnualDate])
+ for row in eachrow(df)
+ croptypes[row.persefone_name] =
+ CropType(name = row.persefone_name,
+ aquacrop_name = row.aquacrop_name,
+ group = row.crop_group,
+ minsowdate = row.min_sow_date,
+ maxsowdate = row.max_sow_date)
+ end
+ return croptypes
+end
+
+mutable struct CropState
+ croptype::CropType
height::Tlength # TODO: remove height field, supply from cropstate
+ soiltype::SoilType
cropstate::AquaCrop.AquaCropField
- function AquaCropState(croptype::AquaCropType, height::Tlength=0.0cm)
- runtype = AquaCrop.TomlFileRun()
- parentdir = AquaCrop.test_toml_dir # TODO: hardcoded croptype
- ac_cropfield, all_ok = AquaCrop.start_cropfield(; parentdir, runtype)
- if ! all_ok.logi
- error("AquaCrop.start_cropfield() failed, status = $all_ok")
- end
- AquaCrop.setup_cropfield!(ac_cropfield, all_ok; runtype=runtype)
+ function CropState(crop_type::CropType, soiltype::SoilType,
+ model::SimulationModel, height::Tlength=0.0cm)
+ # TODO: get this mapping for soil types from a CSV file?
+ soiltype_to_aquacrop = Dict(soiltype => replace(string(soiltype), "_" => " ")
+ for soiltype in instances(SoilType))
+ aquacrop_soiltype = soiltype_to_aquacrop[soiltype]
+
+ start_date = model.date # TODO: maybe `today(model)`
+ end_date = model.weather.lastdate # TODO: maybe `lastdate(model)`
+ start_daynum = daynumber(model.weather, start_date)
+ end_daynum = daynumber(model.weather, end_date)
+ dayrange = start_daynum:end_daynum
+ input_Tmin = @view model.weather.mintemp[dayrange]
+ input_Tmax = @view model.weather.maxtemp[dayrange]
+ input_ETo = @view model.weather.evapotranspiration[dayrange]
+ input_Rain = @view model.weather.precipitation[dayrange]
+
+ # Generate the keyword object for the AquaCrop simulation
+ kwargs = (
+ ## Necessary keywords
+ # runtype
+ runtype = AquaCrop.NoFileRun(),
+ # project input
+ Simulation_DayNr1 = start_date,
+ Simulation_DayNrN = end_date,
+ Crop_Day1 = start_date, # TODO: originally start_date + Week(1)
+ Crop_DayN = end_date,
+ # soil
+ soil_type = aquacrop_soiltype,
+ # crop
+ crop_type = crop_type.aquacrop_name,
+ # Climate
+ InitialClimDate = start_date,
+ ## Optional keyworkds
+ # Climate
+ Tmin = input_Tmin,
+ Tmax = input_Tmax,
+ ETo = input_ETo,
+ Rain = input_Rain,
+ # change soil properties
+ soil_layers = Dict("Thickness" => 5.0)
+ )
+
+ aquacrop_cropfield, all_ok = AquaCrop.start_cropfield(; kwargs...)
if ! all_ok.logi
- error("AquaCrop.setup_cropfield!() failed, status = $all_ok")
+ @error "Failed calling AquaCrop.start_cropfield()\nAquaCrop error: $(all_ok.msg)"
end
- return new(croptype, height, ac_cropfield)
+
+ return new(crop_type, height, soiltype, aquacrop_cropfield)
end
end
-croptype(cs::AquaCropState) = cs.croptype
-cropname(cs::AquaCropState) = cropname(croptype(cs))
-cropheight(cs::AquaCropState) = cs.height # TODO: use AquaCrop state info
-cropcover(cs::AquaCropState) = AquaCrop.canopycover(cs.cropstate)
-cropyield(cs::AquaCropState) = AquaCrop.dryyield(cs.cropstate) # TODO: there is also freshyield
-isharvestable(cs::AquaCropState) = true # TODO: implement this correctly
+croptype(cs::CropState) = cs.croptype
+cropname(cs::CropState) = cropname(croptype(cs))
+cropheight(cs::CropState) = cs.height # TODO: calculate from AquaCrop state info
+cropcover(cs::CropState) = AquaCrop.canopycover(cs.cropstate)
+cropyield(cs::CropState) = AquaCrop.dryyield(cs.cropstate) # TODO: there is also freshyield
+isharvestable(cs::CropState) = true # TODO: implement this correctly
"""
stepagent!(cropstate, model)
Update a crop state by one day.
"""
-function stepagent!(cs::AquaCropState, model::SimulationModel)
+function stepagent!(cs::CropState, model::SimulationModel)
AquaCrop.dailyupdate!(cs.cropstate)
end
@@ -76,7 +180,10 @@ end
Change the cropstate to sow the specified crop.
"""
-function sow!(cs::AquaCropState, model::SimulationModel, cropname::String)
+function sow!(cs::CropState, model::SimulationModel, cropname::String)
+ if cropname ∉ keys(model.crops)
+ @error "cropname \"$cropname\" not found"
+ end
cs.croptype = model.crops[cropname]
cs.height = 0.0cm
# TODO: other things needed for AquaCrop?
@@ -87,9 +194,9 @@ end
Harvest the crop of this cropstate.
"""
-function harvest!(cs::AquaCropState, model::SimulationModel)
+function harvest!(cs::CropState, model::SimulationModel)
# TODO: implement this
return 0.0u"g/m^2"
end
-end # module WrapAquaCrop
+end # module
diff --git a/src/crop/cropmodels.jl b/src/crop/cropmodels.jl
index 4a67eba..715e74b 100644
--- a/src/crop/cropmodels.jl
+++ b/src/crop/cropmodels.jl
@@ -23,12 +23,9 @@ function initcropmodel(cropmodel::AbstractString, cropdirectory::AbstractString)
crops = Dict(name => SimpleCrop.CropType(ct.name, ct.group, ct.minsowdate, ct.maxsowdate)
for (name, ct) in crops_almass)
elseif cropmodel == "aquacrop"
- Tcroptype = WrapAquaCrop.AquaCropType
- Tcropstate = WrapAquaCrop.AquaCropState
- # TODO: temporarily using ALMaSS crop types
- crops_almass = ALMaSS.readcropparameters(cropdirectory)
- crops = Dict(name => WrapAquaCrop.AquaCropType(ct.name, ct.group, ct.minsowdate, ct.maxsowdate)
- for (name, ct) in crops_almass)
+ Tcroptype = AquaCropWrapper.CropType
+ Tcropstate = AquaCropWrapper.CropState
+ crops = AquaCropWrapper.readcropparameters(cropdirectory)
else
error("initcropmodel: no implementation for crop model '$cropmodel'")
end
@@ -54,8 +51,8 @@ function initfields!(model::SimulationModel)
model.landscape[x,y].fieldid = objectid
push!(model.farmplots[objectid].pixels, (x,y))
else
- cropstate = makecropstate(model)
soiltype = model.landscape[x,y].soiltype
+ cropstate = makecropstate(model, soiltype)
fp = FarmPlot(length(model.farmplots) + 1, [(x, y)], -1, soiltype, cropstate)
push!(model.farmplots, fp)
model.landscape[x,y].fieldid = fp.id
@@ -64,20 +61,35 @@ function initfields!(model::SimulationModel)
end
end
- # Adjust soil type to most common soil type of the pixels
- # belonging to a farmplot
+ # Adjust farmplot soil type to most common soil type of its
+ # pixels, and set the cropstate soil type to this as well.
for fp in model.farmplots
fp.soiltype = mode(map(p -> model.landscape[p[1], p[2]].soiltype, fp.pixels))
+ set_cropstate_soiltype!(fp, model)
end
+
@info "Initialised $(length(model.farmplots)) farm plots."
end
+# TODO: this function should be moved to the individual crop models,
+# and overloaded on the type of cropstate
+"""
+ set_cropstate_soiltype!(farmplot, model)
+"""
+function set_cropstate_soiltype!(farmplot::FarmPlot, model::SimulationModel)
+ if @param(crop.cropmodel) == "aquacrop"
+ farmplot.cropstate.soiltype = farmplot.soiltype
+ else
+ # do nothing for other cropmodels
+ end
+end
+
"""
- makecropstate(model, cropmodel)
+ makecropstate(model, soiltype)
An internal utility function to initialise one instance of the configured crop growth model.
"""
-function makecropstate(model::SimulationModel)
+function makecropstate(model::SimulationModel, soiltype::SoilType)
if @param(crop.cropmodel) == "almass"
phase = (month(model.date) < 3 ? ALMaSS.janfirst : ALMaSS.marchfirst)
cs = ALMaSS.CropState(
@@ -90,7 +102,8 @@ function makecropstate(model::SimulationModel)
0.0m
)
elseif @param(crop.cropmodel) == "aquacrop"
- cs = WrapAquaCrop.AquaCropState(model.crops["natural grass"], 0.0cm)
+ croptype = model.crops["natural grass"]
+ cs = AquaCropWrapper.CropState(croptype, soiltype, model)
else
Base.error("Unhandled crop model '$(@param(crop.cropmodel))' in makecropstate().")
end
diff --git a/src/parameters.toml b/src/parameters.toml
index 743343a..c3d8aff 100644
--- a/src/parameters.toml
+++ b/src/parameters.toml
@@ -23,7 +23,7 @@ enddate = 2023-12-31 # last day of the simulation
#enddate = 2022-01-02 # last day of the simulation (test value)
[world]
-mapdirectory = "data/regions/jena" # the directory in which all geographic data are stored
+mapdirectory = "data/regions/jena/" # the directory in which all geographic data are stored
mapresolution = 10 # map resolution in meters
landcovermap = "landcover.tif" # name of the landcover map in the map directory
farmfieldsmap = "fields.tif" # name of the field geometry map in the map directory
@@ -46,5 +46,5 @@ indoutfreq = "end" # output frequency individual-level data, daily/monthly/yearl
insectmodel = ["season", "habitat", "pesticides", "weather"] # factors affecting insect growth
[crop]
-cropmodel = "almass" # crop growth model to use, "almass", "aquacrop", or "simple"
-cropdirectory = "data/crops/almass" # the directory storing all data files for the selected crop model
+cropmodel = "aquacrop" # crop growth model to use, "almass", "aquacrop", or "simple"
+cropdirectory = "data/crops/aquacrop/" # the directory storing all data files for the selected crop model
diff --git a/test/crop_tests.jl b/test/crop_tests.jl
index 17ec180..e80ed59 100644
--- a/test/crop_tests.jl
+++ b/test/crop_tests.jl
@@ -57,15 +57,18 @@ end
end
@testset "Submodule AquaCrop" begin
- PsAC = Ps.WrapAquaCrop
- ct = PsAC.AquaCropType("olive tree", "no_crop_group", AnnualDate(4, 1), AnnualDate(5, 1))
+ model = initialise(; configfile=TESTPARAM_AQUACROP)
+
+ PsAC = Ps.AquaCropWrapper
+ ct = PsAC.CropType("olive tree", "alfalfaGDD", "no_crop_group", AnnualDate(4, 1), AnnualDate(5, 1))
id = 0
pixels = [(0, 0)]
farmer = 0
- fp = FarmPlot(id, pixels, farmer, PsAC.AquaCropState(ct))
+ soiltype = Ps.sand
+ fp = FarmPlot(id, pixels, farmer, soiltype, PsAC.CropState(ct, soiltype, model))
@test fp isa FarmPlot
- @test fp isa FarmPlot{PsAC.AquaCropState}
- @test croptype(fp) isa PsAC.AquaCropType
+ @test fp isa FarmPlot{PsAC.CropState}
+ @test croptype(fp) isa PsAC.CropType
@test cropname(fp) isa String
@test cropheight(fp) isa Length{Float64}
@test cropcover(fp) isa Float64
diff --git a/test/test_parameters_aquacrop.toml b/test/test_parameters_aquacrop.toml
index 82bc805..9f29fdc 100644
--- a/test/test_parameters_aquacrop.toml
+++ b/test/test_parameters_aquacrop.toml
@@ -42,5 +42,4 @@ insectmodel = ["season", "habitat", "pesticides"] # which factors affect insect
[crop]
cropmodel = "aquacrop" # crop growth model to use, "almass", "aquacrop", or "simple"
-# TODO: dir is for ALMaSS parameters
-cropdirectory = "data/crops/almass/"
+cropdirectory = "data/crops/aquacrop/"
--
GitLab