Add initial support for soil types

- `Pixel` now has a field `soiltype`

- read soil type information from tiff file

src/core/simulation.jl

@@ -134,7 +134,8 @@ function initmodel(settings::Dict{String, Any})
     with_logger(logger) do
         landscape = initlandscape(settings["world.mapdirectory"],
                                   settings["world.landcovermap"],
-                                  settings["world.farmfieldsmap"])
+                                  settings["world.farmfieldsmap"],
+                                  settings["world.soiltypesmap"])
         weather = initweather(joinpath(settings["world.mapdirectory"],
                                        settings["world.weatherfile"]),
                               settings["core.startdate"],

src/parameters.toml

@@ -27,6 +27,7 @@ mapdirectory = "data/regions/jena" # the directory in which all geographic data
 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
+soiltypesmap = "soil.tif" # name of the soil type map in the map directory
 weatherfile = "weather.csv" # name of the weather data file in the map directory
 	
 [farm]

src/world/landscape.jl

@@ -9,9 +9,31 @@ using Printf
 "The land cover classes encoded in the Mundialis Sentinel data."
 @enum LandCover nodata forest grass water builtup soil agriculture
 
+# TODO: check names and order of enum
+# TODO: note where values come from
+## IMPORTANT: do not change the order of this enum, or initlandscape() will break!
+"The soil type"
+@enum SoilType begin
+    soiltype_nodata = 0       # 0
+    soiltype_sand             # 1
+    soiltype_loamy_sand       # 2
+    soiltype_sandy_loam       # 3
+    soiltype_loam             # 4
+    soiltype_silt_loam        # 5
+    soiltype_silt             # 6
+    soiltype_sandy_clay_loam  # 7
+    soiltype_clay_loam        # 8
+    soiltype_silty_clay_loam  # 9
+    soiltype_sandy_clay       # 10
+    soiltype_silty_clay       # 11
+    soiltype_clay             # 12
+    soiltype_unknown          # 13  TODO: invented this to get a type 13
+end
+
 "The types of management event that can be simulated"
 @enum Management tillage sowing fertiliser pesticide harvesting
 
+
 """
     Pixel
 
@@ -21,14 +43,17 @@ in a single object. The model landscape consists of a matrix of pixels.
 """
 mutable struct Pixel
     landcover::LandCover          # land cover class at this position
+    soiltype::SoilType            # soil type class at this position
     fieldid::Union{Missing,Int64} # ID of the farmplot (if any) at this position
     events::Vector{Management}    # management events that have been applied to this pixel
     animals::Vector{Int64}        # IDs of animals currently at this position
     territories::Vector{String}   # IDs of animals that claim this pixel as part of their territory
 end
 
+Pixel(landcover::LandCover, soiltype::SoilType, fieldid::Union{Missing,Int64}) =
+    Pixel(landcover, soiltype, fieldid, Vector{Management}(), Vector{Int64}(), Vector{String}())
 Pixel(landcover::LandCover, fieldid::Union{Missing,Int64}) =
-    Pixel(landcover, fieldid, Vector{Management}(), Vector{Int64}(), Vector{Int64}())
+    Pixel(landcover, soiltype_nodata, fieldid)
 Pixel(landcover::LandCover) = Pixel(landcover, missing)
 
 showlandscape(mat::T) where T <: AbstractMatrix{Pixel} = showlandscape(stdout, mat)
@@ -102,22 +127,26 @@ mutable struct FarmEvent
 end
 
 """
-    initlandscape(directory, landcovermap, farmfieldsmap)
+    initlandscape(directory, landcovermap, farmfieldsmap, soiltypesmap)
 
 Initialise the model landscape based on the map files specified in the
 configuration. Returns a matrix of pixels.
 """
-function initlandscape(directory::String, landcovermap::String, farmfieldsmap::String)
+function initlandscape(directory::String, landcovermap::String, farmfieldsmap::String, soiltypesmap::String)
     @debug "Initialising landscape"
     landcovermap = joinpath(directory, landcovermap)
     farmfieldsmap = joinpath(directory, farmfieldsmap)
-    #TODO replace errors with exception
-    !(isfile(landcovermap)) && Base.error("Landcover map $(landcovermap) doesn't exist.")
-    !(isfile(farmfieldsmap)) && Base.error("Field map $(farmfieldsmap) doesn't exist.")
+    soiltypesmap = joinpath(directory, soiltypesmap)
+    !(isfile(landcovermap)) && error("Landcover map $(landcovermap) doesn't exist.")
+    !(isfile(farmfieldsmap)) && error("Field map $(farmfieldsmap) doesn't exist.")
+    !(isfile(soiltypesmap)) && error("Soil type map $(soiltypesmap) doesn't exist.")
     # read in TIFFs, replacing missing values with 0
     landcover = coalesce(GeoArrays.read(landcovermap), 0)
     farmfields = coalesce(GeoArrays.read(farmfieldsmap), 0)
-    (size(landcover) != size(farmfields)) && Base.error("Input map sizes don't match.")
+    soiltypes = coalesce(GeoArrays.read(soiltypesmap), 0)
+    if size(landcover) != size(farmfields) || size(farmfields) != size(soiltypes)
+        error("Input map sizes don't match.")
+    end
     width, height = size(landcover)[1:2]
     landscape = Matrix{Pixel}(undef, width, height)
     for x in 1:width

test/nature_tests.jl

@@ -47,7 +47,7 @@ end) # end eval
 @testset "Habitat macros" begin
     # set up the testing landscape
     model = inittestmodel()
-    model.landscape[6,6] = Pixel(Ps.agriculture, 1, [], [], [])
+    model.landscape[6,6] = Pixel(Ps.agriculture, Ps.soiltype_sand, 1, [], [], [])
     fp = Ps.FarmPlot(
         1, [(6,6)], 1,
         Ps.ALMaSS.CropState(

test/runtests.jl

@@ -37,7 +37,8 @@ function inittestmodel(smallmap=true)
     else
         landscape = Ps.initlandscape(TESTSETTINGS["world.mapdirectory"],
                                      TESTSETTINGS["world.landcovermap"],
-                                     TESTSETTINGS["world.farmfieldsmap"])
+                                     TESTSETTINGS["world.farmfieldsmap"],
+                                     TESTSETTINGS["world.soiltypesmap"])
     end
     weather = Ps.initweather(joinpath(TESTSETTINGS["world.mapdirectory"],
                                       TESTSETTINGS["world.weatherfile"]),

test/test_parameters.toml

@@ -26,6 +26,7 @@ mapresolution = 10 # map resolution in meters
 landcovermap = "landcover.tif" # location of the landcover map
 farmfieldsmap = "fields.tif" # location of the field geometry map
 weatherfile = "weather.csv" # location of the weather data file
+soiltypesmap = "soil.tif" # name of the soil type map in the map directory
 
 [farm]
 farmmodel = "BasicFarmer" # which version of the farm model to use

test/test_parameters_almass.toml

@@ -26,6 +26,7 @@ mapresolution = 10 # map resolution in meters
 landcovermap = "landcover.tif" # location of the landcover map
 farmfieldsmap = "fields.tif" # location of the field geometry map
 weatherfile = "weather.csv" # location of the weather data file
+soiltypesmap = "soil.tif" # name of the soil type map in the map directory
 
 [farm]
 farmmodel = "BasicFarmer" # which version of the farm model to use

test/test_parameters_simplecrop.toml

@@ -26,6 +26,7 @@ mapresolution = 10 # map resolution in meters
 landcovermap = "landcover.tif" # location of the landcover map
 farmfieldsmap = "fields.tif" # location of the field geometry map
 weatherfile = "weather.csv" # location of the weather data file
+soiltypesmap = "soil.tif" # name of the soil type map in the map directory
 
 [farm]
 farmmodel = "BasicFarmer" # which version of the farm model to use