diff --git a/src/core/simulation.jl b/src/core/simulation.jl
index e57f075b9b722ae78c920bd058b4a40d50413523..5c3628eb62e5e7b6970094e1df3bc77bba149448 100644
--- a/src/core/simulation.jl
+++ b/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"],
diff --git a/src/parameters.toml b/src/parameters.toml
index b743ca5882c71730f481206826029e56320d45c5..743343aeacda665e22907746aa20ffedb3afa230 100644
--- a/src/parameters.toml
+++ b/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]
diff --git a/src/world/landscape.jl b/src/world/landscape.jl
index 4d3e2b67572c494dd45a6f264822adee56da03be..cfc26691f349dd300526e255742db7515e94aa6d 100644
--- a/src/world/landscape.jl
+++ b/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
diff --git a/test/nature_tests.jl b/test/nature_tests.jl
index 0a786cabd2cd81eacd45fd12d9930ee561d57566..1e829fc9d7b6f9e30c542dce8c2741759ba57560 100644
--- a/test/nature_tests.jl
+++ b/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(
diff --git a/test/runtests.jl b/test/runtests.jl
index a655c01849151bd01c0fdffea176676d1c11ac1e..906be637f6631b120b328d326a8b7fc2cc122f18 100644
--- a/test/runtests.jl
+++ b/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"]),
diff --git a/test/soil_jena.tif b/test/soil_jena.tif
new file mode 100644
index 0000000000000000000000000000000000000000..0ffb7317a7702a1da2289f4805c780d0d199221e
Binary files /dev/null and b/test/soil_jena.tif differ
diff --git a/test/test_parameters.toml b/test/test_parameters.toml
index abc3dab62641bf5518f6f65d497d7bf64cd883e2..fa539d902d2257572de15218d4bb2456822072fc 100644
--- a/test/test_parameters.toml
+++ b/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
diff --git a/test/test_parameters_almass.toml b/test/test_parameters_almass.toml
index be08e4af053ca7368bd0418b937f624d3ccc6781..6d8fa7efb4d6e821fdb875018460b93c91f856dc 100644
--- a/test/test_parameters_almass.toml
+++ b/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
diff --git a/test/test_parameters_simplecrop.toml b/test/test_parameters_simplecrop.toml
index 1bcf23b8d91bba31d8840cea91d00425e75c8773..6540f38f3552df351b8c367e5fe20a9404f73963 100644
--- a/test/test_parameters_simplecrop.toml
+++ b/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