- Email: [email protected]
Digital soil mapping of arable land in Sweden – Validation of performance at multiple scales
Geoderma xxx (xxxx) xxx–xxx
Contents lists available at ScienceDirect
Geoderma journal homepage: www.elsevier.com/locate/geoderma
Digital soil mapping of arable land in Sweden – Validation of performance at multiple scales Kristin Piikki⁎, Mats Söderström Department of Soil and Environment, Swedish University of Agricultural Sciences (SLU), P.O. Box 234, SE-53223 Skara, Sweden
A R T I C L E I N F O
A B S T R A C T
Handling Editor: A.B. McBratney
In this study, we produced a detailed digital soil map of topsoil texture and soil organic matter (SOM) content for 2.4 million ha of arable land in Sweden (DSMS). Three spatially exhaustive datasets (a laser-scanned digital elevation model, airborne gamma radiation scanning data and a legacy Quaternary deposit map) were calibrated against topsoil texture and SOM content in around 13,500 soil samples, using multivariate adaptive regression splines (MARSplines) modelling. We then deployed the MARSplines models to produce raster maps (50 m × 50 m) of clay, sand and SOM content. The modelling procedure was validated by an independent dataset of about 24,000 samples clustered on 544 farms (with a local sample density of one per 3 ha). The error in clay content was < 8% in 75% of the validation samples, while for sand content and SOM content it was 13% and 2%, respectively. Corresponding values for the farm-average level were 6%, 11% and 2%, respectively. Modelling efficiency values revealed that the clay content map was a considerable improvement over the mean of the reference values at national level, regional level and, in many cases, also farm level. However, SOM content predictions showed little or no improvement over the mean of the reference samples (at any scale), due to poor correlation with the exhaustive predictor variables at all three scales investigated. The DSMS soil geodatabase will continue to be improved and have more layers added, e.g. derived layers calculated from the primary clay, sand and SOM layers by use of pedotransfer functions. Practical use of DSMS is exemplified here in an internet application for deriving prescription files for precision agriculture.
Keywords: Texture Organic matter Digital elevation model Gamma radiation MARSplines Quaternary deposit
1. Introduction Ever since McBratney et al. (2003) outlined an explicit framework for digital soil mapping (DSM), there has been a virtual explosion in the availability of spatial data and pedometricians world-wide creatively applying new computational techniques in order to produce detailed soil maps (e.g. Arrouays et al., 2014; De Brogniez et al., 2015; Hengl et al., 2017; Stoorvogel et al., 2017). DSM involves determining soil variations in relation to the landscape, finding measurable covariates for the soil property of interest and developing quantitative models for prediction of the target property. However, while a number of DSM products provide sufficiently high spatial resolution to give the impression that they are applicable for management decisions at farm level, validations with local data often show very large errors, indicating that the data are not suitable for use in precision agriculture (Söderström et al., 2016, 2017). An effort was recently launched in Sweden to produce a freely available digital map of topsoil texture (Digital Soil Map of Sweden, DSMS) for most of the country's arable land (about 90% or ~2.4 million ha) at an accuracy sufficient for use in
⁎
precision agriculture. This might be possible, since three meticulous national datasets describing the landscape are available and could serve as covariates in modelling. These are: a digital elevation model (DEM) based on airborne laser scanning, a national airborne gamma radiation survey and a Quaternary deposit (QD) map. There is also a newly available reference dataset of soil samples (national 1-km grid) funded by the government in order to improve estimates of nutrient loads to the Baltic Sea (Paulsson et al., 2015). In addition, a large dataset of geographically positioned soil analyses on 544 farms is available for validations at farm scale. A recent pilot study covering about 100,000 ha of intensively cultivated land in south-west Sweden reported root mean square error of prediction of 5.7% for clay, 12.5% for sand and 1.3% for soil organic matter (SOM) (Söderström et al., 2016). Compared with maps of soil texture produced from farmers' own data using contemporary techniques (interpolated maps based on one sample per 3 ha), the regional DSM produced in the pilot study was equally good or better. It was based on the DEM, gamma radiation survey, QD map and reference soil samples mentioned above. The model type used was multivariate adaptive regression splines
Corresponding author. E-mail addresses: [email protected] (K. Piikki), [email protected] (M. Söderström).
http://dx.doi.org/10.1016/j.geoderma.2017.10.049 Received 26 April 2017; Accepted 27 October 2017 0016-7061/ © 2017 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/BY/4.0/).
Please cite this article as: Piikki, K., Geoderma (2017), http://dx.doi.org/10.1016/j.geoderma.2017.10.049
Geoderma xxx (xxxx) xxx–xxx
K. Piikki, M. Söderström
of within 5 m (various brands). Data on these soil samples are spatially complementary to, and were merged with, data obtained in a previous, spatially sparser, soil sampling campaign undertaken during 2001–2007 by the Swedish Environmental Protection Agency (SNV; Naturvårdsverket, Stockholm, Sweden) (Eriksson et al., 2010). Only soil samples with SOM ≤20% were included in the present modelling work. In total, 13,485 soil samples constituted the grid dataset used for calibration of DSM prediction models (11,734 SJV samples and 1751 SNV samples). For local-scale validation, we used an entirely independent dataset collected on farms located mostly in intensely cultivated districts of Sweden, within an environmental support subsidies programme funded through SJV (years 2007–2009). These samples were distributed within farms in a grid pattern, with a density of one sample per three hectares. In this dataset, sample positions were available but farm identities were not. Therefore, synthetic farms were generated through a spatial clustering procedure. Samples within 1000 m of each other were grouped into ‘farms’ and only those with at least five samples were kept for validation. This resulted in a validation set of 24,323 analyses on 544 farms. The median size of these farms was 126 ha. All soil samples in the grid dataset and the farm dataset were taken from the topsoil (0–20 cm; 7–10 topsoil cores within a 3–5 m radius circle at each sampling point were bulked). The samples were air-dried at 35–40 °C, milled and sieved through a 2-mm mesh. The content of sand (0.2–2 mm) was determined by sieving and weighing and the content of clay (< 2 μm) by the sedimentation method (method ISO 11277; Gee and Bauder, 1986). The sand and clay content were both expressed as a percentage of the fine soil fraction (< 2 mm). The SOM content was determined from the clay content and the loss on ignition (Eq. (1); original method by Ekström, 1927). The laboratory analyses were carried out by Eurofins laboratory (Kristianstad, Sweden).
(MARSplines). Appropriate validation of DSM is a challenge. The conventional scale concept used for printed maps is not directly applicable for digital soil maps. Malone et al. (2013) list i) extent, ii) resolution and iii) support as appropriate scale specifications of digital soil maps in raster format. In the present study, validations were carried out at two support (soil samples and farms) and three extents (farm extent, regional extent and national extent). The aims of this study were to: 1) Describe the methods used for production of DSMS, which shows topsoil texture and SOM for arable land. 2) Present the results of validations of DSMS at multiple scales, focusing on assessing local, independent within-farm map performance. 3) Suggest practical user applications based on DSMS. Aim (2) deviates from that for most DSM products, which are generally validated through various types of cross-validation using subsets of the soil mapping dataset. 2. Materials and methods 2.1. Area of interest The mapped area comprises approximately the southern third of Sweden, which includes about 90% of the agricultural land in the country (Fig. 1b). A spatial database of cultivated arable land derived from the EU subsidies database of the Swedish Board of Agriculture (SJV; Jordbruksverket, Jönköping, Sweden), the so-called ‘Block Map’ (from 2013), was used to delineate the area of interest map in detail (example shown in Fig. 1a). Areas of organic soils identified from the 1:50,000 QD map by the Geological Survey of Sweden (SGU, Uppsala, Sweden) and of wetlands identified from the Block Map 2013 were classified as organic, and were excluded from the modelling. The reason for excluding organic soils from the map is that soil samples with SOM content > 20% are only rarely analysed for soil texture, so there were insufficient reference data for mapping areas with this type of soil.
SOM = Loss on ignition − 0.46 − 0.047 × clay content
(1)
2.3. Covariates Three types of spatially exhaustive ancillary information sources were used in the predictive modelling: airborne gamma-ray spectrometry data, elevation data acquired through airborne laser-scanning and information from legacy QD maps. All predictor data were calculated for (or simply extracted to) the centre point of a 50 m × 50 m raster mesh. The pre-processing of covariate data followed the procedure adopted in the pilot study to this DSM work (Söderström et al., 2016). In order to handle the very large amount of data, the map area
2.2. Soil samples and soil analysis methods During 2011–2012, SJV implemented a national soil sampling campaign, resulting in a more or less regular grid of one soil sample per km2 across Swedish agricultural land. Positions were determined by global positioning system (GPS) with an estimated positional accuracy
Fig. 1. Spatial distribution of the soil datasets used. Map a) shows the spatial distribution of the farm and grid datasets over the agricultural fields (blocks) in a small area and map b) shows the entire extent of the digital soil mapping area in southern Sweden and the distribution of agricultural land (2.4 million ha) and lakes.
2
Geoderma xxx (xxxx) xxx–xxx
K. Piikki, M. Söderström
recorded concentrations of K < 0% or Th < 1 ppm were omitted, since these often indicate data affected by water bodies or peat. A small number of abnormally high values of K and Th were also present in the dataset and are detrimental for construction of robust calibration models, even if they are correct. In many cases it is not variations in soil texture that cause such extremes, but rather anomalies in the geological conditions (e.g. certain bedrock outcrops). Therefore, K concentrations > 4.5% were set to 4.5%, and Th concentrations > 30 ppm were set to 30 ppm. In some rare cases, abrupt discrepancies in values were noted between flight lines, but in most cases it was possible to adjust the general level between such measurements manually. If not, such data were deleted. The filtered K and Th data were then interpolated to the 50 m × 50 m mesh by ordinary block kriging. 2.3.2. Elevation data The national digital elevation model derived from airborne laser scanning and supplied by LMV as a raster with 2-m spatial resolution (Grid 2+, LMV, Gävle, Sweden) was used as a data source for elevation and landform. The data were resampled (averaged) to 10 m × 10 m, which was considered more appropriate for the present purpose. Various measures of landform were estimated as topographic position index (TPI; Weiss, 2001), i.e. the relative elevation of raster cells compared with circular neighbourhoods of 5 ha, 50 ha and 500 ha (denoted TPI5, TPI50 and TPI500, respectively) (Fig. 3).
Fig. 2. Sub-regions used for modelling.
2.3.3. Legacy soil data on Quaternary deposits The Weichselian glaciation, the melting of the inland icecap and the subsequent sea-level fluctuations around the Scandinavian Peninsula in the Holocene have determined the overall distribution of the regolith in Sweden. The genetics of these young sediments in Sweden are described by geologists in QD maps, scale 1:50,000, which are available as vector polygon maps from SGU. More than 100 different deposit types have been assigned. To transfer these data into useable covariates, in this work we simplified the deposit classification into two different versions:
was manually divided into 22 sub-regions based on expert knowledge (Fig. 2), in a process that was a trade-off in terms of uniformity of geomorphology and geology and amount of soil observations. This process followed an initial test in which data from six different parts of the study area with different physical characteristics (for example one hilly region with granites known for their high level of radiation, another area mostly consisting of fine post-glacial sediments in a relatively flat landscape) were analysed separately and also jointly. The separate models performed somewhat better, indicating that such an approach was preferable. In some cases, when the physical characteristics were less distinct and the soil observations otherwise were judged to be too few or spatially skewed in a sub-region, or to reduce the effect of slight differences in predicted values between two adjacent regions, the sub-regions overlapped to some extent.
• Seven classes (QD7): Clay, Sand, Silt, Till, Till clay, Organic and Bedrock • Three classes (QD3): Clay, Organic and Other The ‘Organic’ class, which was the same in the QD7 and the QD3 classifications, was used to mask organic areas from the map extent and no samples from areas classified as organic were included in the model parameterisation. We also simplified the geometry of the QD map, by dissolving polygons smaller than 1 ha with the polygon with which they shared their longest border. The reason for this is that the size of such small features is often exaggerated for cartographic reasons in QD maps and this, in combination with positional uncertainty about the extent, makes them unsuitable to include in modelling.
2.3.1. Gamma-ray spectrometry data Natural emissions of gamma radiation in Sweden have been measured by airborne gamma-ray spectrometry since the 1960s by SGU. Data are collected along flight lines normally separated by 200 m, from an altitude of 30 m or 60 m, by recording the concentration of the naturally occurring isotopes 232Th (ppm; eTh), 40K (ppm) and 238U (%; eU) every 16–17 m along the flight lines. The data are then normalised by SGU to take into account variations in atmospheric radiation, topographical influence and variations in altitude (SGU, personal communication). The footprint radius of recorded values on the ground is approximately four times the altitude, and about 80% of the response originates from the upper 0.3 m of soil (IAEA, 2003). The usefulness of gamma-ray data in predictive mapping of soil properties in northern Europe has been demonstrated in several earlier studies, e.g. by Piikki et al. (2013), Söderström et al. (2016) and Van der Klooster et al. (2011). The eU shows considerable short-range variation and smallscale noise in parts of the dataset (Söderström and Eriksson, 2013; Söderström et al., 2016) and was omitted from the calibrations in the present analysis to reduce the risk of erroneous and over-fitted prediction models. Further filtering of the gamma-ray data included removal of observations < 100 m from water bodies (as delimited by the 1:50,000 topographical map of the Land Survey of Sweden (LMV; Lantmäteriet, Gävle, Sweden)) and > 300 m from arable land. To remove outliers and data considered non-representative for arable land,
2.4. Spatial predictive modelling MARSplines, which is a data mining method used to parameterise statistical relationships between the variables in a dataset (X; predictors Y; response variables), was used for modelling. A detailed description of MARSplines models is given by Hastie et al. (2009). In brief, such a model consists of a set of piece-wise, univariate linear regression functions (so-called basis functions) that are defined above or below threshold values (breakpoints) of the predictor in question. MARSplines models can be additive or allow for interactions among the basic functions. We chose to use only additive models, in order to minimise the risk of overfitting associated with more complex regression models. Before the model parameterisation, factor expansion was performed on the categorical data (QD3 and QD7). Each class was transformed to one binary variable, where 1 indicated that the sample or raster cell belonged to the QD class in question, while 0 indicated the opposite. In 3
Geoderma xxx (xxxx) xxx–xxx
K. Piikki, M. Söderström
Fig. 3. Covariates used in modelling. a) Elevation z above sea level; b–d) relative topography TPI500, TPI50 and TPI5, i.e. the elevation relative to the average of a neighbourhood area of 500 ha, 50 ha and 5 ha, respectively; e) concentration of 232Th; f) concentration of 40K; g) Quaternary deposits simplified into seven classes (QD7); and h) Quaternary deposits simplified into three classes (QD3). In e–h), data are shown only for areas of arable land, with elevation as background. Blue areas indicate lakes in all maps. The overview map show the extent of the digital soil map of Sweden and the filled circle indicate the location of the 10 km × 10 km area in maps a–h. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
farm dataset. The following statistics were calculated: mean absolute error (MAE; Janssen and Heuberger, 1995) and Nash-Sutcliffe model efficiency (E; Nash and Sutcliffe, 1970). The latter indicates how well a plot of observed values and values predicted by a model fits the 1:1line, and is a useful index for testing whether the map values are more accurate than the average of the reference samples. The value of E varies between −∞ and 1, where a negative value indicates that the map is less accurate than the average of the reference samples, E = 0 means that the map is just as accurate as the reference sample average, a positive value of E means that the map is an improvement over the average and a value of 1 indicates that the map values at the soil sampling locations are equal to the soil analysis results, i.e. perfect correlation.
order to construct as robust models as possible, we also used a built-in cropping function to remove all basis functions that did not significantly contribute to the prediction. In summary, the following eight covariate predictors (Fig. 3) were used to derive models for clay, sand and SOM content:
• Concentrations of naturally occurring radioactive isotopes: and K • Elevation and relative topography: z, TPI , TPI and TPI • Quaternary deposits: QD3 and QD7
232
Th
40
5
50
500
Models were calibrated by the reference data (the grid soil dataset combined with covariate predictors) in each sub-region (Fig. 2) and then deployed onto the 50 m × 50 m mesh for that area. The predictions in all 22 sub-regions were then merged to a final raster dataset covering the entire map area. In some cases, neighbouring sub-regions overlapped. In the overlapping areas the final prediction was a distance-weighted average of the two predicted values. The distance was calculated from the raster cell in question to the edge of the overlapping area. That yielded a smooth transition from one model area to the next.
2.7. Software All spatial analyses were conducted in ArcGIS with the extensions Spatial Analyst and Geostatistical Analyst (ESRI, Redlands, CA, USA), while the MARSplines modelling was performed using R (R Core Team, 2017), package Earth (Milborrow, 2017). Statistica 13 (Dell Inc., Tulsa, OK, USA) was also used for data management and analyses.
2.5. Multiple-scale validation 3. Results Validation statistics were calculated for all orthogonal combinations of the two support levels and the three spatial extents, except at the farm extent, for which map accuracy could not be tested at farm support (since n = 1 in each validation set). This means that at the national extent one value of E and one value of MAE was calculated. In the regional case, 22 values of MAE and E were calculated and subsequently averaged across all modelling areas, while at farm level 544 values of MAE and E were calculated and then averaged across all farms.
3.1. Descriptive statistics The quartiles of the frequency distributions were similar for the calibration data (grid dataset) and the validation data (farm dataset) for all three map variables (Table 1.) As can be inferred from the interquartile ranges, the sand content in the laboratory-analysed samples was right-skewed, while the clay content was somewhat left-skewed and the frequency distribution of SOM had an extreme skewness to the left. At this national level, the frequency distributions of the predicted soil properties were very similar to the frequency distributions of their laboratory-analysed counterparts. However, there were deviations at the right tail of the distributions (especially for the SOM content), indicating that under-prediction of high values was more common than
2.6. Validation statistics Data from the prediction rasters (clay content, sand content and SOM content) were extracted for all validation sample locations in the 4
Geoderma xxx (xxxx) xxx–xxx
K. Piikki, M. Söderström
Table 1 Descriptive statistics on the response variables: values for individual analyses of the grid dataset (n = 13,485), used for calibration, and for the farm dataset, used for validation. For the latter, statistics for the individual sample values (n = 24,323) are given before the slash and for farm-average values (n = 544) after the slash. SOM = soil organic matter. Dataset
Value
Variable
Minimum
1st Quartile
Median
3rd Quartile
Maximum
Grid Grid Grid Farm Farm Farm Farm Farm Farm
Measured Measured Measured Measured Measured Measured Predicted Predicted Predicted
Clay content [%] Sand content [%] SOM content [%] Clay content [%] Sand content [%] SOM content [%] Clay content [%] Sand content [%] SOM content [%]
2 0 0 1/1 1/4 0/1 0/3 0/4 0/1
10 14 3 12/12 14/20 2/3 13/13 18/22 1/1
19 38 4 20/18 37/43 3/3 21/18 37/44 2/2
33 59 5 36/28 58/56 4/4 35/28 54/55 2/2
80 96 20 73/58 96/91 19/17 63/56 93/81 7/4
Table 2 Cross-tables of Pearson correlation coefficient (r) among all pairs of continuous variables. The r values were calculated and averaged at: a) national level; b) sub-region level and c) farm level. 232
Th
40
TPI50
TPI500
Clay
Sand
SOM
– 0.69 −0.27 0.26 −0.20
– − 0.23 0.27 − 0.18
– − 0.85 0.16
– − 0.22
–
b) Average of the pair-wise correlation coefficients calculated among all samples in each of the model areas 232 Th – 40 K 0.32 – z − 0.18 −0.05 – RT5 0.02 0.08 0.02 – RT50 − 0.05 0.12 0.10 0.62 – RT500 − 0.13 0.09 0.26 0.22 0.57 Clay 0.60 0.12 − 0.26 − 0.05 −0.13 Sand − 0.53 −0.12 0.23 0.00 0.09 SOM 0.06 −0.25 0.09 − 0.15 −0.17
– − 0.21 0.22 − 0.14
– − 0.79 0.20
– − 0.24
–
c) Average of the pair-wise correlation coefficients calculated among all samples on each of the farms Th – K 0.33 – z − 0.12 0.11 – RT5 0.05 0.09 0.17 – RT50 − 0.01 0.16 0.39 0.61 – RT500 − 0.08 0.14 0.66 0.29 0.66 Clay 0.35 0.03 − 0.19 − 0.09 −0.17 Sand − 0.35 −0.01 0.21 0.05 0.15 SOM 0.03 −0.24 − 0.16 − 0.18 −0.22
– − 0.18 0.21 − 0.19
– − 0.76 0.23
– − 0.30
–
a) Pair-wise correlation coefficients 232 Th – 40 K 0.41 z − 0.26 TPI5 − 0.16 TPI50 − 0.25 TPI500 − 0.25 Clay 0.83 Sand − 0.76 SOM 0.08
K
z
TPI5
calculated among all samples in the national dataset – 0.01 0.00 0.02 0.00 0.27 −0.26 −0.18
– 0.08 0.15 0.19 − 0.24 0.13 0.11
– 0.65 0.31 − 0.19 0.17 − 0.16
232 40
z = elevation above sea level, TPI = topographic position index.
soils than in soils with a high clay content because the clay particles protect the SOM from microbial action (e.g. Christensen, 1987; van Veen and Kuikman, 1990). When quantifying and interpreting the correlations between SOM and the predictor variables, one should bear in mind that SOM values > 20% were excluded from the dataset, which means that the correlations are evaluated for mineral soils only and not for organic soils.
over-prediction. Correlation matrices of Pearson correlation coefficients (r) were calculated for three different spatial extents (Table 2). It is possible that different processes govern covariance between two variables at different scales. A weaker correlation (smaller absolute value of r) at finer scale may also be due to the range of variation and the number of observations being smaller when a smaller area is considered. Therefore, two variables can be well correlated at the national scale, but show a weaker correlation for regional model areas or farms. In the present case, for example, there was much lower correlation between Th and clay content at farm scale (Table 2c) compared with the entire area (Table 2a) or compared with the sub-region level (Table 2b), probably due to the fact that farms with small standard deviation in Th and clay content reduce the average r value at farm scale. The SOM content showed weak correlations to all predictors used in this study, but was most strongly related to the content of sand. This correlation had a negative sign at all three extents. This may actually be a causal relationship, since mineralisation of organic matter is faster in sandy
3.2. Absolute map errors The graphs in Fig. 4 visualise the MAE distributions in the predicted soil properties. As expected, the higher percentiles were smaller for the farm-wise aggregated data (solid lines) than for the individual samples (dotted lines). The error in clay content was < 8% in 75% of the validation samples, while for sand content it was 13% and for SOM content it was 2%. Corresponding values for farm averages were 6%, 11% and 2%, respectively. This kind of validation at two support levels is important for assessment of possible practical applications of the DSMS 5
Geoderma xxx (xxxx) xxx–xxx
K. Piikki, M. Söderström
Fig. 4. Absolute error. Distribution of mean absolute error (MAE) of predicted soil properties: a) clay content, b) sand content and c) soil organic matter (SOM) content for individual soil samples (n = 24,323) and farm averages of these soil samples (n = 544). Relative frequency on the y-axis is the fraction of samples within discrete intervals of the x-axis (width of 5%units for clay and sand and 1%-unit for SOM). The integrated area under the each curve is equal to 1. Fig. 5. Maps of predictions and prediction errors. Predicted a) clay content and b) sand content of the digital soil map of Sweden (DSMS) and (c, d) geographical distribution of the magnitude of error calculated for 25 km × 25 km cell areas. As 70% of the area of Sweden is forested, the coverage appears sparse in the map, even though most of the arable land (~90%) is mapped. MAE = mean absolute error.
land that was not mapped in the current version of DSMS (one small 10 km × 1 km area) is located there. The reason for this was missing airborne gamma data. In the agricultural districts in the south and south-west of the country, the clay content prediction errors were relatively small (< 5%). The smallest clay prediction errors (in general about 2.5%) were found in the central part of southern Sweden. In this region, both the clay content and the magnitude of its variation are low. For sand content (Fig. 5d), the variation did not follow a similarly clear spatial pattern. In the districts with large areas of contiguous arable
data. For precision agriculture, point location validation is necessary, whereas e.g. for regional modelling of nutrient leaching, farm-level aggregated data may be adequate. Fig. 5 shows the geographical distribution of MAE for the clay content and sand content. Large clay content errors coincided generally with areas of high clay content in the topsoil, e.g. in the east and north-east parts of the map area (Fig. 5c). Particularly high errors for clay were recorded in the east-central region, where the airborne gamma radiation scanning also showed some quality problems in places. In fact, the only area of arable (non-organic)
6
Geoderma xxx (xxxx) xxx–xxx
K. Piikki, M. Söderström
support). In all cases, the E value was higher for farm average values than for individual samples. However, as Fig. 4 shows, the MAE for individual samples did not differ much from that of farm-averaged values, although the error was more heterogeneous and showed a larger span at sample support level than at farm support level. Based on the MAE of SOM content (Fig. 4c), the DSMS SOM layer could be deemed acceptable for regional applications, but the scatter plots (Fig. 6) and the E values (Table 3) indicate a weak correlation between the map and the reference samples, indicating that an average of all samples in a region may be almost as accurate as the DSMS SOM layer. The reason that SOM content could not be mapped well enough to be considered very useful at the local scale is most likely that it was poorly correlated to the exhaustive predictors (Table 2). Even if advanced data mining methods are used, it is not possible to parameterise relationships that are not present in the dataset (causal or not). Plausible reasons why no useful relationship could be parameterised between the available predictors and SOM content are that the small-scale (between-field) variation in SOM is highly influenced by management, that the local magnitude of SOM variation is often small, that the laboratory analysis uncertainty is relatively high in relation to the variation, and that topsoil gamma radiation is not very much affected by the SOM content unless the SOM content is high. As an example of the latter, Berglund and Berglund (2010) described a method for delineating peat areas using the activity concentration of 40K as a predictor. Since all samples for calibration and validation in the present study were obtained during a 12-year period, it is possible that the SOM content had changed to some extent, at least in fields where management had changed, adding further difficulty in predicting variation. At national and sub-region level, DSMS was an improvement over using the mean of the reference samples (both farm-average values and point location data, SOM content excluded; Table 3). For individual farms the average E values were closer to 0, which indicates that DSMS is not always an improvement over the soil sample mean. However, there was considerable standard deviation in the E values, leading to the interpretation that for many farms, but not all, DSMS is an improvement over the reference sample mean. DSMS data can be used as input data in national, overall modelling of phosphorus (P) leaching from agricultural soils, which is one of the reasons why donors funded this project. It can also be used to target control measures against P leaching to parts of fields where they would be most effective. This would be a considerable improvement for environmental authorities, which otherwise can often only decide on very general mitigation actions for large areas to control e.g. leaching of nutrients. Such measures are likely to be more efficient if they can be targeted at high-risk areas and hotspots. For example, structure liming is recommended to be applied to soils with > 15% clay (Albertsson et al., 2016). DSMS data can also be used as the basis for generating variable-rate application (VRA) files aimed at practical precision agriculture. In an interactive internet application (www.markdata.se; a free-to-use web application in Swedish language), DSMS clay content is classified by the user, and outputs such as VRA files of seed rate can be produced and downloaded. Most users choose to reduce the seed density in areas with a coarser texture, where germination is often better, in order to achieve an optimal canopy in all parts of the field. Products such as DSMS can also be expanded to include derivative products through pedotransfer functions. Bulk density, water content at field capacity and at wilting point and porosity are examples of properties that can be calculated from DSMS data (Rawls et al., 2003; Kätterer et al., 2006). In 2017, the DSMS SOM layer was used as input to the Swedish contribution to the Global Soil Organic Carbon (GSOC) map (FAO, 2017). In that case, pedotransfer functions were used to estimate bulk density from SOM (Eriksson et al., 2010), and the SOC stock data were aggregated to a 1 km × 1 km raster. Another application of DSMS involving pedotransfer functions is in development of map layers useful for determining the amount of lime required for increasing soil pH in acid soils. According to guidelines from the Swedish Board of
Table 3 Modelling efficiency (E) of texture and soil organic matter (SOM) predictions calculated for point locations and points aggregated to farm averages at three scales (farm, subregion and entire digital soil map of Sweden (DSMS) area) (mean ± standard deviation among farms or modelling areas).
Clay Sand SOM
Farms
Model areas
Entire DSMS area
Point locations
Point locations
Farm averages
Point locations
Farm averages
0.18 ± 0.22 0.18 ± 0.21 0.05 ± 0.07
0.40 ± 0.22 0.37 ± 0.23 0.07 ± 0.06
0.60 ± 0.26 0.54 ± 0.28 0.07 ± 0.07
0.76 0.75 0.12
0.90 0.86 0.14
land, prediction errors of sand content were moderate (10–13%), whereas the largest errors (up to 15% and above) were found in areas where agricultural land is interspersed with forested areas and fields in general are smaller. 3.3. Relative map accuracy An alternative to using regional digital soil maps is to use averages of samples at the farm level or possibly averages for larger regions. Table 3 shows the effects of this on the performance of texture and SOM modelling of aggregation of predictions made at a finer scale. In general, models for clay content performed better than those for sand content. For SOM, the models performed less well. It is clear that aggregation resulted in much more accurate predictions, e.g. the E value for clay content of farms in the entire validation dataset (E = 0.90) indicates a very accurate and precise map. However, looking at finer scale, it is obvious that the predictions were not always successful in describing the variation within the mapped area. The E values consistently increased with increasing extent and also with increasing support. 4. Discussion It is a challenge to present relative validation statistics for a largeextent, high-resolution dataset like DSMS in a manner that is useful for different possible applications of the data (in terms of resolution and extent). This is important to consider when reporting accuracies from validations, and also for users who wish to understand how well a prediction model works, i.e. how accurate a digital soil map is. Correlation analyses at one extent are not directly transferrable to another (Table 2). The same is true for E (Table 3; Fig. 6), which was always smaller at farm level than for the entire modelling area, probably because of the generally much smaller ranges of variation in soil properties and a smaller number of observations in smaller areas (Söderström et al., 2016). Moreover, testing for different extents might indicate that causal and/or indirect relationships between the predictors and the response variables may be governed by different processes operating at different scales. Hence, the digital soil maps produced may accurately describe this at one spatial level, but not at another. For example, the coarse variation at the national scale (e.g. governed by large scale geological variation) may be well described by DSMS data, but the data may fail to describe the small-scale variation within farms or sub-regions (e.g. governed by agricultural management). The rationale behind using two aggregations of the validation soil samples (individual samples and farm average values) was that some covariate data used in the prediction models were coarse. For example, the footprint of the gamma radiation data was much larger than the support of an individual soil sample. Individual samples will therefore not necessarily be the most appropriate resolution for applications. Here we tested whether the accuracy of the predictions was better for farm averages than for describing within-field variation (sample 7
Geoderma xxx (xxxx) xxx–xxx
K. Piikki, M. Söderström
Fig. 6. Scatterplots of predictions against laboratory values. The plots show individual soil samples (a–c; n = 24,323; a random subset of 2000 samples are plotted) and farm averages of these soil samples (d–f; n = 544). SOM = soil organic matter.
include local adaptation through interactive functionality where users can upload their own data, and the application of pedotransfer functions to produce important map layers of soil functional properties.
Agriculture, the models applied in Sweden for determining both the recommended suitable soil pH (target-pH) and the buffering capacity are a function of soil clay content and SOM (Gustafsson, 1999). A test of producing such map layers in DSMS was reported by Söderström et al. (2016) for use at the field level and a method was presented to improve the accuracy of digital soil maps locally through addition of locally available data. Such a procedure was recently further developed and assessed (Piikki et al., 2013), and an automated interactive method for users to upload their own data in a web application to remodel the DSMS at farm scale was presented. With such a system in place, conventional soil mapping (with one sample per 3 ha) was better than the digital soil map for less than one-third of the validated farms.
Acknowledgements This work was financially supported by the Swedish National Space Board (SNSB) [dnr 214/13] and the Geological Survey of Sweden (SGU) [project agreement: 15-03-31]. We acknowledge valuable contributions from the Swedish Board of Agriculture (JV), The Rural Economy and Agricultural Societies (HS), Precision Agriculture Sweden (POS) and many interested farmers. References
5. Conclusions
Albertsson, B., Börling, K., Kvarmo, P., Listh, U., Malgeryd, J., Stenberg, M., 2016. Rekommendationer för gödsling och kalkning 2017. Report JO16 24 Swedish Board of Agriculture, Jönköping, Sweden 102 pp. Available online: http://webbutiken. jordbruksverket.se/sv/artiklar/jo1624.html (verified: October 26, 2017). Arrouays, D., Grundy, M.G., Hartemink, A.E., Hempel, J.W., Heuvelink, G.B., Hong, S.Y., et al., 2014. Chapter three - GlobalSoilMap: toward a fine-resolution global grid of soil properties. Adv. Agron. 125, 93–134. Berglund, Ö., Berglund, K., 2010. Distribution and cultivation intensity of agricultural peat and gyttja soils in Sweden and estimation of greenhouse gas emissions from cultivated peat soils. Geoderma 154, 173–180. Christensen, B.T., 1987. Decomposability of organic matter in particle size fractions from field soils with straw incorporation. Soil Biol. Biochem. 19, 429–435. De Brogniez, D., Ballabio, C., Stevens, A., Jones, R.J.A., Montanarella, L., van Wesemael, B., 2015. A map of the topsoil organic carbon content of Europe generated by a generalized additive model. Eur. J. Soil Sci. 66, 121–134. Ekström, G., 1927. Klassifikation av Svenska Åkerjordar (Classification of Swedish Arable Soils) Ser C. No. 345 (Årsbok 20). Sveriges Geologiska Undersökning, Uppsala (in Swedish). (161 pp). Eriksson, J., Mattsson, L., Söderström, M., 2010. Tillståndet i svensk åkermark och gröda. Naturvårdsverket Rapport 6349 (in Swedish) 131 pp. Available online: https://www. naturvardsverket.se/Documents/publikationer/978-91-620-6349-8.pdf (verified: October 26, 2017). FAO, 2017. GSP Soil Data Policy. Available online: http://www.fao.org/3/a-bs975e.pdf (verified: October 26, 2017). Gee, G.W., Bauder, J.W., 1986. Particle-size analysis. In: Klute, A. (Ed.), Physical and Mineralogical Methods. Soil Science Society of America, Madison. Gustafsson, K., 1999. Models for precision application of lime. In: Stafford, J. (Ed.), Precision Agriculture ‘99 Odense. Denmark, SCI.
We describe production of the digital soil map of Sweden (DSMS), a map of topsoil texture and SOM content in arable land created through the use of data mining techniques whereby a set of detailed covariates were combined with a large number of reference samples. We show that transparency in the validation procedure is crucial if users are to understand the accuracy of a map product. Access to a large validation dataset consisting of soil analyses at farm level throughout southern Sweden made it possible to demonstrate the performance of DSMS at multiple scales. However, although excellent validation statistics were obtained at the overall scale, these figures reveal very little on the performance on farms or within fields. The clay content map produced was more accurate than the sand content map. The SOM content could not be calibrated with acceptable accuracy for precision agriculture from the available predictor data, but DSMS SOM has been used to produce 1 km × 1 km grid data for the FAO GSOC map. Through the validation procedure, we showed that DSMS is a versatile product that can be applied at different scales and for different objectives. For example, it can be used within farms as decision support for site-specific management (i.e. precision agriculture) and it can be used at municipality, county and national level as input to environmental impact scenarios and modelling of soil processes. Further developments 8
Geoderma xxx (xxxx) xxx–xxx
K. Piikki, M. Söderström
artiklar/ra1519.html (verified: October 26, 2017). Piikki, K., Söderström, M., Stenberg, B., 2013. Sensor data fusion for topsoil clay mapping. Geoderma 199, 106–116. R Core Team, 2017. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria Available online: https:// www.R-project.org (verified: October 26, 2017). Rawls, W.J., Pachepsky, Y.A., Ritchie, J.C., Sobecki, T.M., Bloodworth, H., 2003. Effect of organic carbon on soil water retention. Geoderma 116, 61–76. Söderström, M., Eriksson, J., 2013. Gamma-ray spectrometry and geological maps as tools for cadmium risk assessment in arable soils. Geoderma 192, 323–334. Söderström, M., Sohlenius, G., Rodhe, L., Piikki, K., 2016. Adaptation of regional digital soil mapping for precision agriculture. Precis. Agric. 17, 588–607. Söderström, M., Piikki, K., Cordingley, J., 2017. Improved usefulness of continental soil databases for agricultural management through local adaptation. S. Afr. J. Plant Soil 34, 35–45. Stoorvogel, J.J., Bakkenes, M., Temme, A.J., Batjes, N.H., Brink, B.J., 2017. S-world: a global soil map for environmental modelling. Land Degrad. Dev. 28, 22–33. Van der Klooster, E., Van Egmond, F.M., Sonneveld, M.P.W., 2011. Mapping soil clay contents in Dutch marine districts using gamma-ray spectrometry. Eur. J. Soil Sci. 62, 743–753. Van Veen, J.A., Kuikman, P.J., 1990. Soil structural aspects of decomposition of organic matter by micro-organisms. Biogeochemistry 11, 213–233. Weiss, A., 2001. Topographic Position and Landforms Analysis. Poster presentation, ESRI User Conference, San Diego, CA. Available at: http://www.jennessent.com/ downloads/tpi_documentation_online.pdf (Verified: October 10, 2017).
Hastie, T., Tibshirani, R., Friedman, J., 2009. The Elements of Statistical Learning: Data Mining Inference and Prediction, 2nd edition. Springer Series in Statistics, New York, pp. 383–411. Hengl, T., de Jesus, J.M., Heuvelink, G.B., Gonzalez, M.R., Kilibarda, M., Blagotić, A., et al., 2017. SoilGrids250m: global gridded soil information based on machine learning. PLoS One 12, e0169748. IAEA, 2003. Guidelines for Radioelement Mapping Using Gamma Ray Spectrometry Data. International Atomic Energy Agency, Vienna Austria IAEA-TECDOC-1363. Available online. http://www-pub.iaea.org/mtcd/publications/pdf/te_1363_web.pdf (verified: October 26, 2017). Janssen, P., Heuberger, P., 1995. Calibration of process-oriented models. Ecol. Model. 83, 55–66. Kätterer, T., Andrén, O., Jansson, P.-E., 2006. Pedotransfer functions for estimating plant available water and bulk density in Swedish agricultural soils. Acta Agr. Scand. B-S P 56, 263–276. Malone, B.P., McBratney, A.B., Minasny, B., 2013. Spatial scaling for digital soil mapping. Soil Sci. Soc. Am. J. 77, 890–902. McBratney, A.B., Santos, M.L.M., Minasny, B., 2003. On digital soil mapping. Geoderma 117, 3–52. Milborrow, M., 2017. Package ‘Earth’. Available online: https://cran.r-project.org/web/ packages/earth/earth.pdf (verified: October 26, 2017). Nash, J.E., Sutcliffe, J.V., 1970. River flow forecasting through conceptual models part I -a discussion of principles. J. Hydrol. 10, 282–290. Paulsson, R., Djodjic, D., Carlsson, Ross C., Hjerpe, K., 2015. Nationell jordartskartering Matjordens egenskaper i åkermarken. Swedish Board of Agriculture Report: RA15:19 (in Swedish). 44pp. Available online: http://webbutiken.jordbruksverket.se/sv/
9
Contents lists available at ScienceDirect
Geoderma journal homepage: www.elsevier.com/locate/geoderma
Digital soil mapping of arable land in Sweden – Validation of performance at multiple scales Kristin Piikki⁎, Mats Söderström Department of Soil and Environment, Swedish University of Agricultural Sciences (SLU), P.O. Box 234, SE-53223 Skara, Sweden
A R T I C L E I N F O
A B S T R A C T
Handling Editor: A.B. McBratney
In this study, we produced a detailed digital soil map of topsoil texture and soil organic matter (SOM) content for 2.4 million ha of arable land in Sweden (DSMS). Three spatially exhaustive datasets (a laser-scanned digital elevation model, airborne gamma radiation scanning data and a legacy Quaternary deposit map) were calibrated against topsoil texture and SOM content in around 13,500 soil samples, using multivariate adaptive regression splines (MARSplines) modelling. We then deployed the MARSplines models to produce raster maps (50 m × 50 m) of clay, sand and SOM content. The modelling procedure was validated by an independent dataset of about 24,000 samples clustered on 544 farms (with a local sample density of one per 3 ha). The error in clay content was < 8% in 75% of the validation samples, while for sand content and SOM content it was 13% and 2%, respectively. Corresponding values for the farm-average level were 6%, 11% and 2%, respectively. Modelling efficiency values revealed that the clay content map was a considerable improvement over the mean of the reference values at national level, regional level and, in many cases, also farm level. However, SOM content predictions showed little or no improvement over the mean of the reference samples (at any scale), due to poor correlation with the exhaustive predictor variables at all three scales investigated. The DSMS soil geodatabase will continue to be improved and have more layers added, e.g. derived layers calculated from the primary clay, sand and SOM layers by use of pedotransfer functions. Practical use of DSMS is exemplified here in an internet application for deriving prescription files for precision agriculture.
Keywords: Texture Organic matter Digital elevation model Gamma radiation MARSplines Quaternary deposit
1. Introduction Ever since McBratney et al. (2003) outlined an explicit framework for digital soil mapping (DSM), there has been a virtual explosion in the availability of spatial data and pedometricians world-wide creatively applying new computational techniques in order to produce detailed soil maps (e.g. Arrouays et al., 2014; De Brogniez et al., 2015; Hengl et al., 2017; Stoorvogel et al., 2017). DSM involves determining soil variations in relation to the landscape, finding measurable covariates for the soil property of interest and developing quantitative models for prediction of the target property. However, while a number of DSM products provide sufficiently high spatial resolution to give the impression that they are applicable for management decisions at farm level, validations with local data often show very large errors, indicating that the data are not suitable for use in precision agriculture (Söderström et al., 2016, 2017). An effort was recently launched in Sweden to produce a freely available digital map of topsoil texture (Digital Soil Map of Sweden, DSMS) for most of the country's arable land (about 90% or ~2.4 million ha) at an accuracy sufficient for use in
⁎
precision agriculture. This might be possible, since three meticulous national datasets describing the landscape are available and could serve as covariates in modelling. These are: a digital elevation model (DEM) based on airborne laser scanning, a national airborne gamma radiation survey and a Quaternary deposit (QD) map. There is also a newly available reference dataset of soil samples (national 1-km grid) funded by the government in order to improve estimates of nutrient loads to the Baltic Sea (Paulsson et al., 2015). In addition, a large dataset of geographically positioned soil analyses on 544 farms is available for validations at farm scale. A recent pilot study covering about 100,000 ha of intensively cultivated land in south-west Sweden reported root mean square error of prediction of 5.7% for clay, 12.5% for sand and 1.3% for soil organic matter (SOM) (Söderström et al., 2016). Compared with maps of soil texture produced from farmers' own data using contemporary techniques (interpolated maps based on one sample per 3 ha), the regional DSM produced in the pilot study was equally good or better. It was based on the DEM, gamma radiation survey, QD map and reference soil samples mentioned above. The model type used was multivariate adaptive regression splines
Corresponding author. E-mail addresses: [email protected] (K. Piikki), [email protected] (M. Söderström).
http://dx.doi.org/10.1016/j.geoderma.2017.10.049 Received 26 April 2017; Accepted 27 October 2017 0016-7061/ © 2017 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/BY/4.0/).
Please cite this article as: Piikki, K., Geoderma (2017), http://dx.doi.org/10.1016/j.geoderma.2017.10.049
Geoderma xxx (xxxx) xxx–xxx
K. Piikki, M. Söderström
of within 5 m (various brands). Data on these soil samples are spatially complementary to, and were merged with, data obtained in a previous, spatially sparser, soil sampling campaign undertaken during 2001–2007 by the Swedish Environmental Protection Agency (SNV; Naturvårdsverket, Stockholm, Sweden) (Eriksson et al., 2010). Only soil samples with SOM ≤20% were included in the present modelling work. In total, 13,485 soil samples constituted the grid dataset used for calibration of DSM prediction models (11,734 SJV samples and 1751 SNV samples). For local-scale validation, we used an entirely independent dataset collected on farms located mostly in intensely cultivated districts of Sweden, within an environmental support subsidies programme funded through SJV (years 2007–2009). These samples were distributed within farms in a grid pattern, with a density of one sample per three hectares. In this dataset, sample positions were available but farm identities were not. Therefore, synthetic farms were generated through a spatial clustering procedure. Samples within 1000 m of each other were grouped into ‘farms’ and only those with at least five samples were kept for validation. This resulted in a validation set of 24,323 analyses on 544 farms. The median size of these farms was 126 ha. All soil samples in the grid dataset and the farm dataset were taken from the topsoil (0–20 cm; 7–10 topsoil cores within a 3–5 m radius circle at each sampling point were bulked). The samples were air-dried at 35–40 °C, milled and sieved through a 2-mm mesh. The content of sand (0.2–2 mm) was determined by sieving and weighing and the content of clay (< 2 μm) by the sedimentation method (method ISO 11277; Gee and Bauder, 1986). The sand and clay content were both expressed as a percentage of the fine soil fraction (< 2 mm). The SOM content was determined from the clay content and the loss on ignition (Eq. (1); original method by Ekström, 1927). The laboratory analyses were carried out by Eurofins laboratory (Kristianstad, Sweden).
(MARSplines). Appropriate validation of DSM is a challenge. The conventional scale concept used for printed maps is not directly applicable for digital soil maps. Malone et al. (2013) list i) extent, ii) resolution and iii) support as appropriate scale specifications of digital soil maps in raster format. In the present study, validations were carried out at two support (soil samples and farms) and three extents (farm extent, regional extent and national extent). The aims of this study were to: 1) Describe the methods used for production of DSMS, which shows topsoil texture and SOM for arable land. 2) Present the results of validations of DSMS at multiple scales, focusing on assessing local, independent within-farm map performance. 3) Suggest practical user applications based on DSMS. Aim (2) deviates from that for most DSM products, which are generally validated through various types of cross-validation using subsets of the soil mapping dataset. 2. Materials and methods 2.1. Area of interest The mapped area comprises approximately the southern third of Sweden, which includes about 90% of the agricultural land in the country (Fig. 1b). A spatial database of cultivated arable land derived from the EU subsidies database of the Swedish Board of Agriculture (SJV; Jordbruksverket, Jönköping, Sweden), the so-called ‘Block Map’ (from 2013), was used to delineate the area of interest map in detail (example shown in Fig. 1a). Areas of organic soils identified from the 1:50,000 QD map by the Geological Survey of Sweden (SGU, Uppsala, Sweden) and of wetlands identified from the Block Map 2013 were classified as organic, and were excluded from the modelling. The reason for excluding organic soils from the map is that soil samples with SOM content > 20% are only rarely analysed for soil texture, so there were insufficient reference data for mapping areas with this type of soil.
SOM = Loss on ignition − 0.46 − 0.047 × clay content
(1)
2.3. Covariates Three types of spatially exhaustive ancillary information sources were used in the predictive modelling: airborne gamma-ray spectrometry data, elevation data acquired through airborne laser-scanning and information from legacy QD maps. All predictor data were calculated for (or simply extracted to) the centre point of a 50 m × 50 m raster mesh. The pre-processing of covariate data followed the procedure adopted in the pilot study to this DSM work (Söderström et al., 2016). In order to handle the very large amount of data, the map area
2.2. Soil samples and soil analysis methods During 2011–2012, SJV implemented a national soil sampling campaign, resulting in a more or less regular grid of one soil sample per km2 across Swedish agricultural land. Positions were determined by global positioning system (GPS) with an estimated positional accuracy
Fig. 1. Spatial distribution of the soil datasets used. Map a) shows the spatial distribution of the farm and grid datasets over the agricultural fields (blocks) in a small area and map b) shows the entire extent of the digital soil mapping area in southern Sweden and the distribution of agricultural land (2.4 million ha) and lakes.
2
Geoderma xxx (xxxx) xxx–xxx
K. Piikki, M. Söderström
recorded concentrations of K < 0% or Th < 1 ppm were omitted, since these often indicate data affected by water bodies or peat. A small number of abnormally high values of K and Th were also present in the dataset and are detrimental for construction of robust calibration models, even if they are correct. In many cases it is not variations in soil texture that cause such extremes, but rather anomalies in the geological conditions (e.g. certain bedrock outcrops). Therefore, K concentrations > 4.5% were set to 4.5%, and Th concentrations > 30 ppm were set to 30 ppm. In some rare cases, abrupt discrepancies in values were noted between flight lines, but in most cases it was possible to adjust the general level between such measurements manually. If not, such data were deleted. The filtered K and Th data were then interpolated to the 50 m × 50 m mesh by ordinary block kriging. 2.3.2. Elevation data The national digital elevation model derived from airborne laser scanning and supplied by LMV as a raster with 2-m spatial resolution (Grid 2+, LMV, Gävle, Sweden) was used as a data source for elevation and landform. The data were resampled (averaged) to 10 m × 10 m, which was considered more appropriate for the present purpose. Various measures of landform were estimated as topographic position index (TPI; Weiss, 2001), i.e. the relative elevation of raster cells compared with circular neighbourhoods of 5 ha, 50 ha and 500 ha (denoted TPI5, TPI50 and TPI500, respectively) (Fig. 3).
Fig. 2. Sub-regions used for modelling.
2.3.3. Legacy soil data on Quaternary deposits The Weichselian glaciation, the melting of the inland icecap and the subsequent sea-level fluctuations around the Scandinavian Peninsula in the Holocene have determined the overall distribution of the regolith in Sweden. The genetics of these young sediments in Sweden are described by geologists in QD maps, scale 1:50,000, which are available as vector polygon maps from SGU. More than 100 different deposit types have been assigned. To transfer these data into useable covariates, in this work we simplified the deposit classification into two different versions:
was manually divided into 22 sub-regions based on expert knowledge (Fig. 2), in a process that was a trade-off in terms of uniformity of geomorphology and geology and amount of soil observations. This process followed an initial test in which data from six different parts of the study area with different physical characteristics (for example one hilly region with granites known for their high level of radiation, another area mostly consisting of fine post-glacial sediments in a relatively flat landscape) were analysed separately and also jointly. The separate models performed somewhat better, indicating that such an approach was preferable. In some cases, when the physical characteristics were less distinct and the soil observations otherwise were judged to be too few or spatially skewed in a sub-region, or to reduce the effect of slight differences in predicted values between two adjacent regions, the sub-regions overlapped to some extent.
• Seven classes (QD7): Clay, Sand, Silt, Till, Till clay, Organic and Bedrock • Three classes (QD3): Clay, Organic and Other The ‘Organic’ class, which was the same in the QD7 and the QD3 classifications, was used to mask organic areas from the map extent and no samples from areas classified as organic were included in the model parameterisation. We also simplified the geometry of the QD map, by dissolving polygons smaller than 1 ha with the polygon with which they shared their longest border. The reason for this is that the size of such small features is often exaggerated for cartographic reasons in QD maps and this, in combination with positional uncertainty about the extent, makes them unsuitable to include in modelling.
2.3.1. Gamma-ray spectrometry data Natural emissions of gamma radiation in Sweden have been measured by airborne gamma-ray spectrometry since the 1960s by SGU. Data are collected along flight lines normally separated by 200 m, from an altitude of 30 m or 60 m, by recording the concentration of the naturally occurring isotopes 232Th (ppm; eTh), 40K (ppm) and 238U (%; eU) every 16–17 m along the flight lines. The data are then normalised by SGU to take into account variations in atmospheric radiation, topographical influence and variations in altitude (SGU, personal communication). The footprint radius of recorded values on the ground is approximately four times the altitude, and about 80% of the response originates from the upper 0.3 m of soil (IAEA, 2003). The usefulness of gamma-ray data in predictive mapping of soil properties in northern Europe has been demonstrated in several earlier studies, e.g. by Piikki et al. (2013), Söderström et al. (2016) and Van der Klooster et al. (2011). The eU shows considerable short-range variation and smallscale noise in parts of the dataset (Söderström and Eriksson, 2013; Söderström et al., 2016) and was omitted from the calibrations in the present analysis to reduce the risk of erroneous and over-fitted prediction models. Further filtering of the gamma-ray data included removal of observations < 100 m from water bodies (as delimited by the 1:50,000 topographical map of the Land Survey of Sweden (LMV; Lantmäteriet, Gävle, Sweden)) and > 300 m from arable land. To remove outliers and data considered non-representative for arable land,
2.4. Spatial predictive modelling MARSplines, which is a data mining method used to parameterise statistical relationships between the variables in a dataset (X; predictors Y; response variables), was used for modelling. A detailed description of MARSplines models is given by Hastie et al. (2009). In brief, such a model consists of a set of piece-wise, univariate linear regression functions (so-called basis functions) that are defined above or below threshold values (breakpoints) of the predictor in question. MARSplines models can be additive or allow for interactions among the basic functions. We chose to use only additive models, in order to minimise the risk of overfitting associated with more complex regression models. Before the model parameterisation, factor expansion was performed on the categorical data (QD3 and QD7). Each class was transformed to one binary variable, where 1 indicated that the sample or raster cell belonged to the QD class in question, while 0 indicated the opposite. In 3
Geoderma xxx (xxxx) xxx–xxx
K. Piikki, M. Söderström
Fig. 3. Covariates used in modelling. a) Elevation z above sea level; b–d) relative topography TPI500, TPI50 and TPI5, i.e. the elevation relative to the average of a neighbourhood area of 500 ha, 50 ha and 5 ha, respectively; e) concentration of 232Th; f) concentration of 40K; g) Quaternary deposits simplified into seven classes (QD7); and h) Quaternary deposits simplified into three classes (QD3). In e–h), data are shown only for areas of arable land, with elevation as background. Blue areas indicate lakes in all maps. The overview map show the extent of the digital soil map of Sweden and the filled circle indicate the location of the 10 km × 10 km area in maps a–h. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
farm dataset. The following statistics were calculated: mean absolute error (MAE; Janssen and Heuberger, 1995) and Nash-Sutcliffe model efficiency (E; Nash and Sutcliffe, 1970). The latter indicates how well a plot of observed values and values predicted by a model fits the 1:1line, and is a useful index for testing whether the map values are more accurate than the average of the reference samples. The value of E varies between −∞ and 1, where a negative value indicates that the map is less accurate than the average of the reference samples, E = 0 means that the map is just as accurate as the reference sample average, a positive value of E means that the map is an improvement over the average and a value of 1 indicates that the map values at the soil sampling locations are equal to the soil analysis results, i.e. perfect correlation.
order to construct as robust models as possible, we also used a built-in cropping function to remove all basis functions that did not significantly contribute to the prediction. In summary, the following eight covariate predictors (Fig. 3) were used to derive models for clay, sand and SOM content:
• Concentrations of naturally occurring radioactive isotopes: and K • Elevation and relative topography: z, TPI , TPI and TPI • Quaternary deposits: QD3 and QD7
232
Th
40
5
50
500
Models were calibrated by the reference data (the grid soil dataset combined with covariate predictors) in each sub-region (Fig. 2) and then deployed onto the 50 m × 50 m mesh for that area. The predictions in all 22 sub-regions were then merged to a final raster dataset covering the entire map area. In some cases, neighbouring sub-regions overlapped. In the overlapping areas the final prediction was a distance-weighted average of the two predicted values. The distance was calculated from the raster cell in question to the edge of the overlapping area. That yielded a smooth transition from one model area to the next.
2.7. Software All spatial analyses were conducted in ArcGIS with the extensions Spatial Analyst and Geostatistical Analyst (ESRI, Redlands, CA, USA), while the MARSplines modelling was performed using R (R Core Team, 2017), package Earth (Milborrow, 2017). Statistica 13 (Dell Inc., Tulsa, OK, USA) was also used for data management and analyses.
2.5. Multiple-scale validation 3. Results Validation statistics were calculated for all orthogonal combinations of the two support levels and the three spatial extents, except at the farm extent, for which map accuracy could not be tested at farm support (since n = 1 in each validation set). This means that at the national extent one value of E and one value of MAE was calculated. In the regional case, 22 values of MAE and E were calculated and subsequently averaged across all modelling areas, while at farm level 544 values of MAE and E were calculated and then averaged across all farms.
3.1. Descriptive statistics The quartiles of the frequency distributions were similar for the calibration data (grid dataset) and the validation data (farm dataset) for all three map variables (Table 1.) As can be inferred from the interquartile ranges, the sand content in the laboratory-analysed samples was right-skewed, while the clay content was somewhat left-skewed and the frequency distribution of SOM had an extreme skewness to the left. At this national level, the frequency distributions of the predicted soil properties were very similar to the frequency distributions of their laboratory-analysed counterparts. However, there were deviations at the right tail of the distributions (especially for the SOM content), indicating that under-prediction of high values was more common than
2.6. Validation statistics Data from the prediction rasters (clay content, sand content and SOM content) were extracted for all validation sample locations in the 4
Geoderma xxx (xxxx) xxx–xxx
K. Piikki, M. Söderström
Table 1 Descriptive statistics on the response variables: values for individual analyses of the grid dataset (n = 13,485), used for calibration, and for the farm dataset, used for validation. For the latter, statistics for the individual sample values (n = 24,323) are given before the slash and for farm-average values (n = 544) after the slash. SOM = soil organic matter. Dataset
Value
Variable
Minimum
1st Quartile
Median
3rd Quartile
Maximum
Grid Grid Grid Farm Farm Farm Farm Farm Farm
Measured Measured Measured Measured Measured Measured Predicted Predicted Predicted
Clay content [%] Sand content [%] SOM content [%] Clay content [%] Sand content [%] SOM content [%] Clay content [%] Sand content [%] SOM content [%]
2 0 0 1/1 1/4 0/1 0/3 0/4 0/1
10 14 3 12/12 14/20 2/3 13/13 18/22 1/1
19 38 4 20/18 37/43 3/3 21/18 37/44 2/2
33 59 5 36/28 58/56 4/4 35/28 54/55 2/2
80 96 20 73/58 96/91 19/17 63/56 93/81 7/4
Table 2 Cross-tables of Pearson correlation coefficient (r) among all pairs of continuous variables. The r values were calculated and averaged at: a) national level; b) sub-region level and c) farm level. 232
Th
40
TPI50
TPI500
Clay
Sand
SOM
– 0.69 −0.27 0.26 −0.20
– − 0.23 0.27 − 0.18
– − 0.85 0.16
– − 0.22
–
b) Average of the pair-wise correlation coefficients calculated among all samples in each of the model areas 232 Th – 40 K 0.32 – z − 0.18 −0.05 – RT5 0.02 0.08 0.02 – RT50 − 0.05 0.12 0.10 0.62 – RT500 − 0.13 0.09 0.26 0.22 0.57 Clay 0.60 0.12 − 0.26 − 0.05 −0.13 Sand − 0.53 −0.12 0.23 0.00 0.09 SOM 0.06 −0.25 0.09 − 0.15 −0.17
– − 0.21 0.22 − 0.14
– − 0.79 0.20
– − 0.24
–
c) Average of the pair-wise correlation coefficients calculated among all samples on each of the farms Th – K 0.33 – z − 0.12 0.11 – RT5 0.05 0.09 0.17 – RT50 − 0.01 0.16 0.39 0.61 – RT500 − 0.08 0.14 0.66 0.29 0.66 Clay 0.35 0.03 − 0.19 − 0.09 −0.17 Sand − 0.35 −0.01 0.21 0.05 0.15 SOM 0.03 −0.24 − 0.16 − 0.18 −0.22
– − 0.18 0.21 − 0.19
– − 0.76 0.23
– − 0.30
–
a) Pair-wise correlation coefficients 232 Th – 40 K 0.41 z − 0.26 TPI5 − 0.16 TPI50 − 0.25 TPI500 − 0.25 Clay 0.83 Sand − 0.76 SOM 0.08
K
z
TPI5
calculated among all samples in the national dataset – 0.01 0.00 0.02 0.00 0.27 −0.26 −0.18
– 0.08 0.15 0.19 − 0.24 0.13 0.11
– 0.65 0.31 − 0.19 0.17 − 0.16
232 40
z = elevation above sea level, TPI = topographic position index.
soils than in soils with a high clay content because the clay particles protect the SOM from microbial action (e.g. Christensen, 1987; van Veen and Kuikman, 1990). When quantifying and interpreting the correlations between SOM and the predictor variables, one should bear in mind that SOM values > 20% were excluded from the dataset, which means that the correlations are evaluated for mineral soils only and not for organic soils.
over-prediction. Correlation matrices of Pearson correlation coefficients (r) were calculated for three different spatial extents (Table 2). It is possible that different processes govern covariance between two variables at different scales. A weaker correlation (smaller absolute value of r) at finer scale may also be due to the range of variation and the number of observations being smaller when a smaller area is considered. Therefore, two variables can be well correlated at the national scale, but show a weaker correlation for regional model areas or farms. In the present case, for example, there was much lower correlation between Th and clay content at farm scale (Table 2c) compared with the entire area (Table 2a) or compared with the sub-region level (Table 2b), probably due to the fact that farms with small standard deviation in Th and clay content reduce the average r value at farm scale. The SOM content showed weak correlations to all predictors used in this study, but was most strongly related to the content of sand. This correlation had a negative sign at all three extents. This may actually be a causal relationship, since mineralisation of organic matter is faster in sandy
3.2. Absolute map errors The graphs in Fig. 4 visualise the MAE distributions in the predicted soil properties. As expected, the higher percentiles were smaller for the farm-wise aggregated data (solid lines) than for the individual samples (dotted lines). The error in clay content was < 8% in 75% of the validation samples, while for sand content it was 13% and for SOM content it was 2%. Corresponding values for farm averages were 6%, 11% and 2%, respectively. This kind of validation at two support levels is important for assessment of possible practical applications of the DSMS 5
Geoderma xxx (xxxx) xxx–xxx
K. Piikki, M. Söderström
Fig. 4. Absolute error. Distribution of mean absolute error (MAE) of predicted soil properties: a) clay content, b) sand content and c) soil organic matter (SOM) content for individual soil samples (n = 24,323) and farm averages of these soil samples (n = 544). Relative frequency on the y-axis is the fraction of samples within discrete intervals of the x-axis (width of 5%units for clay and sand and 1%-unit for SOM). The integrated area under the each curve is equal to 1. Fig. 5. Maps of predictions and prediction errors. Predicted a) clay content and b) sand content of the digital soil map of Sweden (DSMS) and (c, d) geographical distribution of the magnitude of error calculated for 25 km × 25 km cell areas. As 70% of the area of Sweden is forested, the coverage appears sparse in the map, even though most of the arable land (~90%) is mapped. MAE = mean absolute error.
land that was not mapped in the current version of DSMS (one small 10 km × 1 km area) is located there. The reason for this was missing airborne gamma data. In the agricultural districts in the south and south-west of the country, the clay content prediction errors were relatively small (< 5%). The smallest clay prediction errors (in general about 2.5%) were found in the central part of southern Sweden. In this region, both the clay content and the magnitude of its variation are low. For sand content (Fig. 5d), the variation did not follow a similarly clear spatial pattern. In the districts with large areas of contiguous arable
data. For precision agriculture, point location validation is necessary, whereas e.g. for regional modelling of nutrient leaching, farm-level aggregated data may be adequate. Fig. 5 shows the geographical distribution of MAE for the clay content and sand content. Large clay content errors coincided generally with areas of high clay content in the topsoil, e.g. in the east and north-east parts of the map area (Fig. 5c). Particularly high errors for clay were recorded in the east-central region, where the airborne gamma radiation scanning also showed some quality problems in places. In fact, the only area of arable (non-organic)
6
Geoderma xxx (xxxx) xxx–xxx
K. Piikki, M. Söderström
support). In all cases, the E value was higher for farm average values than for individual samples. However, as Fig. 4 shows, the MAE for individual samples did not differ much from that of farm-averaged values, although the error was more heterogeneous and showed a larger span at sample support level than at farm support level. Based on the MAE of SOM content (Fig. 4c), the DSMS SOM layer could be deemed acceptable for regional applications, but the scatter plots (Fig. 6) and the E values (Table 3) indicate a weak correlation between the map and the reference samples, indicating that an average of all samples in a region may be almost as accurate as the DSMS SOM layer. The reason that SOM content could not be mapped well enough to be considered very useful at the local scale is most likely that it was poorly correlated to the exhaustive predictors (Table 2). Even if advanced data mining methods are used, it is not possible to parameterise relationships that are not present in the dataset (causal or not). Plausible reasons why no useful relationship could be parameterised between the available predictors and SOM content are that the small-scale (between-field) variation in SOM is highly influenced by management, that the local magnitude of SOM variation is often small, that the laboratory analysis uncertainty is relatively high in relation to the variation, and that topsoil gamma radiation is not very much affected by the SOM content unless the SOM content is high. As an example of the latter, Berglund and Berglund (2010) described a method for delineating peat areas using the activity concentration of 40K as a predictor. Since all samples for calibration and validation in the present study were obtained during a 12-year period, it is possible that the SOM content had changed to some extent, at least in fields where management had changed, adding further difficulty in predicting variation. At national and sub-region level, DSMS was an improvement over using the mean of the reference samples (both farm-average values and point location data, SOM content excluded; Table 3). For individual farms the average E values were closer to 0, which indicates that DSMS is not always an improvement over the soil sample mean. However, there was considerable standard deviation in the E values, leading to the interpretation that for many farms, but not all, DSMS is an improvement over the reference sample mean. DSMS data can be used as input data in national, overall modelling of phosphorus (P) leaching from agricultural soils, which is one of the reasons why donors funded this project. It can also be used to target control measures against P leaching to parts of fields where they would be most effective. This would be a considerable improvement for environmental authorities, which otherwise can often only decide on very general mitigation actions for large areas to control e.g. leaching of nutrients. Such measures are likely to be more efficient if they can be targeted at high-risk areas and hotspots. For example, structure liming is recommended to be applied to soils with > 15% clay (Albertsson et al., 2016). DSMS data can also be used as the basis for generating variable-rate application (VRA) files aimed at practical precision agriculture. In an interactive internet application (www.markdata.se; a free-to-use web application in Swedish language), DSMS clay content is classified by the user, and outputs such as VRA files of seed rate can be produced and downloaded. Most users choose to reduce the seed density in areas with a coarser texture, where germination is often better, in order to achieve an optimal canopy in all parts of the field. Products such as DSMS can also be expanded to include derivative products through pedotransfer functions. Bulk density, water content at field capacity and at wilting point and porosity are examples of properties that can be calculated from DSMS data (Rawls et al., 2003; Kätterer et al., 2006). In 2017, the DSMS SOM layer was used as input to the Swedish contribution to the Global Soil Organic Carbon (GSOC) map (FAO, 2017). In that case, pedotransfer functions were used to estimate bulk density from SOM (Eriksson et al., 2010), and the SOC stock data were aggregated to a 1 km × 1 km raster. Another application of DSMS involving pedotransfer functions is in development of map layers useful for determining the amount of lime required for increasing soil pH in acid soils. According to guidelines from the Swedish Board of
Table 3 Modelling efficiency (E) of texture and soil organic matter (SOM) predictions calculated for point locations and points aggregated to farm averages at three scales (farm, subregion and entire digital soil map of Sweden (DSMS) area) (mean ± standard deviation among farms or modelling areas).
Clay Sand SOM
Farms
Model areas
Entire DSMS area
Point locations
Point locations
Farm averages
Point locations
Farm averages
0.18 ± 0.22 0.18 ± 0.21 0.05 ± 0.07
0.40 ± 0.22 0.37 ± 0.23 0.07 ± 0.06
0.60 ± 0.26 0.54 ± 0.28 0.07 ± 0.07
0.76 0.75 0.12
0.90 0.86 0.14
land, prediction errors of sand content were moderate (10–13%), whereas the largest errors (up to 15% and above) were found in areas where agricultural land is interspersed with forested areas and fields in general are smaller. 3.3. Relative map accuracy An alternative to using regional digital soil maps is to use averages of samples at the farm level or possibly averages for larger regions. Table 3 shows the effects of this on the performance of texture and SOM modelling of aggregation of predictions made at a finer scale. In general, models for clay content performed better than those for sand content. For SOM, the models performed less well. It is clear that aggregation resulted in much more accurate predictions, e.g. the E value for clay content of farms in the entire validation dataset (E = 0.90) indicates a very accurate and precise map. However, looking at finer scale, it is obvious that the predictions were not always successful in describing the variation within the mapped area. The E values consistently increased with increasing extent and also with increasing support. 4. Discussion It is a challenge to present relative validation statistics for a largeextent, high-resolution dataset like DSMS in a manner that is useful for different possible applications of the data (in terms of resolution and extent). This is important to consider when reporting accuracies from validations, and also for users who wish to understand how well a prediction model works, i.e. how accurate a digital soil map is. Correlation analyses at one extent are not directly transferrable to another (Table 2). The same is true for E (Table 3; Fig. 6), which was always smaller at farm level than for the entire modelling area, probably because of the generally much smaller ranges of variation in soil properties and a smaller number of observations in smaller areas (Söderström et al., 2016). Moreover, testing for different extents might indicate that causal and/or indirect relationships between the predictors and the response variables may be governed by different processes operating at different scales. Hence, the digital soil maps produced may accurately describe this at one spatial level, but not at another. For example, the coarse variation at the national scale (e.g. governed by large scale geological variation) may be well described by DSMS data, but the data may fail to describe the small-scale variation within farms or sub-regions (e.g. governed by agricultural management). The rationale behind using two aggregations of the validation soil samples (individual samples and farm average values) was that some covariate data used in the prediction models were coarse. For example, the footprint of the gamma radiation data was much larger than the support of an individual soil sample. Individual samples will therefore not necessarily be the most appropriate resolution for applications. Here we tested whether the accuracy of the predictions was better for farm averages than for describing within-field variation (sample 7
Geoderma xxx (xxxx) xxx–xxx
K. Piikki, M. Söderström
Fig. 6. Scatterplots of predictions against laboratory values. The plots show individual soil samples (a–c; n = 24,323; a random subset of 2000 samples are plotted) and farm averages of these soil samples (d–f; n = 544). SOM = soil organic matter.
include local adaptation through interactive functionality where users can upload their own data, and the application of pedotransfer functions to produce important map layers of soil functional properties.
Agriculture, the models applied in Sweden for determining both the recommended suitable soil pH (target-pH) and the buffering capacity are a function of soil clay content and SOM (Gustafsson, 1999). A test of producing such map layers in DSMS was reported by Söderström et al. (2016) for use at the field level and a method was presented to improve the accuracy of digital soil maps locally through addition of locally available data. Such a procedure was recently further developed and assessed (Piikki et al., 2013), and an automated interactive method for users to upload their own data in a web application to remodel the DSMS at farm scale was presented. With such a system in place, conventional soil mapping (with one sample per 3 ha) was better than the digital soil map for less than one-third of the validated farms.
Acknowledgements This work was financially supported by the Swedish National Space Board (SNSB) [dnr 214/13] and the Geological Survey of Sweden (SGU) [project agreement: 15-03-31]. We acknowledge valuable contributions from the Swedish Board of Agriculture (JV), The Rural Economy and Agricultural Societies (HS), Precision Agriculture Sweden (POS) and many interested farmers. References
5. Conclusions
Albertsson, B., Börling, K., Kvarmo, P., Listh, U., Malgeryd, J., Stenberg, M., 2016. Rekommendationer för gödsling och kalkning 2017. Report JO16 24 Swedish Board of Agriculture, Jönköping, Sweden 102 pp. Available online: http://webbutiken. jordbruksverket.se/sv/artiklar/jo1624.html (verified: October 26, 2017). Arrouays, D., Grundy, M.G., Hartemink, A.E., Hempel, J.W., Heuvelink, G.B., Hong, S.Y., et al., 2014. Chapter three - GlobalSoilMap: toward a fine-resolution global grid of soil properties. Adv. Agron. 125, 93–134. Berglund, Ö., Berglund, K., 2010. Distribution and cultivation intensity of agricultural peat and gyttja soils in Sweden and estimation of greenhouse gas emissions from cultivated peat soils. Geoderma 154, 173–180. Christensen, B.T., 1987. Decomposability of organic matter in particle size fractions from field soils with straw incorporation. Soil Biol. Biochem. 19, 429–435. De Brogniez, D., Ballabio, C., Stevens, A., Jones, R.J.A., Montanarella, L., van Wesemael, B., 2015. A map of the topsoil organic carbon content of Europe generated by a generalized additive model. Eur. J. Soil Sci. 66, 121–134. Ekström, G., 1927. Klassifikation av Svenska Åkerjordar (Classification of Swedish Arable Soils) Ser C. No. 345 (Årsbok 20). Sveriges Geologiska Undersökning, Uppsala (in Swedish). (161 pp). Eriksson, J., Mattsson, L., Söderström, M., 2010. Tillståndet i svensk åkermark och gröda. Naturvårdsverket Rapport 6349 (in Swedish) 131 pp. Available online: https://www. naturvardsverket.se/Documents/publikationer/978-91-620-6349-8.pdf (verified: October 26, 2017). FAO, 2017. GSP Soil Data Policy. Available online: http://www.fao.org/3/a-bs975e.pdf (verified: October 26, 2017). Gee, G.W., Bauder, J.W., 1986. Particle-size analysis. In: Klute, A. (Ed.), Physical and Mineralogical Methods. Soil Science Society of America, Madison. Gustafsson, K., 1999. Models for precision application of lime. In: Stafford, J. (Ed.), Precision Agriculture ‘99 Odense. Denmark, SCI.
We describe production of the digital soil map of Sweden (DSMS), a map of topsoil texture and SOM content in arable land created through the use of data mining techniques whereby a set of detailed covariates were combined with a large number of reference samples. We show that transparency in the validation procedure is crucial if users are to understand the accuracy of a map product. Access to a large validation dataset consisting of soil analyses at farm level throughout southern Sweden made it possible to demonstrate the performance of DSMS at multiple scales. However, although excellent validation statistics were obtained at the overall scale, these figures reveal very little on the performance on farms or within fields. The clay content map produced was more accurate than the sand content map. The SOM content could not be calibrated with acceptable accuracy for precision agriculture from the available predictor data, but DSMS SOM has been used to produce 1 km × 1 km grid data for the FAO GSOC map. Through the validation procedure, we showed that DSMS is a versatile product that can be applied at different scales and for different objectives. For example, it can be used within farms as decision support for site-specific management (i.e. precision agriculture) and it can be used at municipality, county and national level as input to environmental impact scenarios and modelling of soil processes. Further developments 8
Geoderma xxx (xxxx) xxx–xxx
K. Piikki, M. Söderström
artiklar/ra1519.html (verified: October 26, 2017). Piikki, K., Söderström, M., Stenberg, B., 2013. Sensor data fusion for topsoil clay mapping. Geoderma 199, 106–116. R Core Team, 2017. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria Available online: https:// www.R-project.org (verified: October 26, 2017). Rawls, W.J., Pachepsky, Y.A., Ritchie, J.C., Sobecki, T.M., Bloodworth, H., 2003. Effect of organic carbon on soil water retention. Geoderma 116, 61–76. Söderström, M., Eriksson, J., 2013. Gamma-ray spectrometry and geological maps as tools for cadmium risk assessment in arable soils. Geoderma 192, 323–334. Söderström, M., Sohlenius, G., Rodhe, L., Piikki, K., 2016. Adaptation of regional digital soil mapping for precision agriculture. Precis. Agric. 17, 588–607. Söderström, M., Piikki, K., Cordingley, J., 2017. Improved usefulness of continental soil databases for agricultural management through local adaptation. S. Afr. J. Plant Soil 34, 35–45. Stoorvogel, J.J., Bakkenes, M., Temme, A.J., Batjes, N.H., Brink, B.J., 2017. S-world: a global soil map for environmental modelling. Land Degrad. Dev. 28, 22–33. Van der Klooster, E., Van Egmond, F.M., Sonneveld, M.P.W., 2011. Mapping soil clay contents in Dutch marine districts using gamma-ray spectrometry. Eur. J. Soil Sci. 62, 743–753. Van Veen, J.A., Kuikman, P.J., 1990. Soil structural aspects of decomposition of organic matter by micro-organisms. Biogeochemistry 11, 213–233. Weiss, A., 2001. Topographic Position and Landforms Analysis. Poster presentation, ESRI User Conference, San Diego, CA. Available at: http://www.jennessent.com/ downloads/tpi_documentation_online.pdf (Verified: October 10, 2017).
Hastie, T., Tibshirani, R., Friedman, J., 2009. The Elements of Statistical Learning: Data Mining Inference and Prediction, 2nd edition. Springer Series in Statistics, New York, pp. 383–411. Hengl, T., de Jesus, J.M., Heuvelink, G.B., Gonzalez, M.R., Kilibarda, M., Blagotić, A., et al., 2017. SoilGrids250m: global gridded soil information based on machine learning. PLoS One 12, e0169748. IAEA, 2003. Guidelines for Radioelement Mapping Using Gamma Ray Spectrometry Data. International Atomic Energy Agency, Vienna Austria IAEA-TECDOC-1363. Available online. http://www-pub.iaea.org/mtcd/publications/pdf/te_1363_web.pdf (verified: October 26, 2017). Janssen, P., Heuberger, P., 1995. Calibration of process-oriented models. Ecol. Model. 83, 55–66. Kätterer, T., Andrén, O., Jansson, P.-E., 2006. Pedotransfer functions for estimating plant available water and bulk density in Swedish agricultural soils. Acta Agr. Scand. B-S P 56, 263–276. Malone, B.P., McBratney, A.B., Minasny, B., 2013. Spatial scaling for digital soil mapping. Soil Sci. Soc. Am. J. 77, 890–902. McBratney, A.B., Santos, M.L.M., Minasny, B., 2003. On digital soil mapping. Geoderma 117, 3–52. Milborrow, M., 2017. Package ‘Earth’. Available online: https://cran.r-project.org/web/ packages/earth/earth.pdf (verified: October 26, 2017). Nash, J.E., Sutcliffe, J.V., 1970. River flow forecasting through conceptual models part I -a discussion of principles. J. Hydrol. 10, 282–290. Paulsson, R., Djodjic, D., Carlsson, Ross C., Hjerpe, K., 2015. Nationell jordartskartering Matjordens egenskaper i åkermarken. Swedish Board of Agriculture Report: RA15:19 (in Swedish). 44pp. Available online: http://webbutiken.jordbruksverket.se/sv/
9