Make AquaCrop usable with weather and soil data and crop types

data/crops/aquacrop/crop_data.csv

+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

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.

src/crop/aquacrop.jl

-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

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

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

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

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/"