AquaCrop: initial implementation

Project.toml

@@ -4,6 +4,7 @@ authors = ["Daniel Vedder <daniel.vedder@idiv.de>"]
 version = "0.6.1"
 
 [deps]
+AquaCrop = "8f16cebd-c0b4-44a3-857f-c101f083256c"
 ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
@@ -14,6 +15,7 @@ Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 GeoArrays = "2fb1d81b-e6a0-5fc5-82e6-8e06903437ab"
+ImageMagick = "6218d12a-5da1-5696-b52f-db25d2ecc6d1"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 LoggingExtras = "e6f89c97-d47a-5376-807f-9c37f3926c36"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
@@ -28,4 +30,5 @@ TiffImages = "731e570b-9d59-4bfa-96dc-6df516fadf69"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
+ImageMagick = "1.2.1"
 julia = ">= 1.10"

src/Persefone.jl

@@ -147,6 +147,7 @@ include("crop/farmplot.jl")
 include("crop/cropmodels.jl")
 include("crop/almass.jl")
 include("crop/simplecrop.jl")
+include("crop/aquacrop.jl")
 
 include("farm/farm.jl")
 include("farm/farmdata.jl")

src/core/input.jl

@@ -20,7 +20,7 @@ const PARAMFILE = joinpath(pkgdir(Persefone), "src/parameters.toml")
 """
 The crop models that can be used in the simulation.
 """
-const AVAILABLE_CROPMODELS = ["almass", "simple"]
+const AVAILABLE_CROPMODELS = ["almass", "simple", "aquacrop"]
 
 """
     @param(domainparam)

src/crop/aquacrop.jl

+import AquaCrop
+
+# We can't use Length{Float64} as that is a Union type
+const Tlength = typeof(1.0cm)
+
+struct AquaCropType
+    name::String
+    group::String
+    minsowdate::Union{Missing,AnnualDate}
+    maxsowdate::Union{Missing,AnnualDate}
+end
+
+cropname(ct::AquaCropType) = ct.name
+
+mutable struct AquaCropState
+    croptype::AquaCropType
+    height::Length{Float64}  # TODO: remove height field, supply from cropstate
+    cropstate::AquaCrop.AquaCropField
+
+    function AquaCropState(croptype::AquaCropType, height::Length{Float64}=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)
+        if ! all_ok.logi
+            error("AquaCrop.setup_cropfield!() failed, status = $all_ok")
+        end
+        return new(croptype, height, ac_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
+
+"""
+    stepagent!(cropstate, model)
+
+Update a crop state by one day.
+"""
+function stepagent!(cs::AquaCropState, model::SimulationModel)
+    AquaCrop.dailyupdate!(cs.cropstate)
+end
+
+"""
+    sow!(cropstate, model, cropname)
+
+Change the cropstate to sow the specified crop.
+"""
+function sow!(cs::AquaCropState, model::SimulationModel, cropname::String)
+    cs.croptype = model.crops[cropname]
+    cs.height = 0.0cm
+    # TODO: other things needed for AquaCrop?
+end
+
+"""
+    harvest!(cropstate, model)
+
+Harvest the crop of this cropstate.
+"""
+function harvest!(cs::AquaCropState, model::SimulationModel)
+    # TODO: implement this
+    return 0.0u"g/m^2"
+end

src/crop/cropmodels.jl

@@ -20,6 +20,13 @@ function initcropmodel(cropmodel::AbstractString, cropdirectory::AbstractString)
         crops_almass = ALMaSS.readcropparameters(cropdirectory)
         crops = Dict(name => SimpleCrop.CropType(ct.name, ct.group, ct.minsowdate, ct.maxsowdate)
                      for (name, ct) in crops_almass)
+    elseif cropmodel == "aquacrop"
+        Tcroptype = AquaCropType
+        Tcropstate = AquaCropState
+        # TODO: temporarily using ALMaSS crop types
+        crops_almass = ALMaSS.readcropparameters(cropdirectory)
+        crops = Dict(name => AquaCropType(ct.name, ct.group, ct.minsowdate, ct.maxsowdate)
+                     for (name, ct) in crops_almass)
     else
         error("initcropmodel: no implementation for crop model '$cropmodel'")
     end
@@ -73,6 +80,8 @@ function makecropstate(model::SimulationModel)
             model.crops["natural grass"],
             0.0m
         )
+    elseif @param(crop.cropmodel) == "aquacrop"
+        cs = AquaCropState(model.crops["natural grass"], 0.0cm)
     else
         Base.error("Unhandled crop model '$(@param(crop.cropmodel))' in makecropstate().")
     end

test/crop_tests.jl

@@ -3,10 +3,13 @@
 ### These are the tests for the crop growth models.
 ###
 
+import Unitful
+
 const TESTPARAM_ALMASS = joinpath(pkgdir(Persefone), "test", "test_parameters_almass.toml")
 const TESTPARAM_SIMPLECROP = joinpath(pkgdir(Persefone), "test", "test_parameters_simplecrop.toml")
+const TESTPARAM_AQUACROP = joinpath(pkgdir(Persefone), "test", "test_parameters_aquacrop.toml")
 
-@testset for configfile in (TESTPARAM_ALMASS, TESTPARAM_SIMPLECROP)
+@testset for configfile in (TESTPARAM_ALMASS, TESTPARAM_SIMPLECROP, TESTPARAM_AQUACROP)
     @testset "Model initialisation" begin
         model = initialise(; configfile)
         @test model isa AgricultureModel
@@ -15,6 +18,7 @@ const TESTPARAM_SIMPLECROP = joinpath(pkgdir(Persefone), "test", "test_parameter
         model = initialise(; configfile)
         stepsimulation!(model)
         @test model isa AgricultureModel
+        # TODO: test stepagent!(), sow!(), harvest!()
     end
 end
 
@@ -51,3 +55,18 @@ end
     @test cropcover(fp) isa Float64
     @test cropyield(fp) isa Float64
 end
+
+@testset "Submodule AquaCrop" begin
+    ct = Ps.AquaCropType("olive tree", "no_crop_group", AnnualDate(4, 1), AnnualDate(5, 1))
+    id = 0
+    pixels = [(0, 0)]
+    farmer = 0
+    fp = FarmPlot(id, pixels, farmer, Ps.AquaCropState(ct))
+    @test fp isa FarmPlot
+    @test fp isa FarmPlot{Ps.AquaCropState}
+    @test croptype(fp) isa Ps.AquaCropType
+    @test cropname(fp) isa String
+    @test cropheight(fp) isa Length{Float64}
+    @test cropcover(fp) isa Float64
+    @test Unitful.dimension(cropyield(fp)) == Unitful.𝐌 * Unitful.𝐋^-2
+end

test/test_parameters_aquacrop.toml

+### Persefone.jl - a model of agricultural landscapes and ecosystems in Europe.
+###
+### This configuration file is used for the test suite.
+###
+
+#XXX remember that changes here may break tests!
+
+[core]
+configfile = "test_parameters_aquacrop.toml" # location of the configuration file
+outdir = "results_testsuite" # location and name of the output folder
+logoutput = "both"
+overwrite = true # overwrite the output directory? (true/false/"ask")
+csvoutput = true # save collected data in CSV files
+visualise = true # generate result graphs
+storedata = true # keep collected data in memory
+loglevel = "warn" # verbosity level: "debug", "info", "warn"
+processors = 6 # number of processors to use on parallel runs
+seed = 1 # seed value for the RNG (0 -> random value)
+# dates to start and end the simulation
+startdate = 2022-02-01
+enddate = 2022-03-31
+
+[world]
+mapdirectory = "." # the directory in which all geographic data are stored
+mapresolution = 10 # map resolution in meters
+landcovermap = "landcover_jena.tif" # location of the landcover map
+farmfieldsmap = "fields_jena.tif" # location of the field geometry map
+weatherfile = "weather_jena.csv" # location of the weather data file
+
+[farm]
+farmmodel = "BasicFarmer" # which version of the farm model to use
+setaside = 0.04 # proportion of farm area set aside as fallow
+fieldoutfreq = "daily" # output frequency for crop/field data, daily/monthly/yearly/end/never
+
+[nature]
+targetspecies = ["Wolpertinger", "Wyvern"] # list of target species to simulate - example species
+#targetspecies = ["Skylark"] # list of target species to simulate
+popoutfreq = "daily" # output frequency population-level data, daily/monthly/yearly/end/never
+indoutfreq = "daily" # output frequency individual-level data, daily/monthly/yearly/end/never
+insectmodel = ["season", "habitat", "pesticides"] # which factors affect insect growth ("weather" is not yet implemented)
+
+[crop]
+cropmodel = "aquacrop" # crop growth model to use, "almass", "aquacrop", or "simple"
+# TODO: these two files below are ALMaSS parameters
+cropfile = "crop_data_general.csv" # file with general crop parameters
+growthfile = "almass_crop_growth_curves.csv" # file with crop growth parameters