This manuscript is a
non-peer reviewed
preprint that has been submitted for publi-
cation. Subsequent versions of this manuscript may have updated content. Feedback
and comments are welcomed, feel free to contact the corresponding author:
Alexandre Wadoux
Machine learning for digital soil mapping: applications,
challenges and suggested solutions
Alexandre M.J-C. Wadoux
a,
∗
, Budiman Minasny
a
, Alex B. McBratney
a
a
Sydney Institute of Agriculture & School of Life and Environmental Sciences, The University of
Sydney, Australia
Abstract
The uptake of machine learning (ML) algorithms in digital soil mapping (DSM)
is transforming the way soil scientists produce their maps. Machine learning is cur-
rently applied to mapping soil properties or classes much in the same way as other
unrelated fields of science. Mapping of soil, however, has unique aspects which require
adaptations of the ML algorithms. These features are for example, but not limited
to, the inclusion of pedological knowledge into the ML algorithm, the accounting
of spatial structure present in the soil data, or the desire to increase our scientific
understanding of the distribution and genesis of soil from a calibrated ML model.
Tackling these challenges is critical for machine learning to gain credibility and sci-
entific consistency in soil science. In this article, we review the current applications
of machine learning in digital soil mapping and suggest improvements. We found a
growing interest of the use of ML in DSM. Most studies focus on obtaining accurate
maps and disregard the characteristics of soil data, such as spatial autocorrelation.
Only a few studies account for existing soil knowledge or quantify the uncertainty
of the predicted maps. We then discuss the challenges related to the application of
ML for soil mapping and offer solutions from existing studies in the natural sciences.
The challenges are organized as follows: sampling, resampling, accounting for the
spatial information, multivariate mapping, uncertainty analysis, validation, integra-
tion of pedological knowledge and, interpretation of the models. We conclude that
for future developments, machine learning should incorporate three core elements:
plausibility
,
interpretability
, and
explainability
, which will trigger soil scientists to
move beyond model prediction and towards explanation of soil processes.
∗
Corresponding author: Sydney Institute of Agriculture & School of Life and Environmental
Sciences, The University of Sydney, New South Wales, Australia
Email address:
alexandre.wadoux@sydney.edu.au
(Alexandre M.J-C. Wadoux)
Preprint (Review article)
February 6, 2020
Keywords:
Soil science, Pedometrics, Data mining, Spatial data, Geostatistics,
Random forest
1. Introduction
In recent years, soil science has witnessed a considerable increase in digital soil
1
mapping activities. This is caused by the convergence of several timely factors which
2
are, among others, a huge demand for quantitative and spatial soil information, the
3
accumulation of databases of measured or inferred soil properties coupled with ex-
4
haustively known environmental variables and the development of numerical models
5
combined with computer resources to mine these stores of soil data. The digital soil
6
mapping (DSM) framework was formalized by the publication of
7
) which builds on Jenny’s
S
=
clorpt
model (
) of soil formation,
8
where
S
is the soil and the acronym
clorpt
stands for climate, organisms, relief, par-
9
ent material and time, respectively. In short,
clorpt
is a list of variables which, if
10
they are known without error, are likely to explain the soil variation over a region.
11
) supplemented Jenny’s formulation with
n
, which stands for
12
spatial position, and advocated the
scorpan
model for soil spatial variation. This
13
updated equation provides a spatial model to express quantitatively the relationship
14
between a soil property or class and environmental variables, for a given spatial lo-
15
cation.
16
17
Conventionally, spatial prediction of soil has been embedded in the geostatisti-
18
cal framework (
) in which a sample of a soil property is
19
modelled as a sum of a linear combination of environmental covariates and a spa-
20
tially autocorrelated (stochastic) residual, and prediction at unobserved locations is
21
made by kriging. Geostatistical models are often used in soil mapping because they
22
have several advantages (
). First, a statistically sound model is assumed
23
for spatial variation. This enables interpretation of the underlying physical processes
24
conveyed (inferred) by the model. Secondly, spatial autocorrelation is explicitly mod-
25
elled. This is relevant for environmental variables such as soil which vary from place
26
to place, but exhibit correlation between places. Thirdly, an explicit measure of the
27
uncertainty is associated with the prediction. In many circumstances such as in a
28
decision making process, the prediction is not the only interest and uncertainty maps
29
are required for the evaluation of the map quality or modelling risk.
30
31
Geostatistical mapping of soil has, conversely, several limitations which have only
32
partially been resolved in the current literature. To begin, the residuals are as-
33
2
sumed normally distributed, stationary (with constant mean and unit variance) and
34
isotropic. Next, modelling the non-linear relation between a soil property or class
35
and numerous cross-correlated covariates is not straightforward and introduces addi-
36
tional challenges (e.g. many parameters have to be estimated). Finally, geostatistical
37
models are computationally demanding if the sample size and/or the number of pre-
38
diction locations are large (
39
40
As an alternative, machine learning (ML) emerged in the 1990s as a tool for
41
spatial prediction and digital soil mapping (
). Machine learning
42
techniques refer to a large class of non-linear data-driven algorithms employed pri-
43
marily for data mining and pattern recognition purposes, and now frequently used for
44
regression and classification tasks in all fields of science. ML algorithms do not make
45
an assumption of the observations’ distribution, unlike geostatistical methods where
46
transformation of the original observations is often required to satisfy the assump-
47
tions. ML algorithms can also handle a large number of cross-correlated covariates
48
as predictor.
49
50
In parallel, there has been a tremendous increase in the production and availabil-
51
ity of regional and global soil databases. For example, the Soil and Terrain Digital
52
Database (SOTER,
)) made by FAO-UNESCO com-
53
piled quantitative information on soil and terrain for different parts of the world while
54
WoSIS is a harmonised database of more than 6 million geo-referenced soil records
55
). Additionally, numerous spatially exhaustive
scorpan
covariates
56
are available at global scale for climate (
57
). Further potential
58
covariates are provided by remote sensing such as by the MODIS (
59
satellite or Sentinel-2A hyperspectral sensor (
). Soil mappers
60
are now confronted with an increasing complexity in both soil data and covariates.
61
Conventional regression techniques seem, to some extent, outdated to accommodate
62
the increased complexity of soil datasets. This justifies the increasing use of machine
63
learning algorithms for digital soil mapping.
64
65
An essential distinction between conventional (statistical and geostatistical) mod-
66
els and ML algorithms applied in DSM is their purpose. Machine learning algorithms
67
mostly emphasize prediction accuracy whereas statistical models infer the process
68
which generated the data through a pre-defined model of spatial variation. In the
69
latter case, any interpretation is made in light of the model functions and the value of
70
the covariates or input data. In machine learning, a predictive model is constructed
71
3
to predict a set of input values to output values using an error-minimization proce-
72
dure. Since ML algorithms are not conditioned to follow any statistical assumptions,
73
they often appear more accurate than conventional models. The exact path between
74
input and output is ignored, and may not resemble an actual process described by
75
the existing knowledge. In soil science, the explosion of articles using ML algorithms
76
have made difficult to see the difference between model fitting and inference, and,
77
as a result between data science and soil science. Research seems to be driven by
78
the technique rather than by the hypothesis to be tested. This seems a poor bet
79
for the advancement of knowledge since “almost invariably the technician’s skill is a
80
solution looking for a problem” (
81
82
In DSM, the use of ML algorithms has led to an increasing number of publica-
83
tions where prediction (viz. mapping) of a soil property or class is the main interest.
84
Many “easy-to-follow” software implementations have supported this increase. Dig-
85
ital soil mapping, however, has unique characteristics which require adaptation of
86
the ML algorithms. These features are for example, but not limited to, the inclusion
87
of pedological knowledge in the ML algorithm, the accounting of spatial structure
88
present in the raw soil data, or the need to increase our scientific understanding of
89
the soil from a calibrated ML model.
90
91
This article aims to review the development of ML applied to digital soil mapping
92
by identifying key challenges and opportunities to solve them from the literature. In
93
this review, we define ML as the computer assisted practice of using data-driven (and
94
mostly non-linear) algorithms which resort to a large amount of calibration data to
95
learn a pattern and make a prediction. We start by reviewing and summarizing the
96
current use of machine learning in DSM. Based on this summary, we identify gaps in
97
the knowledge and define areas in which adapting ML algorithms would be beneficial
98
for their use in DSM. We propose solutions and a framework based on the literature
99
from different fields of natural science. Finally, we define three core elements that
100
should trigger soil scientists to move from model prediction to explanation of soil
101
processes.
102
2. A summary of applications
103
2.1. Extent, resolution, depths
104
Table
summarizes some recent case studies of digital soil maps that have been
105
produced using a ML algorithm. There is a large range of case studies, mapping soil
106
properties or classes from the plot (
<
1 km
2
) to the global (
>
10
7
km
2
) scale. Most
107
4
studies in our literature review predict at a local to regional scale. The mean extent
108
of the study area is 3,900 km
2
, but most (90%) studies consider a study area smaller
109
than 650,000 km
2
(equivalent to the size of metropolitan France). Few studies map
110
at plot or global scales. For example,
) make a quantitative map
111
over a 10 ha (0.1 km
2
) field in Denmark while
) produce quanti-
112
tative and categorical maps for the whole world.
113
114
We found a clear correlation between the spatial extent of the study area and the
115
grid spacing (i.e. the spacing between point predictions) at which the soil property or
116
class is mapped: the larger the study area, the coarser the resolution. The resolution
117
spans between 2 m
×
2 m (
) to 1 km
×
1 km for large, regional
118
or continental study areas (e.g.
). Most studies, however, map at a
119
standard spatial resolution of 30, 90 or 250 m.
120
121
While most of the studies (70%) predict a soil property or class for a single depth
122
(topsoil), a number of studies accounts for the soil variation at multiple depths.
123
) follow the GlobalSoilMap project specifications (
124
) to produce a quantitative three dimensional map of several soil
125
properties for six depths intervals, namely 0-0.05 m, 0.05-0.15 m, 0.15-0.30 m, 0.30-
126
0.60 m, 0.60-1.00 m and 1.00-2.00 m. Similar depth intervals are used in
127
) and
) for soil organic prediction in France or Denmark,
128
respectively. Several other studies (e.g.
) use
129
standard depth intervals for prediction, based on national mapping requirements or
130
suitable for their specific case study.
131
2.2. Sampling design, sample size and density
132
The sampling design is the spatial location of the sampling units used to cal-
133
ibrate or validate the ML algorithm.
Most studies do not specify the sampling
134
design used to generate the observations. It is speculated that the sample originates
135
from multiple sources, e.g. legacy data, expert-based designs, and combination of
136
several surveys, each of which had a different sampling design. When specified, non-
137
probability sampling such as grid-based sampling designs are by far the most used
138
(e.g. by
Pahlavan-Rad & Akbarimoghaddam
139
). Another non-probability sampling design is conditioned Latin Hyper-
140
cube (cLHS), used to collect a sample in
141
Probability sampling is used in about one fourth of the studies. For example, simple
142
random sampling is used in
), while a sample is collected based
143
on stratified random sampling in
) using land use and topog-
144
5
raphy as stratifying variables.
145
146
In our literature review, we found that the sample size varies considerably be-
147
tween studies. While the average sample is composed of 1,000 units, about one third
148
of the studies use a sample with less than 150 units, mostly for local or small-scale
149
regional areas. For example,
) use a sample of size 47 for mapping
150
soil water retention in a 93 km
2
area while
) observed 33 soil
151
profiles to calibrate a ML algorithm and to predict soil taxa over a 11,600 km
2
area.
152
As expected, global studies have very large sample sizes.
) and
153
) use a sample composed of more than 150,000 units to make
154
soil property or class maps of the whole world, or of the United States, respectively.
155
156
When the sample size is associated to the extent of the study area, our review
157
shows that large-scale studies have a very coarse sampling density. While the average
158
sampling density in our literature is 0.24 units/km
2
, studies by
159
and
) have both a sampling density smaller than 3 units/10,000 km
2
160
for mapping soil properties in the rangelands of eastern Australia or in the Canadian
161
boreal forests. Small-scale studies have, conversely, high sampling density. All studies
162
with area size less than 50 km
2
have a sampling density larger than 7 units/km
2
.
163
2.3. What is mapped?
164
2.3.1. Quantitative variables
165
ML algorithms have been successfully applied for quantitative mapping of vari-
166
ous soil properties such as soil organic carbon concentration (
167
168
169
), to map soil texture (viz. clay, silt
170
and sand content) (
171
), or cation exchange
172
capacity (
173
174
ML algorithms have also been applied to make maps of soil nutrients such as
175
nitrogen (
176
), potassium, calcium or
177
magnesium (
178
179
A number of studies have also predicted soil attributes and conditions with ma-
180
chine learning such as bulk density (
) or soil pollutants
181
6
) map soil background concentrations of arsenic
182
in the Jiangxi Province in China.
) map soil salin-
183
ity in Iran.
) map soil invertebrate using environmental covariates
184
in a deciduous forest ecosystem in northern Iran while
) map
185
carbon storage and available water capacity in an area in eastern Australia.
186
2.3.2. Categorical variables
187
Compared with continuous soil property mapping, fewer studies apply ML to
188
categorical variables. Digital mapping of soil classes using machine learning started
189
in the 90s. Probably the first of its kind,
) predict soil
190
classes in a regional area while
) predict soil drainage classes
191
using remote sensing and elevation covariates.
) map soil units
192
in a 600 km
2
area of Western Germany. These studies have recently been completed
193
by a number of publications comparing the maps predicted by a ML model to con-
194
ventional soil maps (e.g.
195
) employ machine learning
196
to classify soil taxonomic units.
) map salt-affected
197
areas in irrigation schemes in South Africa. Table
provides an additional summary
198
of case studies.
199
200
A special case of categorical mapping occurs when the map of soil class already
201
exists but needs to be disaggregated.
) and
202
use a decision tree to disaggregate an existing map and obtain a realization of the
203
disaggregated soil class distribution. With multiple realizations, the most probable
204
soil class is obtained for a given location. This is further investigated by
205
) to disaggregate a reconnaissance soil map using a binary decision tree.
206
A similar approach with decision tree is used in
) to downscale
207
soil types within existing map unit boundaries. More recently,
208
use ML to model and disaggregate soil classes and report the probability associated
209
to each soil class at a given location in the area of interest. A growing number of
210
publication exploits the DSMART approach proposed by
) (e.g.
211
212
2.4. Covariates
213
Environmental covariates are used as predictors in ML algorithms. They are
214
supposed to explain part of the physical and chemical process governing soil spa-
215
tial variation. Most studies use about 20 covariates. Only a few use less than five
216
(e.g.
) while other use more than 100 (e.g
217
). Since the covariates represent soil forming fac-
218
tors, numerous studies (e.g.
219
) logically select the covariates to represent the
220
key factors of the
scorpan
model of soil spatial variation. The most common ones
221
are existing soil property or class maps, (long-term) average annual precipitation
222
and temperature, remote sensing images (e.g. SPOT satellite images or vegetation
223
indices derived from satellite images), elevation, terrain attributes (e.g. slope, local
224
curvature, topographic wetness index) and existing geological maps.
225
226
Covariates representing
scorpan
factor of soil variation might not be available or
227
easily obtainable in all case studies. In some cases, covariates are chosen based on
228
expert knowledge. A number of studies therefore calibrate machine learning algo-
229
rithms using sets of climatic variables, remote sensing images or terrain attributes
230
only, or a combination of them. For example,
) use a set of
231
eight climatic and eight terrain attribute variables to map C, N and soil texture in a
232
large area in Canada.
) use six terrain attributes as predictors.
233
There are chosen from a large set of environmental covariates using knowledge on
234
the expected relationship between the covariate and the soil property to be mapped.
235
We note that a few studies (e.g.
) consider
236
that if a sufficiently large (
>
100) number of covariates is used for calibration, the
237
machine learning algorithm learns a representation of the spatial pattern and pre-
238
dicts a realistic spatial pattern. This large amount of covariates relies mostly on
239
remote sensing images, e.g. MODIS land products (long-term averages, several near-
240
or mid-infrared bands) or Landsat products (near-, short-wave near-infrared, or
γ
241
radiometric bands, bare ground images).
242
243
A few studies account for the multi-scale variation of the environmental covari-
244
ates. In other words, terrain derivatives may well be aggregated to account for
245
physical processes in soil that are not visible are finer scale. Examples of studies us-
246
ing multi-scale covariates for mapping with machine learning algorithms are
247
) or more recently
248
), for example, use a total of 412 covariates, several of which are derived
249
from the aggregation of terrain attributes from a fine (i.e. a grid cell size of 2 m
×
250
2 m) elevation map.
251
252
A growing number of studies have advocated the use of spatial surrogate covari-
253
ates as an indicator of spatial position in the
scorpan
model of soil variation The
254
most common surrogate is the use of geographical coordinates (easting and northing)
255
8
as covariates in the model. Maps of distances from observation locations, or group
256
of locations, have been used by
). They are categorized into Eu-
257
clidean, downslopes or “resistance” distances. More recently,
258
use Euclidean distance fields, which are maps of distance from reference locations in
259
the study area such as a corner or a center.
260
2.5. Covariate selection
261
Covariate (
aka
feature) selection aims at reducing the number of covariates used
262
to calibrate the machine learning models. While most ML models are robust to mul-
263
ticolinearity between covariates, there are several reasons for selecting a subset of
264
covariates to calibrate the model. Some of them are: (i) to calibrate the ML model
265
faster, (ii) to reduce complexity, (ii) to increase the prediction accuracy or (iv) to
266
prevent over-fitting of the ML model, i.e. to prevent poor prediction accuracy on
267
unseen data. In our literature review, about one third of the studies apply covariate
268
selection. Two main categories of covariate selection techniques are found. The first
269
applies the covariate selection as a pre-processing step, i.e. before calibrating the ML
270
model. This is the case in
271
) select the covariates to be used in calibration
272
by computing the Pearson’s
r
correlation coefficient between the covariates, and by
273
discarding the ones that were highly correlated, while
) select the
274
covariates based on the Pearson
r
correlation coefficient between the soil property
275
values and the covariates, and select a subset of covariates which are strongly corre-
276
lated with the property. The second type of covariate selection are called “wrapper”
277
methods and rely on the inference made by a calibrated ML model to determine
278
whether covariates are important. By re-calibrating a ML model several times, each
279
time removing the least important covariate, one may expect to reduce considerably
280
the overall number of covariates with little or no decrease in model prediction accu-
281
racy. Examples on the use of “wrapper” methods are found in
282
) or
283
). The most used of “wrapper” methods is an optimization algorithm
284
called recursive feature elimination.
285
286
2.6. Machine learning models
287
A large number of ML algorithms and their variants have been used in the DSM
288
literature.
For quantitative mapping, tree-based algorithms are the most popu-
289
lar ones, the simplest version of which is the regression tree, used for example by
290
). Regression tree is known to be sensitive to the
291
9
calibration sample. To solve this problem, the bagging (bootstrap and aggregating)
292
procedure (
) has been introduced in random forest (RF). Our litera-
293
ture review shows that RF is currently the most popular ML algorithm for regression
294
purposes. Example of case studies using RF for mapping are
295
296
). More recently,
) introduced a variant of
297
random forest, called quantile regression forest, as a method to map the uncertainty
298
associated with the prediction of the soil property. Another tree-based method is
299
cubist, employed in about 10% of the reviewed literature (e.g. by
300
). A few studies (less than five) use
301
boosted regression tree (
). In addition, a num-
302
ber of studies use neural networks (
303
), such as artificial neural networks (
304
). A relatively small number of studies use alternative algorithms such as sup-
305
port vector machines (
k
-nearest neighbours (
306
) or generalized boosted regression (
307
308
For classification purposes, tree-based algorithms are also the most popular ones.
309
About 80% of the case studies used at least one tree-based algorithm such as regres-
310
sion tree (e.g.
), random forest
311
(e.g.
) or boosted regression tree (e.g.
). Al-
312
ternatively, gradient boosting is used by
k
-nearest neighbors by
313
) and compared to support vector machines. The
314
latter algorithm is also used in
). Neural networks
315
is also popular and used in
316
317
Recent studies have proposed to use model ensemble techniques to improve the
318
predicted map of several individual models in terms of accuracy.
319
) combined seven ML model predictions for soil class mapping in a case
320
study in Iran while
) implemented a weighted ensemble learning
321
model to map soil organic carbon in consideration of pedoclimatic zones in China.
322
Ensembles are also considered in
) for global soil mapping.
323
2.7. Parameter tunning
324
The performance of a machine learning model is impacted by the values of its
325
model parameters. While most ML would perform well on default tuning parameter
326
values, almost half of the studies perform a search to find optimal values.
327
) manually decide the artificial neural network neurons number for
328
10
each layer of the network. This manual search is automated by a so-called grid-
329
search process. This is by far the most used technique for parameter tuning. In a
330
grid-search process, a number of parameter values are evaluated based on the model
331
prediction error. The process is computationally intensive (the ML model must be
332
calibrated for each parameter set proposal). Examples of studies using a grid-search
333
to find ML parameter values are
334
Pahlavan-Rad & Akbarimoghaddam
335
). An alternative to the grid search is to apply
336
an optimization algorithm, such as the particle swarm method, to find optimal pa-
337
rameter values. For example,
) compare two genetic algorithms and
338
a grid search process to find the ML parameters. Recently,
339
use Bayesian optimization to optimize the number of layers, the neuron number,
340
the learning rate and the batch size of an artificial neural network for mapping soil
341
organic carbon.
342
2.8. Validation and uncertainty quantification
343
In our literature review, all studies compute at least one validation statistic to
344
assess the quality of the prediction. A list of validation statistics is provided in
345
Table
. About 30% of the studies obtain the validation statistics through cross-
346
validation, while 30% through data-splitting. The remaining studies either repeat
347
data-splitting several times, validate through visual examination or use a grid-based
348
sampling design. Only two studies collect an additional probability sample for vali-
349
dation (
350
351
In addition to the validation statistics, about 30% of the studies quantify the
352
uncertainty associated with the prediction. These studies report confidence inter-
353
val, obtained by bootstrapping the original set of observations (e.g.
354
). A few studies use the kriging
355
variance computed on the residuals of a trend obtained by predicting with a ML al-
356
gorithm (e.g.
), or a combination of bootstrap and kriging variance
357
(e.g.
). In three studies, prediction intervals are obtained
358
through the quantile regression forest.
) obtain the prediction inter-
359
vals following a two-step procedure called mean plus variance estimate for mapping
360
several soil properties using an artificial neural network.
361
11
Table 1: Non-exhaustive list with summary of case studies in which machine learning algorithms are used for digital soil mapping.
Spatial
extent
Sample size
Sampling design
Number of
covariates
Machine learning
model
Covariate
selection
Parameter
tuning
Validation statis-
tics
Uncertainty
quantification
Reference
Quantitative maps
Plot
285
grid-based
19
cubist, RF
no
no
R
2
, RMSE
no
Local
47
stratified random
41
RF
yes
yes
RMSE, IQR
yes
Local
70
cLHS
19
cubist
no
no
MAE, RMSE, R
2
,
CCC
yes
Local
75
grid-based
9
ANN
no
no
R
2
, MSE
no
Local
98
varied sources
173
RF
yes
no
RMSE, R
2
no
Local
116
simple random
20
RF
no
no
R
2
, RMSE, CCC
no
Local
117
not specified
13
GBM
yes
yes
R
2
, RMSE, MAE
yes
Local
117
not specified
412
cubist
yes
no
ME, MAE, R
2
,
R
2
adj
no
Local
120
stratified random
not speci-
fied
RF
no
no
ME, RMSE, R
2
,
MSE
no
Local
120
stratified random
22
ANN, BRT
yes
yes
R
2
, RMSE, ME
no
Local
137
systematic random
20
ANN, GEP
yes
yes
RMSE, R
2
, MBE
no
Local
138
not specified
15
RF
yes
no
RMSE, R
2
, CCC
no
Local
150
grid-based
not speci-
fied
ANN
no
yes
correlation coeffi-
cient, R
2
, RMSE,
Willmott’s in-
dex of agreement,
RPIQ
no
Local
151
not specified
not speci-
fied
no
no
R
2
, NRMSD
no
Local
153
grid-based
26
RF
yes
no
RMSE, R
2
no
Local
159/34
not specified
37
RF, cubist, QRF, NN,
avNNet, ctree, evtree,
GBM,
k
-NN, RT,
SVM
yes
no
R
2
, RMSE, MAE,
MARE
yes
Local
165
stratified random
18
RF
no
yes
MSE, NMSE
no
Local
173 profiles
cLHS
19
Rf
no
no
ME, RMSE, R
2
no
Local
188 profiles
cLHS
16
ANN, SVR,
k
-NN,
RF, RT
no
yes
RMSE, CCC
no
Local
234
not specified
410
cubist
yes
no
MAE, R
2
yes
1
Plot: 0-1 km
2
; Local:
>
1 km
2
–10
4
km
2
; Regional:
>
10
4
km
2
–10
7
km
2
; Global:
>
10
7
km
2
.
2
RF: random forest; ANN: artificial neural networks; CNN: convolutional neural networks; GBM: gradient boosting machine; BRT: boosted
regression tree; GEP: gene expression programming; QRF: quantile regression forest; avNNet: neural networks using model averaging; ctree: conditional
inference trees; evtree: evolutionary algorithm for classification and regression tree; NN: neural networks; GBM: generalized boosted regression;
k
-NN:
k
-nearest neighbors; RT: regression tree; SVM: support vector machine; MARS: multivariate adaptive regression splines; SGB: stochastic gradient
boosting; CART: classification and regression tree; NSC: nearest shrunken centroids; CT: classification tree; BCT: bagged classification tree; DT:
decision tree; LMT: logistic model tree; EGB: extreme gradient boosting.
3
R
2
: coefficient of determination; R
2
adj
: adjusted coefficient of determination; RMSE: root mean square error; IQR: interquartile range; MAE:
mean absolute error; CCC: Lin’s concordance correlation coefficient; MSE: mean square error; ME: mean error; MBE: mean bias error; RPIQ: ratio of
performance to interquartile distance; NRMSD; normalized root mean squared deviation; MARE: median absolute relative error; NMSE: normalized
mean square error; sMAPE: symmetric mean absolute percentage error; SS: skill score; RMSD: minimum root mean square deviation; RPD: residual
prediction deviation; SDE: standard deviation of the error; EC: overall ratio; OA: overall accuracy; PA: producer accuracy; UA; user accuracy; AUROC:
area under receiver operating characteristic curve; AUC: area under the curve.
12
Local
330 profiles
not specified
12
BRT, ANN, least-
square SVM
no
yes
R
2
, R
2
adj
,RMSE,
relative RMSE
no
Local
330
simple random
10
RF, GBM
no
yes
ME, MAE, RMSE,
R
2
Local
334
cLHS
16
cubist, RF, RT
yes
no
R
2
, RMSE
no
Local
342/321
-
14
MARS, SVR, RF,
Cubist, NN
-
yes
R
2
no
Local
399
not specified
12
RF
no
no
R
2
, RMSE
no
Local
440
varied sources
19
RF, SVM, ANN
no
yes
RMSE, ME
no
Local
460
grid-based
21
RF
no
yes
ME, MAE, RMSE
no
Pahlavan-Rad & Akbarimoghad-
dam
Local
568
simple random
26
QRF
no
no
R
2
, RMSE, range-
normalized RMSE,
Moran’s I
yes
Local
1104
expert
29
RF, SVM, SGB
no
yes
RMSE, sMAPE
no
Local
≤
1052/2050/
2379
varied sources
300-500
BRT, RF
yes
yes
bias, RMSE, SS,
R
2
no
Local
2388
varied sources
3
CNN, RF
no
yes
ME, RMSE, R
2
,
CCC
no
Regional
not specified
not specified
20
cubist
no
no
R
2
, RMSE, bias,
CCC
yes
Regional
125 profiles
purposive
12
BRT, RF
no
no
MAE, RMSE,R
2
,
CCC
no
Regional
244
grid-based
4
ANN
no
yes
ME, MAE, RMSE,
CCC
no
Regional
339/961
varied sources
40
QRF
no
no
R
2
, RMSE
yes
Regional
485 profiles
not specified
5
CNN
no
yes
R
2
, RMSE
yes
Regional
500
not specified
12
RF, BRT
yes
no
R
2
, RMSE
no
Regional
528
subset from a
systematic grid
18
k
-NN
yes
no
RMSE, R
2
, Bias,
coefficient of vari-
ance
no
Regional
705
simple random
16
RF, BRT, SVM
yes
yes
R
2
, MAE, RMSE,
CCC
yes
Regional
978 profiles
not specified
24
RF
no
no
R
2
, ME, RMSE,
CCC
no
Regional
1,014
stratified random
327
CART, BRT, BRT,
RF, SVM
yes
no
R
2
, RMSD, RPD,
RPIQ
no
Regional
1,134
not specified
81
NN
no
no
R
2
, ME, MAE,
RMSE
no
Regional
1,300 profiles
not specified
6
RF
no
no
CCC, RMSE
yes
Regional
1,626
not specified
40
SVM
no
yes
R
2
, MSE
no
Regional
2,024 profiles
legacy data
16
QRF
no
no
ME, RMSE, R
2
,
accuracy plot
yes
Regional
2,024 profiles
legacy data
16
no
yes
MSE, R
2
no
Regional
2,943
two-stage system-
atic
37
CNN, RF
no
yes
ME, RMSE, R
2
,
CCC
yes
Regional
4,859
not specified
26
QRF
no
no
ME, RMSE, accu-
racy plot
yes
Regional
4,859
not specified
32
QRF
no
no
ME, RMSE, accu-
racy plot
yes
Regional
5,386
varied sources
6
cubist, SVM
no
no
R
2
, MSE, CCC
Regional
13,000
not specified
18
RF
no
no
R
2
yes
Regional
19,790
two-stage system-
atic
197
RF
no
no
ME
no
Regional
37,693
legacy soil data
74
RF, Cubist, SVM
yes
yes
R
2
, RMSE, MAE
yes
Regional
- Global
2,268-27,262
varied sources
34
cubist
no
yes
CCC, RMSE,
SDE, ME
yes
13
Regional
- Global
366,034
varied sources
>
200
RF, GBM
no
yes
R
2
,ME, RMSE,
MAE
yes
Global
11,268
legacy soil data
118
SVM, kernel weighted
NN, RF
yes
no
EC, RMSE, R
2
yes
Global
150,000
legacy soil data
>
200
RF, GBM
no
yes
R
2
no
Categorical maps
Local
-
not specified
125
ANN
no
no
Accuracy, recall,
precision
no
Local
33 profiles
not specified
16
RF, J48
no
no
not specified
no
Local
103/297/ 57
cLHS
130
k
-NN, NSC, CT,
BCT, RF, linear SVM,
radial-basis SVM, NN,
ANN
yes
yes
Kappa analysis,
Brier scores, vi-
sual inspection,
confusion index
no
Local
125 profiles
cLHS
17
RF
no
no
map purity, Co-
hen’s kappa, Shan-
non entropy index,
relative purity,
relative diversity
no
Local
151
not specified
not speci-
fied
SVM
no
no
NRMSD, mi-
cro averaged F1
measure, kappa
statistics
no
Local
175, 63 profiles
varied sources
27
k
-NN, SVm, DT, RF
no
no
OA, PA, UA,
kappa coefficient,
AUROC
no
Local
452 profiles
regular grid
6
DT, RF
yes
no
OA, UA, PA,
Kappa coefficient
of agreement
no
Local
917
grid-based
33
RF
yes
no
Kappa index
no
Local
3,121
by-polygon,
equal-class, area-
weighted, and
area-weighted
with random over
sampling
20
CART, CART with
bagging, RF,
k
-NN,
NSC, ANN, LMT,
SVM
no
yes
overall agreement,
quantity disagree-
ment, allocation
disagreement, total
disagreement
no
Regional
89,323
random sampling
26
k
-NN, RF
yes
no
recall, accuracy
no
Regional
366,034
varied sources
>
200
RF, GBM
no
yes
OA, regional
dataset
yes
Regional
7,664 profiles
varied sources
110
DT, RF, EGB, SVM,
k
-NN
yes
no
OA, precision,
recall, F-score,
K-index
no
Regional
9,924
not specified
23
RF
yes
no
error matrix
no
Global
150,000
legacy data
>
200
RF, GBM
no
yes
map purity,
weighted kappa
metrics, AUC,
True positive rate,
scaled Shannon’s
entropy index
no
14
3. Challenges and opportunities
362
Based on the review, here we identify some knowledge gaps and challenges in
363
the current use of ML algorithms for DSM. We will outline some opportunities for
364
research.
365
3.1. Sampling
366
Despite abundant evidence that the sampling design and sample size play a key
367
role in the resulting map accuracy (
), sampling designs suit-
368
able for mapping with machine learning are yet to be uncovered. The impact of the
369
sample size for mapping with ML is discussed in
) where
370
the efficiency of several ML algorithms are compared for the spatial prediction of
371
soil carbon. The study shows that having a sufficiently large sample size is more
372
important than choosing a sophisticated ML algorithm, and that when the sample
373
size is small, it it best to use simple models. About sampling designs,
374
speculates that machine learning algorithms would benefit from a spread of the sam-
375
pling units in the feature (covariate) space, and suggests the use of feature space
376
coverage sampling (FSCS) using
k
-means clustering or conditioned Latin Hypercube
377
sampling (cLHS). Both sampling designs aim at covering the space spanned by the
378
covariates, but in different ways. Experimental results are provided by
379
) in a study comparing five sampling designs (viz. simple random sampling,
380
cLHS, spatial coverage sampling (SCS), FSCS and a design optimized in terms of
381
mean square error) for soil property mapping with random forest. The results show
382
large differences in mapping accuracy between the designs, and that a FSCS de-
383
sign optimized in the most important covariates of the random forest model had
384
the closest match to an optimized design. By performing further diagnostics, the
385
study concludes that RF does not benefit from a uniform spread of the units in the
386
geographic/feature space, nor from reproducing the marginal distribution of the co-
387
variates (as it is done in cLHS). These results apply for RF but there is a need to
388
further investigate sampling designs for other machine learning algorithms. While
389
most studies in our literature review (Table
) use a grid-based sampling or cLHS,
390
there is now evidence that most conventional sampling designs (e.g. spatial coverage
391
sampling) are not effective for the purpose of mapping with machine learning.
392
393
To discover what makes a good design for mapping with machine learning, one
394
should ideally derive optimal designs. More importantly, one should investigate the
395
characteristics of these designs, so that future research can generate simple designs
396
that resemble the optimal ones (
). It is likely that optimal designs
397
differ between machine learning algorithms. We speculate that a somewhat uniform
398
15
spread in the feature (i.e. covariate) space remains important for all ML algorithms
399
since they all link the covariates and the sample values in a non-linear way, but that
400
additional considerations might outweigh or overtake this uniform spread. An ex-
401
ample of optimal design is given by the studies of
402
and
) where the sampling configurations are optimized with active
403
learning for mapping with support vector machines. In the first study, the selected
404
sampling units are the most beneficial for the algorithm, avoiding mis-classification
405
between temperature below or above 20Cs (categorical mapping) by becoming sup-
406
port vectors. In
), a similar methodology is adopted and tested in
407
three case studies to subsample an existing sample for quantitative mapping, to add
408
optimally new sampling units in a continuous map or to define suitable areas for
409
sampling. In all case studies, the authors obtained a design optimal for the purpose
410
of mapping with support vectors machine. They conclude that while a sampling
411
design can be representative of the geographical space, the latter can be judged un-
412
representative if other dimensions are considered. These results encourage the use
413
of new methods for sampling design optimization such as active learning. Active
414
learning is a model-based sequential re-design algorithm. In active learning, the
415
objective function (e.g. the spatially averaged prediction uncertainty) is explicitly
416
quantified and used to define the additional sampling units that are the most ben-
417
eficial for the model (e.g. the boundary between two classes). In this sense, active
418
learning is similar to optimization with spatial simulated annealing routinely used
419
in geostatistical sampling design optimization. Besides the optimization algorithms,
420
a set of objective functions needs to be tested.
) defined an objective
421
function that searches for the optimal units in the space spanned by the predictors
422
(i.e. covariates) for prediction using a neural network algorithm. Taking the latter
423
considerations and testing active learning for sampling design optimization would
424
certainly make a valuable contribution to digital soil mapping research.
425
3.2. Resampling
426
Regional and global scale studies almost invariably use legacy soil data (
427
). Legacy soil samples provide valuable information on soil classes and
428
properties but are often highly clustered in areas of specific interest. In modelling
429
with machine learning, it is assumed that the sample is composed of independent
430
and identically distributed sampling units whereas soil observations within an area
431
typically exhibit spatial autocorrelation (i.e close observations are more similar than
432
remote ones) This has important implications in terms of sampling, resampling of
433
the observations and validation of the models. A ML algorithm calibrated with a
434
spatially clustered sample may lead to biased predictions over the area because of
435
16
the over-representation in the calibration process of regions of high sampling density.
436
Despite being critical, this has yet been disregarded in DSM studies. In geostatis-
437
tics, spatial declustering has been applied to reduce the effect of clustered data in the
438
calculation of experimental variogram (
). One form, called cell
439
declustering involves overlaying a grid over the area and assigning a weight to the
440
sampling units based on the inverse of the number of units in the cell. In ecology, a
441
first attempt was made by
) and later
) to decluster the
442
sampling units used in the calibration of a CART model. In
), weights
443
are given to the sampling units, where the weights are obtained from a kriging of
444
the spatial mean.
) elaborate a more complex procedure in which
445
all quantities involved in the CART algorithm (e.g. the proportion of leaves) have
446
a spatial estimate. This has been further considered by
) for
447
both categorical and quantitative mapping of ecological variables.
448
applied polygonal declustering technique to spatially clustered samples by assigning
449
weights on the units based on Voronoi’s area proportion.
450
451
We point out that the clustering may also occur in the feature (i.e. covariates)
452
space and speculate that this may also affect the prediction if most units are clustered
453
at some specific areas of the feature space. For example, a model trained to predict
454
organic carbon in a mountainous area will exhibit biased prediction if most sampling
455
units originate from valley, and that elevation is used as predictors. Similar to
456
), weights can be assigned to the units so down-weight the importance of
457
over-sampled areas in the feature space. An example method is provided by
458
). The authors assume that a good sample have a uniform spread in the
459
feature space and thus covers all strata of a hypercube based on the covariates. A
460
weight is assigned to each unit in the sample based on the density of the units in each
461
stratum. The larger the density within the stratum, the smaller the weight assigned
462
to a single unit.
463
464
The nature of the legacy soil data in categorical mapping also poses additional
465
challenges. ML algorithms for categorical mapping rely on balanced sets of units.
466
In other words, all classes shall comprise a comparable number of sampling units.
467
Legacy soil samples are considered imbalanced in that all classes are not represented
468
equally. Most ML algorithms are calibrated by maximizing the average (classifica-
469
tion) accuracy on an independent validation sample. This often results in very low
470
predictive accuracy for under-sampled classes, and models biased toward the over-
471
sampled classes (
). In the ML literature, several approaches have
472
been developed to handle class imbalanced samples. At the higher level, one may
473
17
distinguish between cost function and resampling based approaches. In the first ap-
474
proach, the model is penalized for miss-classification to under-represented classes.
475
This stems from the calibration of ML algorithms, which minimize a loss function
476
to find optimal parameter values (e.g. in neural networks). In the second approach,
477
resampling of the sample is performed by either adding units in the under-sampled
478
class, removing units from the over-sampled class, or a mix of the two. The second
479
approach has been recently been applied in soil mapping studies, in particular by
480
) and
481
tested eight resampling approaches and their effect on the prediction accuracy of five
482
ML algorithms in two large-scale case studies. However, to date resampling tech-
483
niques are applied the same way as in other disciplines while soil data often presents
484
spatial autocorrelation which may impact the resampling strategies. This has not yet
485
been investigated in the literature. The integration of resampling strategies within
486
a general framework for mapping with ML is provided in Fig.
487
3.3. Accounting for spatial information
488
Machine learning algorithms do not account for spatial autocorrelation contained
489
in the raw soil data, unless explicitly specified.
) have tested ran-
490
dom forest for different scenarios of spatial autocorrelation in the observations and
491
confirmed that the presence of spatial autocorrelation leads to high variance of the
492
residuals. ML algorithms accounting for autocorrelated observations have recently
493
been formulated, such as geographical random forest (
), or
494
spatial ensemble techniques (
). The two methods boil down to
495
geographically weighted regression by fitting spatially local sub-models using only
496
neighbouring observations.
) decomposed the area into geographic
497
disjoint sub-areas, and fitted a local model in each sub-area.
498
fitted a sub-model to each observation using random forest, accounting for both non-
499
stationarity and spatial autocorrelation.
500
501
Applying a non-spatial model for digital soil mapping is not a problem in itself.
502
This is corroborated by the definition of DSM given in
503
), which gives provision for mapping using “non-spatial soil inference systems”.
504
In theory, if one includes all relevant environmental variables to model the soil prop-
505
erty or class, there should be no spatial autocorrelation in the residuals of the fitted
506
models. If this happens, some important predictors are likely to be missing. More
507
importantly, this also means that predictions made by the ML algorithm might be
508
biased or the model underfitted because this is a violation of the assumption of inde-
509
pendence between data points that is implicitly assumed.
510
18
recommend mapping the spatial distribution of the residual autocorrelation to facili-
511
tate the identification of a missing spatial process. In some cases, a map of residuals
512
exhibits a clear pattern (e.g. increasing residuals with distance from the river) and
513
might help to generate a new hypothesis or to refine the existing model (see Fig.
514
515
Despite the availability of datasets and care made during modelling, residual au-
516
tocorrelation is still likely to occur. Several authors have advocated the use of spatial
517
surrogate covariates as an indicator of spatial position in the
scorpan
model of soil
518
variation or to account for spatial autocorrelation contained in the data. The most
519
common surrogate is the use of geographical coordinates (easting and northing) as
520
covariate in the model. This has led to maps with visible artefacts, in particular
521
when used in combination with tree-based algorithms. Alternatively, maps of dis-
522
tances from observation locations, or a group of locations, have been proposed by
523
). They are categorized into Euclidean, downslopes or “resistance”
524
distances. Maps of distance to observation locations generally have no direct mean-
525
ing in terms of soil process over an area (e.g. distance from the river).
526
) propose to use Euclidean distance fields, which are maps of distance from
527
reference locations in the study area such as the corner or the centre. The studies
528
using distance maps as covariates have shown for several case studies an important
529
reduction of the residual autocorrelation, when compared to a model without dis-
530
tance maps in the set of covariates.
531
532
In the context of digital soil mapping, we infer that the current use of distance
533
maps is not satisfactory for several reasons. Including pseudo-covariates with the set
534
of pedologically relevant covariates can be harmful because it precludes analysis of the
535
residuals and the generation of new hypotheses from these residuals (
536
It also hampers the interpretation of the most important predictors (
537
), which is key in several studies on soil mapping. Finally, pseudo-covariates
538
of distance may well integrate over several pedologically relevant covariates, making
539
them better predictors or masking the effect of pedologically relevant covariates. In
540
spatial ecology, alternatives to distance maps are found in the use of spatial eigenvec-
541
tor maps, spatial filters or trend-surface regression computed on, or optimized for,
542
the residuals of a model calibrated using ecologically relevant covariates (
543
). The process is generally in three steps (Fig.
). In the first step, the variable
544
of interest is fitted using ecologically relevant covariates, and the (autocorrelated)
545
residuals are mapped to investigate whether there is an obvious missing spatial pro-
546
cess in the model. In the second step, spatial surrogate covariates are computed on,
547
or optimized for, the residuals. Finally, a model is calibrated using the covariates
548
19
from steps 1 and 2.
549
550
Step 1 [a + b]
Step 2 [b + c]
Step 3 [a + b + c]
y
= f
(X)
+
ε
y
= f(
W
) +
ε
y
= f(
XW
) +
ε
soil property
of interest
regression
f
on
pedologically relevant
covariates
X
residuals,
possibly spatially
autocorrelated
soil property
of interest
regression
f
on spatial
covariates
W
residuals,
possibly spatially
autocorrelated
soil property
of interest
uncorrelated
residuals
regression
f
on
pedologically relevant
covariates
X
and spatial
covariates
W
variation of
y =
[a]
[c]
[b]
[d]
Variation explained by
X
Variation explained by
W
Unexplained
variation
Figure 1: The three steps of variation partitioning between environmental
X
and spatial covariates
W
. The variation of
y
is partitioned into four fractions (
) which are [a] the
variation due to the environmental covariates, [b] the variation due to the spatial component of the
environmental variables, [c] the spatial component and [d] the unexplained residual variation. Each
component is estimated using the amount of variance explained. All [a + b + c + d] sum to 1.
The main advantage is to enable subsequent interpretation of the role of environ-
551
mental covariates, spatial covariates (most often in the form of Moran’s eigenvector
552
maps) and unexplained (uncorrelated) residuals (the “ignorance”) using variation
553
partitioning techniques (
). Figure
shows that [a + b] is the
554
relative influence of environmental variables to the model prediction while [b + c]
555
is the relative influence of spatial covariates. The component [b] is the shared vari-
556
ation of [a] and [c] because environmental covariates are spatially structured. The
557
remaining component [d] computed by 1 - [a + b + c] is the residual fraction of
558
the variation. Another benefit of this approach is to have spatial surrogate covari-
559
ates with little or no correlation with the meaningful environmental covariates. This
560
approach has not yet been tested in DSM, but it would certainly make a valuable
561
contribution to increase the interpretability of the ML models and their account of
562
the spatial autocorrelation contained in soil data.
563
3.4. Multivariate mapping
564
Several authors (e.g.
565
) have shown that it is possible to calibrate a single ML model
566
20
to predict either multiple soil properties or a single soil property at multiple depths.
567
This reduces the risk of overfitting, computational resources that would be other-
568
wise required to calibrate several disjoint models (
), and increases
569
prediction accuracy if there is correlation between the variables to predict.
570
) use a multivariate CNN model to predict SOC at multiple soil depths
571
and report a significant increase of prediction accuracy for the deeper soil depths,
572
compared to predictions made for each depth separately by a cubist model.
573
) have shown that for a NN model, it was feasible to constrain the prediction
574
to avoid inconsistent prediction between compositional soil properties, in particular
575
soil texture. It was done by adding an additional layer to the model, but we spec-
576
ulate that this could also be realized by modifying the objective function used to
577
calibrate the model. Despite a few recent studies, there has been little interest in
578
multivariate soil mapping using ML algorithms. In the ML literature, it appears
579
that almost all conventional ML algorithms have a multivariate counterpart. Multi-
580
variate NNs have already been tested in soil mapping studies. An adaptation of the
581
RF algorithm for multivariate mapping is proposed by
) but has
582
several limitations. For example, the calibrated model size increases dramatically
583
when the number of soil properties to predict also increases and it does not allow to
584
separate the contribution of the covariates to each predicted property separately. A
585
theoretical framework for multivariate RF is described by
) and
586
was further implemented in the R language by
). For support
587
vector machines, a multivariate extension is described in
588
589
One objective when mapping soil properties or classes is to learn from the cal-
590
ibrated model. A calibrated multivariate model can provide insights on the soil
591
property and horizon interrelations. Regrettably, in a multivariate machine learning
592
model, the correlation between soil properties or depths is not modelled explicitly
593
(e.g. using a cross-covariance matrix between soil properties). As a result, the corre-
594
lation between properties or depths cannot be assessed internally and no pedological
595
interpretation can be derived from the calibrated model. More research is needed on
596
whether the correlation between original and predicted soil properties (or depths) is
597
preserved in a multivariate ML model. To model the correlation between properties
598
explicitly, two solutions are possible. The first is to calibrate additional stochastic
599
parameters together with the ML parameters (e.g. in a neural network algorithm).
600
This can take the form of an auto-regressive model between the predictions (
601
). Another straightforward solution is to calibrate the model with a crite-
602
rion related to the absolute difference in correlation between the measured properties
603
and predicted properties. While this is easy to implement in ML calibration based
604
21
on an objective function (e.g. neural network), this is not straightforward for models
605
such as RF. Overall, including correlation between properties or depths when predict-
606
ing with a ML algorithm requires further investigation so as to build pedologically
607
realistic and interpretable models.
608
3.5. Uncertainty analysis
609
Uncertainty analysis in digital soil mapping is crucial to deciding whether the
610
predicted soil map is reliable to be used for agricultural production systems or de-
611
cision making. Uncertainty analysis is also about knowing better the limits of the
612
models and is therefore one step towards model interpretability. At the higher level,
613
the machine learning literature distinguishes two sources of uncertainty: aleatoric
614
and epistemic uncertainties (Fig.
). Aleatoric uncertainty is the data noise variance
615
(in other terms, the data error), and arises from noise in the data and measurement
616
error. Epistemic uncertainty refers to model and model parameter uncertainty and
617
represents our ignorance about a true model that generated the data. While epis-
618
temic uncertainty is easy to reduce (e.g. by collecting more data at areas of low
619
sampling density), aleatoric uncertainty is rather difficult to assess (one must repeat
620
the measurement several times) and even more to reduce. Methods to quantify epis-
621
temic uncertainty are bootstrapping, or Bayesian modelling. Quantifying epistemic
622
uncertainty enables to obtain confidence intervals of the prediction. Aleatoric uncer-
623
tainty is mainly quantified by quantile regression methods, but Monte-Carlo simu-
624
lation from the probability distribution of the observations might also be a possible
625
approach. The quantification of both aleatoric and epistemic uncertainty provides
626
prediction intervals with methods such as quantile regression forest (QRF), the Delta
627
or Bayesian methods and the mean plus variance estimate (MVE) for neural network
628
algorithms.
629
630
The recent development of conditional generative adversarial networks (cGAN)
631
) to generate possible realizations of the observations with
632
specific conditions or characteristics seem to be of particular interest to include mea-
633
surement error in DSM. Including measurement error is considered by
634
) for mapping soil organic carbon using uncertain measurement of the soil
635
property. However, the authors do not propose a method to quantify the uncer-
636
tainty of the measurements, nor propagate the measurement error to the predicted
637
map. With cGAN, a probability distribution of the observations is built, which might
638
be used for Monte Carlo simulations. Each Monte Carlo sample is used as input in
639
the ML algorithm, and the final map is the integration of all these simulations. This
640
would effectively tackle the aleatoric uncertainty of the ML model. More impor-
641
22
tantly, this would also quantify the uncertainty present in the measurements, which
642
is currently one of the most important challenges in DSM.
643
644
0
2
4
6
8
10
−
10
−
5
0
5
10
15
x
f(x)
aleatoric
epistemic
Figure 2: Transect with location of the sampling units in red, the true (solid line) and predicted
(dash line) value of the variable of interest, the aleatoric uncertainty (grey shade) and epistemic
uncertainty (blue shape). When no observations are present, the epistemic uncertainty increases.
The aleatoric uncertainty remains somewhat constant across the transect.
Most studies to date do not provide estimate of the uncertainty (Table
). Suc-
645
cessful attempts have been made by
) and
646
to report prediction intervals for random forest and neural networks models, respec-
647
tively. Confidence intervals are reported is several studies (e.g.
648
) and are obtained by training multiple disjoint models
649
using bootstrapped samples of the original data. In a few studies, the variance ob-
650
tained by bootstrapping is averaged by kriging of the residuals (
651
). From Fig.
it follows that if sampling units are selected from a small area
652
in the feature or geographic space, then there will be little uncertainty in this area.
653
Likewise the uncertainty dramatically increases when areas of the feature space are
654
under-sampled, or even worse, ignored. When sampling units are clustered, (spatial)
655
cross-validation might not be sufficient to define realistic prediction accuracy mea-
656
sures because the sampling units used for validation are taken from similar regions
657
of the feature space while the model is biased towards these same regions (
658
). While the (spatial) cross-validation results might show strong agreement
659
between predicted and measured soil property or class and therefore validate a ML
660
23
model with very high predictive abilities, an uncertainty quantification would show
661
unrealistic predictions characterized by a large uncertainty (see right-hand side of
662
Fig.
). This is the results of ML algorithm being very poor predictors for extrapo-
663
lating to areas of the covariate space that is not comprised in the calibration sample.
664
Uncertainty quantification that separates out data and model uncertainties is thus
665
recommended to complete the evaluation of the predicted maps.
666
667
We derive a complementary note about the generation of digital soil maps with
668
ML by the private sector. Companies and commercial software usually do not report
669
measures of the uncertainty associated to the maps and there is no transparency
670
requirement on the methods and quality of soil data. Reporting the uncertainty
671
associated to the prediction is essential to guide decision-making and political action.
672
The danger comes from the generated map which gives the appearance of scientific
673
knowledge where there is none. Making a decision made on maps which are presumed
674
correct but are in fact away from reality, is presumably worse than making a decision
675
made in full appreciation of the limits of the map.
676
3.6. Validation
677
Studies by
) and
) have found that the
678
estimated performance of the machine learning algorithms applied to spatial data
679
depends on the validation strategy. In DSM, model performance is usually assessed
680
using random
k
-fold cross-validation (CV) or single random split of a sample into
681
calibration and validation and/or test subsamples. These strategies give considerably
682
over-optimistic validation statistics estimates because of the presence of autocorre-
683
lation in the observations (
684
). Validation statistics estimated from a random split of the master sample as-
685
sess the ability of the model to reproduce the calibration sample but fail to assess the
686
model performance in terms of spatial mapping (
). As an alterna-
687
tive, several methods (
688
) for spatial cross-validation are proposed to account for spatial
689
autocorrelation of the observations. Two main strategies are adopted.
690
) use a spatial block approach for
k
-fold
691
CV where the master sample is divided into
k
spatially disjoint subsamples using
692
clustering algorithms on the coordinates or by dividing the spatial domain based on
693
k
cells. In
) and
), observations from
694
the calibration subsample that are within a given geographic distance of the valida-
695
tion subsample are omitted from the calibration subsample, after which the model is
696
fitted using the remaining observations from the calibration subsample. While these
697
24
two approaches account for spatial autocorrelation of the observation during valida-
698
tion, further research is required to provide guideline to select the realistic distance
699
from which a validation data point is statistically independent from the calibration
700
sample so as to avoid the opposite effect, i.e. extrapolation and subsequent underop-
701
timistic validation statistics estimates. Spatial-cross validation is integrated in the
702
framework presented in Fig.
703
704
Research on spatial cross-validation has drawn attention to the role of autocor-
705
relation on the calibration on the machine learning algorithms.
706
show that hyperparameter tuning is also impacted by spatial autocorrelation, and
707
that overoptimistic results are reported when the same data are used for performance
708
assessment and parameter tuning. They proposed a nested (block) cross-validation
709
approach for hyperparameter tuning (
) where spatial block are
710
split a second time into spatially disjoint geographic subsamples used to optimize the
711
hyperparameters. The major disadvantage of this method is the dramatic increase in
712
computing time, which is solved by distributed (parallel) computing solutions. Simi-
713
larly to the hyperparameter tuning using nested spatial cross-validation,
714
) showed that autocorrelated covariates lead to overfitting and visible artefacts
715
in the predicted map. The study proposes an iterative procedure for variable se-
716
lection where a group of two variables is first selected based on the error computed
717
with spatial cross-validation, and new variables are iteratively added only if these
718
increase the model performance. The study of
) gives another
719
argument against the use of covariates describing the spatial dependency as these
720
lead to misinterpretation of the model’s important contributors and impossibility for
721
the model to generalize.
722
723
) emphasize the value of visual examination of the predicted
724
maps in addition to the statistical validation. In
), two maps with
725
similar map validation accuracy statistics have a different spatial pattern. The study
726
shows that this is due to the selected covariates, some having strong spatial autocor-
727
relation leading to visible artefacts in the predicted map. This highlights the need
728
for research on the evaluation of predicted maps in terms of spatial pattern.
729
) compare the spatial structure of predicted versus observed values by
730
computing the area under the curve of variograms fitted on the validation locations
731
for both predicted and observed probability of having a peat soil. This relies, how-
732
ever, on the assumption that the variogram of the validation locations represent the
733
mapped area. More research in this direction will be valuable for future DSM stud-
734
ies. To date, visual assessment of the map to detect artefacts, and in consideration
735
25
Does a sampl
e
already exist
?
yes
(legacy data)
no
categorical
continuous
correct for class
imbalanc
e
de-cluster in
geographic or
feature space
Feature space
coverage sampl
e
Does an already
calibrated M
L
algorithm exis
t?
yes
weight by
covariate
importance
Step 2. [b + c]:
y
= f(
W
) +
ε
X
map
of correlated
residuals
scorpan
covariates
W
spatial
covariates
Step 1. [a + b]:
y
= f(
X
) +
ε
uncorrelated
residuals
covariate se
lection us
ing
spatial
cross-validation
covariate se
lection us
ing
spatial
cross-validation on t
he
residuals from Step 1
Step 3. [a + b + c]:
y
= f(
X
W
) +
ε
parameter tuning usi
ng spatial
cross-validation
uncertainty
quantificati
on
prediction error
variance map
prediction
map
spatial
cross
validati
on
+
visual
assess
ement
model-agnostic
interpretation
model-dependent
interpretation
conceptual
V
enn
diagram
[a]
[b]
[c]
unexplained va
riation = [d]
= variation
in
y
validati
on of
uncertainty
estima
tes
Figure
3:
The
recommended
framew
ork
for
digital
soil
mappin
g
with
mac
hine
learning.
The
mo
deller
m
ust
first
decide
whether
a
legacy
soil
sample
or
a
new
sample
is
collected.
He
m
ust
also
decide
whether
the
ob
jectiv
e
is
a
categorical
or
quan
titativ
e
map.
The
recommended
framew
ork
enables
the
separation
b
et
w
een
the
v
ariation
explained
b
y
the
p
edologically
relev
an
t
co
v
ariates
and
b
y
the
spatial
co
v
ariates.
It
is
recommended
to
use
a
spatial
cross-v
alidation
strategy
for
v
alidation,
but
also
for
parameter
tuning
and
co
v
ariate
selection.
26
of our knowledge of soil forming processes, is the best option.
736
3.7. Machine learning and pedological knowledge
737
Accounting for existing expert soil knowledge in DSM with machine learning is a
738
challenging exercise (
). ML algorithms do not build on any existing a
739
priori conceptual model of the soil processes and only processes that are conveyed by
740
the input data are represented in the map (
741
To prevent extrapolation,
) do not provide soil maps in some under-
742
sampled areas of the globe such as deserts and glaciers for global mapping of several
743
soil properties. This stems from incomplete datasets of soil observations for these
744
areas, despite that extensive expert knowledge exists. In
) this is
745
solved by integrating the expert knowledge in the form of expert-based pseudo-points
746
to guide the ML model in areas of evident extrapolation. In
), 600
747
pseudo-points are also added in under-represented areas of the geographic space. The
748
study stresses the importance of consulting an expert when building a ML model.
749
In the same study, meaningful covariates are selected based on existing knowledge
750
on the soil process, and plausibility of the predicted soil map is made in consid-
751
eration of the knowledge of soil forming process. On many occasions, meaningful
752
covariates are selected for mapping soil properties or classes. For example,
753
) used a set of covariates selected a priori by an expert on the area under
754
study. In
) a set of
scorpan
covariates is selected for
755
mapping soil properties in Australia. These examples show that in the literature,
756
adding expert-based pseudo-points and selecting meaningful covariates are, to date,
757
two straightforward options to include existing knowledge into a ML algorithm for
758
DSM.
759
760
The above shows that little is known on how to account for existing knowledge
761
in ML models. Unfortunately, this is the same order in which the complexity of the
762
models increases and our understanding of the model functioning decreases. The
763
increasing caution in the use of predictions made by a complex ML model that one
764
should expect as a result is not evident. A ML model predicting a number based on
765
relationships between covariates that are unknown in the view of existing knowledge,
766
should not be taken with the same seriousness as a number predicted by mechanistic
767
steps or an established theory. Improvement in this situation is made by ensur-
768
ing that the calibrated ML algorithm matches the existing knowledge of the soil
769
processes, for example by reflecting or confirming the current hypothesis or prior
770
knowledge on the soil spatial variation for an area. If the model prediction does not
771
agree with existing maps, this means that the model has instead modelled a different
772
27
process and is thus likely to be invalid. A model is invalid until it is validated, not
773
only against data, but also against the researcher experience and validation of the
774
model creation process (
). In short, pedological knowledge should be
775
integrated to enforce results consistent with the existing scientific principles. This
776
can be done at each step of the model building, calibration and validation. One
777
can incorporate additional knowledge by selecting appropriate covariates or adding
778
pseudo-points. In model building, knowledge takes the form of a hybrid model, a spe-
779
cific model architecture or objective function (in neural networks models) constrain-
780
ing the calibration process according to specific knowledge. For example,
781
) adds the constraint that the prediction of topsoil clay, silt and sand must
782
sum to 100% in a neural network model. Finally, pedological knowledge is used to
783
make
post-hoc
checks on the plausibility of the calibrated model and predicted maps.
784
785
) stress that since ML models (the author used the term “predic-
786
tive process model” in the sense in which “machine learning models” is used in this
787
article) have no connection to established theory, one can never be sure that the
788
outcome is realistic given the real-world processes involved. The problem is that a
789
non-valid model is difficult to recognize and to reject since it is often not interpretable
790
by a human. To ensure that models fit the existing knowledge, they must be opened
791
and understood in their functioning. Opening the “black box” is then necessary but
792
not straightforward (see next section on interpretability), and is often reduced to the
793
analysis of which environmental covariates are the most often used by the model to
794
make a prediction (see for example
) or
795
796
797
Several authors, however, have warned against the use of accuracy metrics for
798
pedological interpretation (e.g.
) or
799
) use meaningless, pseudo-covariates to map soil organic car-
800
bon over a hypothetical area. The authors obtain an accurate map, and conclude
801
that ML algorithm should not be used for obtaining new soil knowledge because the
802
ML algorithm aims at predicting a pattern rather than finding causal relationships.
803
) suggest to use calibrated ML models as a “hypothesis dis-
804
covery” tool, in which the mechanisms conveyed by the calibrated ML model are
805
supplied to the researcher for possible explanations of the soil process, which can
806
then be confronted to experiments and principles of soil genesis. The challenge that
807
then arises, noticed by
) is the conversion of the mechanisms of the
808
ML model (the model “knowledge”) from a data language to a human one. The data
809
language is typically parameters or metrics such as the “mean decrease of purity”
810
28
or “Gini importance index” of a covariate to assess its importance in the prediction
811
of a soil property or class. Such metrics are not interpretable in terms of human
812
explanation and they do not relate to soil processes. Translating the data language
813
to the domain (the human language) requires some attention and further research.
814
More discussion on this issue is found in
).
815
3.8. Interpretation of the models
816
Soil scientists rely on ML algorithms to gain insights into the modelled processes.
817
Despite providing higher prediction accuracy than other conventional models, ML
818
models are considered as a black box. Broadly speaking, we do not learn from the
819
model how the input covariates are related to the output soil property or classes,
820
and what are the underlying mechanisms behind the prediction. This is unfortunate
821
for soil science because in many cases the model itself is considered as a source of
822
knowledge in addition to the collected soil data. Scientific findings remain hidden
823
when the model only gives a prediction without explanations. In this case, the inter-
824
pretability of the model warrants the extraction of the knowledge captured by the
825
calibrated model.
) defines intepretability as the degree to which human
826
can understand the cause of a decision. In general, the need for interpretability of a
827
machine learning algorithm stems from a deficiency in problem formalization (
828
). This means that for a given task (i.e. mapping
829
the spatial distribution of soil organic carbon), the prediction itself does not fully
830
solve the original problem. We suggest three reasons which drive the demand for
831
interpretability in DSM (adapted from
)). The first and
832
most obvious reason is to increase our scientific understanding of the soil system by
833
extracting knowledge from the mechanisms captured by the model. Scientists wish
834
to know which are the drivers of a soil process and, more importantly, whether the
835
mechanisms captured by the model confirm our scientific understanding of the sys-
836
tem (see Section
). The second reason is to audit the calibrated ML algorithm. Is
837
the ML algorithm predicting for the right reasons? If a scientist makes a model for
838
mapping the topsoil nitrogen content of a field, the interpretation might reveal that
839
the model is actually predicting soil clay, that is, a proxy of the initial objective.
840
The third reason is to avoid financial loss or to prevent a safety issue. Take the
841
example of the remediation of the soil due to radioactive fallout after the Fukushima
842
nuclear accident. A map of contaminated soils made by a ML algorithm would typi-
843
cally predict the dominant soil type characteristics, i.e. forest soil (about 75% of the
844
area), for classification into contaminated or not contaminated areas. Interpretation
845
of the model might then reveal that the important features learned by the model are
846
unrealistic for agricultural landscapes and residential areas whose remediation is yet
847
29
?
Nature
Black-box
model
Data
Interpretation
Humans
inform
extract
calibrate
sample
Figure 4: Summary framework for model-agnostic interpretable ML, adapted from
The lowest level is the reality, the unknown real-world soil that one wants to predict. The sec-
ond level is the dataset that is extracted from the reality. We collect a fraction of the reality, a
sample, and link it to exhaustively known environmental covariates. The relationships between the
covariates and the sample is learned by a black-box machine learning model (level 3), on top of
which comes the interpretation level to extract some knowledge from the structure of the calibrated
model. The structure of the model is converted to human understandable knowledge.
31
A straightforward way to increase interpretability is to decrease model complexity,
850
for example by building a single decision tree instead of a random forest composed
851
of several thousand trees. A simple model enables visualization of the important
852
mechanisms of the model and resultant explanations. For DT algorithms, it is possi-
853
ble to map the predicted values for specific rules (if the model is sufficiently simple).
854
Decreasing complexity, however, is done at the expense of model prediction accu-
855
racy. For more complex ML algorithms, built-in features allow the user to retrieve
856
the variable importance. In decision-tree like algorithms, the variable importance
857
is derived from the thresholds used for the splits. For neural networks, the output
858
weights associated with the input layer neurons provides an indication of the impor-
859
tant features (
). One drawback of these techniques is their inability to
860
provide information on whether the covariates have a causal link to the modelled soil
861
property or class, which leads several authors to warn against their use for knowledge
862
discovery (e.g.
). More importantly, these variable importance
863
metrics are summary statistics not always meaningful and they are model-specific,
864
i.e. they preclude comparison between models or parts of the predicted map.
865
) reviews techniques to interpret ML algorithms and define two main categories
866
of interpretation techniques. The first are the model specific ones. They are rou-
867
tinely used in DSM activities (e.g. RF variable importance). The second category
868
falls into model-free techniques, also called model-agnostics (
). It en-
869
ables the users to use any model, thus not restraining themselves to simple models or
870
models with embedded features of interpretation. A summary of how model-agnostic
871
techniques are employed is shown in Fig.
. Examples of model-agnostic techniques
872
are the partial dependence plot (
) if the number of covariate is small
873
(two maximum), individual conditional expectation (
), and
874
global or local model-agnostic explanation (LIME,
). Finally,
875
sensitivity analysis is also a straightforward means of
post-hoc
interpretation of how
876
the model output depends upon the different covariates.
877
4. The way forward
878
Machine learning algorithms are now extensively used in soil mapping for re-
879
gression and classification purposes, much in the same way as routinely employed
880
in other fields of science. There is no doubt that prediction accuracy benefits from
881
these data-driven models because ML algorithms are not constrained by a pre-defined
882
conceptual model of the soil spatial variation, in comparison to mechanistic or even
883
geostatistical models. The question now is how to increase our scientific understand-
884
ing of the soil and how to adapt and guide the use of ML to the challenges pertaining
885
32
to soil mapping and soil science in general. Future research on soil mapping with
886
machine learning should incorporate the three core elements proposed by
887
) and
) which we adapted, as follows:
888
889
Plausibility
: Models should not only be accurate but also valid in light of the
890
current knowledge and scientific theories. A model should predict for the right rea-
891
sons. The plausibility is the solution path taken by the ML algorithm to link the
892
input to the output, and does not depend directly on the data (
). In
893
practical terms, it starts with the model building step, by feeding the model with
894
credible covariates and by accounting for the spatial particularities of soil data. Spa-
895
tial or temporal correlation among data should be modelled, either by using a specific
896
model (e.g. a convolutional neural network), or by using a model architecture that
897
accounts for this particularity (see Section
). Plausibility also takes the form of
898
model constraints, to avoid the prediction of unrealistic proportions or ratios. The
899
plausibility can be further tested in terms of model simulatability (
900
Since ML algorithms can model arbitrary patterns, there should be some attempts
901
to test the model with synthetic data or data from a calibrated mechanistic model
902
representing a large range of dynamics (
). Increasing model
903
plausibility will facilitate the acceptance of ML to a large range of scenarios in soil
904
science.
905
906
Interpretability
: Interpretability is the translation of an abstract model or model
907
output into terms understandable by humans (
). Model in-
908
terpretability pairs with model plausibility and hypothesis discovery. Complex and
909
arbitrary patterns extracted from the data by an algorithm can be understood only
910
by the transparency of the model. Interpretation is obtained by model-specific and
911
model-agnostic methods, described in Section
. Visual examination of the maps is
912
also a means of interpretability. While complex ML models are potentially harmful
913
because they often do not model any real-world process, there is an opportunity to
914
challenge existing knowledge by
post-hoc
comparison of existing maps produced by
915
expert knowledge with the maps predicted by a ML model, and by analysis of the
916
striking differences. This is possible only if the model is interpretable by humans
917
and the physical relationships between variables are realistic (the model is
plausible
).
918
Model interpretation is also an opportunity to generate new hypotheses, by inter-
919
preting the relationships found by the ML algorithm in the stores of soil data. The
920
new hypotheses derived by these interpretations may challenge existing knowledge
921
on the soil spatial variation and genesis.
922
923
33
Explainability
: Modellers should shy away from mindless model fitting and pre-
924
diction and intensify research on models that both predict and explain. Explanations
925
aim to answer the three questions:
what
is the modelled process?,
how
has it been
926
modelled?, and
why
has this process been modelled? (
). In this sense,
927
explaining a process is an interpretation of a ML model plus expert knowledge and
928
contextual information. For example, a different explanation is warranted when one
929
wants to explain the pattern of a predicted soil map or the reason for two close
930
predicted soil classes to be different. To explain, the modeller uses the data, the
931
plausibility of the model and its interpretation using expert knowledge (see Fig.
932
Explainability is helped by model structure providing algorithmic explanations in
933
the form of graphs or equations.
934
935
An example of model structure providing algorithmic explanations in DSM is
936
found with the use of Bayesian belief networks (BBN,
) in
937
) and later
). BBN is a probabilistic graphical model pre-
938
dicting the likely value of a soil property or class given conditional dependencies
939
between covariates. Recent advances in ML have made a step further by discovering
940
the graph structure directly from the data. However, while BNN is an interpretable
941
ML model of conditional dependence between variables, the process that generated
942
these dependencies remains hidden. To discover new processes from data, inductive
943
process modelling (
944
) are the way forward. Both are automated model discovery process,
945
in which equations describing a process are inductively (i.e. using the data) assem-
946
bled into a single predictive model by heuristic search methods (
947
). The calibrated model is a set of equations constrained by existing,
948
verified equations (e.g. differential equations of the water flow) representing causal
949
relationships between variables. The model can be refined using expert knowledge
950
and additional data (
). More importantly, these models produce
951
explanations, which can be refuted or approved in light of scientific principles.
952
34
input data
output
model
pedological
knowledge
plausibility
interpretability
explainability
scientific
outcome
Figure 5: Conceptual framework for the derivation of a scientific outcome from a ML model, adapted
from
). The light grey box represents the conventional use of ML algorithms
in digital soil mapping, in which an output is derived from a calibrated ML model given a set of
input data. A scientific outcome is obtained by explaining the output of a model using pedological
knowledge, but also by ensuring scientific consistency at each link of the chain. Alternatively, a
plausible and interpretable model can be explained using pedological knowledge.
Figure
illustrates the central role played by the three elements
plausibility
,
953
interpretability
and
explainability
in obtaining a scientific outcome from machine
954
learning. Fig.
shows that the three core elements are conditioned to the use of
955
pedological knowledge at each link of the chain. Enforcing pedological knowledge
956
during modelling restricts the solution space to scientifically consistent results and
957
may decrease the overall prediction accuracy. For digital soil mapping purposes,
958
it is not obvious whether an increase of predictive accuracy worth the substantial
959
decrease of model consistency. For this reason, recent studies (e.g.
960
) advocate the use of other criteria to measure the
961
overall performance, such as model complexity or consistency (
962
Including other criteria to assess the overall performance of a ML model would
963
certainly make one step towards “conscious” digital soil mapping, and participate to
964
the uptake of knowledge discovery via machine learning in soil science.
965
35
5. Conclusion
966
In this contribution, we have reviewed the current and prospective use of ML al-
967
gorithms for digital soil mapping. From the existing use of ML in DSM, we identified
968
key challenges and provided partial solutions We draw the following conclusions.
969
•
There has been a large number of studies mapping soil properties or classes
970
using ML algorithms. A wide range of soil properties, attributes and types have
971
been predicted. Likewise, an increasing number of machine learning algorithms
972
have been tested. Case studies are dominated by the use of legacy samples for
973
local to regional scale (about 10
4
km
2
) areas. Ensemble of different algorithms
974
to improve prediction are gaining more attention. All studies reported at least
975
one validation statistics but few reported the uncertainty associated with the
976
prediction.
977
•
The configuration of a good sampling design for mapping with machine learning
978
is largely unknown. The impact of the sampling design on model calibration
979
and prediction has generally been disregarded. More research is needed in this
980
direction.
981
•
A large number of studies have focused solely on achieving a high mapping
982
accuracy. Comparison between models and other studies are made based on
983
validation statistics, while ignoring model complexity or consistency with re-
984
spect to the existing pedological knowledge.
985
•
The benefit of using a large number of covariates, or pseudo-covariates ac-
986
counting for residual spatial autocorrelation for mapping using ML algorithms
987
should be avoided. To build consistent models, we suggested to select a set of
988
pedologically relevant covariates, and to model the potential residual spatial
989
autocorrelation with
post-hoc
fitting of another model using spatial surrogate
990
covariates. This procedure also enables a separate analysis of the variation
991
explained by environmental or spatial covariates.
992
Overall, our review of the literature suggested that in recent studies inference is
993
relegated to the background with the emergence of the mapping accuracy as the sole
994
standard by which progress is measured. While the mapping accuracy is valuable, it
995
should not be the only objective one should pursue. To date, ML is applied to digital
996
soil mapping the same way as other fields such as image detection or pattern recog-
997
nition do. Any prediction can become a soil map, whether it contains soil knowledge
998
or not, and without any assessment on whether the fitted relationships relate to a
999
36
real-world soil process.
1000
1001
We also found, however, that there is opportunity to include pedological knowl-
1002
edge at each step of the modelling chain, to improve or correct the existing dataset,
1003
to design the model architecture, to constrain the model calibration, or to analyse
1004
the output using
post-hoc
checks on the predicted soil maps. Future studies on DSM
1005
should use
plausible
,
interpretable
and
explainable
ML models to extract novel sci-
1006
entific results from soil data. One step towards achieving this goal is to integrate
1007
model consistency in addition to model prediction accuracy to evaluate the overall
1008
performance of the mapping approach. This will ensure that future studies use mod-
1009
els that are not only accurate but also valid in light of the current knowledge and
1010
scientific theories.
1011
Acknowledgement
1012
Budiman Minasny is member of a consortium supported by LE STUDIUM Loire
1013
Valley Institute for Advanced Studies through its LE STUDIUM Research Consor-
1014
tium Programme.
1015
References
1016
Adhikari, K., Hartemink, A. E., Minasny, B., Kheir, R. B., Greve, M. B., & Greve,
1017
M. H. (2014). Digital mapping of soil organic carbon contents and stocks in den-
1018
mark.
PLOS ONE
,
9
, e105519.
1019
Aitkenhead, M. J., & Coull, M. C. (2016). Mapping soil carbon stocks across Scotland
1020
using a neural network model.
Geoderma
,
262
, 187–198.
1021
Akpa, S. I. C., Odeh, I. O. A., Bishop, T. F. A., & Hartemink, A. E. (2014). Digital
1022
mapping of soil particle-size fractions for Nigeria.
Soil Science Society of America
1023
Journal
,
78
, 1953–1966.
1024
Arrouays, D., McKenzie, N., Hempel, J., de Forges, A. R., & McBratney, A. B.
1025
(2014).
GlobalSoilMap: Basis of the Global Spatial Soil Information System
. CRC
1026
press, Boca Raton, USA.
1027
Asgharbeygi, N., Langley, P., Bay, S., & Arrigo, K. (2006). Inductive revision of
1028
quantitative process models.
Ecological Modelling
,
194
, 70–79.
1029
37
Batjes, N. H., Ribeiro, E., van Oostrum, A., Leenaars, J., Hengl, T., & Mendes de
1030
Jesus, J. (2017). WoSIS: providing standardised soil profile data for the world.
1031
Earth System Science Data
,
9
, 1–14.
1032
Beguin, J., Fuglstad, G.-A., Mansuy, N., & Par´
e, D. (2017). Predicting soil properties
1033
in the Canadian boreal forest with limited data: Comparison of spatial and non-
1034
spatial statistical approaches.
Geoderma
,
306
, 195–205.
1035
Behrens, T., F¨
orster, H., Scholten, T., Steinr¨
ucken, U., Spies, E.-D., & Goldschmitt,
1036
M. (2005). Digital soil mapping using artificial neural networks.
Journal of Plant
1037
Nutrition and Soil Science
,
168
, 21–33.
1038
Behrens, T., Schmidt, K., MacMillan, R. A., & Viscarra-Rossel, R. A. (2018a).
1039
Multi-scale digital soil mapping with deep learning.
Scientific Reports
,
8
, 15244.
1040
Behrens, T., Schmidt, K., Viscarra-Rossel, R. A., Gries, P., Scholten, T., & MacMil-
1041
lan, R. A. (2018b). Spatial modelling with Euclidean distance fields and machine
1042
learning.
European Journal of Soil Science
,
69
, 757–770.
1043
Behrens, T., Zhu, A.-X., Schmidt, K., & Scholten, T. (2010). Multi-scale digital
1044
terrain analysis and feature selection for digital soil mapping.
Geoderma
,
155
,
1045
175–185.
1046
Bel, L., Allard, D., Laurent, J. M., Cheddadi, R., & Bar-Hen, A. (2009). CART al-
1047
gorithm for spatial data: Application to environmental and ecological data.
Com-
1048
putational Statistics & Data Analysis
,
53
, 3082–3093.
1049
Bel, L., Laurent, J. M., Bar-Hen, A., Allard, D., & Cheddadi, R. (2005). A spatial
1050
extension of CART: application to classification of ecological data. In P. Renard,
1051
H. Demougeot-Renard, & R. Froidevaux (Eds.),
Geostatistics for Environmental
1052
Applications
(pp. 99–109). Springer, Berlin, Heidelberg.
1053
Bennett, N. D., Croke, B. F. W., Guariso, G., Guillaume, J. H. A., Hamilton, S. H.,
1054
Jakeman, A. J., Marsili-Libelli, S., Newham, L. T. H., Norton, J. P., Perrin, C.
1055
et al. (2013). Characterising performance of environmental models.
Environmental
1056
Modelling & Software
,
40
, 1–20.
1057
Blanco, C. M. G., Gomez, V. M. B., Crespo, P., & Ließ, M. (2018). Spatial prediction
1058
of soil water retention in a P´
aramo landscape: Methodological insight into machine
1059
learning using random forest.
Geoderma
,
316
, 100–114.
1060
Braben, D. W. (1985). Innovation and academic research.
Nature
,
316
, 401–402.
1061
38
Breiman, L. (2017).
Classification and Regression Trees
. Routledge, New York, USA.
1062
Brenning, A. (2012). Spatial cross-validation and bootstrap for the assessment of
1063
prediction rules in remote sensing: The R package sperrorest. In
2012 International
1064
Geoscience and Remote Sensing Symposium
(pp. 5372–5375). IEEE.
1065
Bridewell, W., Langley, P., Todorovski, L., & Dˇ
zeroski, S. (2008). Inductive process
1066
modeling.
Machine learning
,
71
, 1–32.
1067
Brungard, C. W., Boettinger, J. L., Duniway, M. C., Wills, S. A., & Edwards Jr, T. C.
1068
(2015). Machine learning for predicting soil classes in three semi-arid landscapes.
1069
Geoderma
,
239
, 68–83.
1070
Brus, D. J. (2019). Sampling for digital soil mapping: A tutorial supported by R
1071
scripts.
Geoderma
,
338
, 464–480.
1072
Bui, E., Henderson, B., & Viergever, K. (2009). Using knowledge discovery with data
1073
mining from the Australian Soil Resource Information System database to inform
1074
soil carbon mapping in Australia.
Global Biogeochemical Cycles
,
23
, GB4033.
1075
Bui, E. N., Loughhead, A., & Corner, R. (1999). Extracting soil-landscape rules
1076
from previous soil surveys.
Soil Research
,
37
, 495–508.
1077
Carr´
e, F., McBratney, A. B., & Minasny, B. (2007). Estimation and potential im-
1078
provement of the quality of legacy soil samples for digital soil mapping.
Geoderma
,
1079
141
, 1–14.
1080
Chen, S., Mulder, V. L., Martin, M. P., Walter, C., Lacoste, M., Richer-de Forges,
1081
A. C., Saby, N. P. A., Loiseau, T., Hu, B., & Arrouays, D. (2019). Probability
1082
mapping of soil thickness by random survival forest at a national scale.
Geoderma
,
1083
344
, 184–194.
1084
Cialella, A. T., Dubayah, R., Lawrence, W., & Levine, E. (1997). Predicting soil
1085
drainage class using remotely sensed and digital elevation data.
Photogrammetric
1086
Engineering and Remote Sensing
,
63
, 171–177.
1087
Cooper, G. F. (1990). The computational complexity of probabilistic inference using
1088
bayesian belief networks.
Artificial intelligence
,
42
, 393–405.
1089
Coveney, P. V., Dougherty, E. R., & Highfield, R. R. (2016). Big data need big theory
1090
too.
Philosophical Transactions of the Royal Society A: Mathematical, Physical and
1091
Engineering Sciences
,
374
, 20160153.
1092
39
Cressie, N., & Johannesson, G. (2008). Fixed rank kriging for very large spatial data
1093
sets.
Journal of the Royal Statistical Society: Series B (Statistical Methodology)
,
1094
70
, 209–226.
1095
Dai, F., Zhou, Q., Lv, Z., Wang, X., & Liu, G. (2014). Spatial prediction of soil
1096
organic matter content integrating artificial neural network and ordinary kriging
1097
in Tibetan Plateau.
Ecological Indicators
,
45
, 184–194.
1098
Dale, M. B., McBratney, A. B., & Russell, J. S. (1989). On the role of expert systems
1099
and numerical taxonomy in soil classification.
Journal of Soil Science
,
40
, 223–234.
1100
De Gruijter, J. J., Brus, D. J., Bierkens, M. F. P., & Knotters, M. (2006).
Sampling
1101
for Natural Resource Monitoring
. Springer Science & Business Media, Dordrecht,
1102
NL.
1103
Dharumarajan, S., Hegde, R., & Singh, S. K. (2017). Spatial prediction of major
1104
soil properties using Random Forest techniques-A case study in semi-arid tropics
1105
of South India.
Geoderma Regional
,
10
, 154–162.
1106
Doshi-Velez, F., & Kim, B. (2017).
Towards a rigorous science of interpretable
1107
machine learning.
1108
Ellili, Y., Malone, B. P., Michot, D., Minasny, B., Vincent, S., Walter, C., &
1109
Lemercier, B. (2019). Comparing three approaches of spatial disaggregation of
1110
legacy soil maps based on DSMART algorithm.
SOIL Discussions
, .
1111
Evrard, O., Laceby, J. P., & Nakao, A. (2019). Effectiveness of landscape decontam-
1112
ination following the Fukushima nuclear accident: a review.
SOIL
,
5
, 333–350.
1113
Fick, S. E., & Hijmans, R. J. (2017). Worldclim 2: new 1-km spatial resolution
1114
climate surfaces for global land areas.
International Journal of Climatology
,
37
,
1115
4302–4315.
1116
Forkuor, G., Hounkpatin, O. K. L., Welp, G., & Thiel, M. (2017). High resolution
1117
mapping of soil properties using remote sensing variables in south-western Burkina
1118
Faso: a comparison of machine learning and multiple linear regression models.
1119
PLOS ONE
,
12
, e0170478.
1120
Fourcade, Y., Besnard, A. G., & Secondi, J. (2018). Paintings predict the distribution
1121
of species, or the challenge of selecting environmental predictors and evaluation
1122
statistics.
Global Ecology and Biogeography
,
27
, 245–256.
1123
40
Friedman, J. H. (2001). Greedy function approximation: a gradient boosting ma-
1124
chine.
Annals of statistics
,
25
, 1189–1232.
1125
Gahegan, M. (2000). On the application of inductive machine learning tools to
1126
geographical analysis.
Geographical Analysis
,
32
, 113–139.
1127
Gahegan, M. (2019). Fourth paradigm GIScience? prospects for automated discovery
1128
and explanation from data.
International Journal of Geographical Information
1129
Science
,
34
, 1–21.
1130
Gahegan, M., Wachowicz, M., Harrower, M., & Rhyne, T.-M. (2001). The integration
1131
of geographic visualization with knowledge discovery in databases and geocompu-
1132
tation.
Cartography and Geographic Information Science
,
28
, 29–44.
1133
Gasch, C. K., Hengl, T., Gr¨
aler, B., Meyer, H., Magney, T. S., & Brown, D. J.
1134
(2015). Spatio-temporal interpolation of soil water, temperature, and electrical
1135
conductivity in 3D+ T: The Cook Agronomy Farm data set.
Spatial Statistics
,
1136
14
, 70–90.
1137
Gascon, F., Bouzinac, C., Th´
epaut, O., Jung, M., Francesconi, B., Louis, J., Lonjou,
1138
V., Lafrance, B., Massera, S., Gaudel-Vacaresse, A. et al. (2017). Copernicus
1139
Sentinel-2A calibration and products validation status.
Remote Sensing
,
9
, 584.
1140
Georganos, S., Grippa, T., Gadiaga, A. N., Linard, C., Lennert, M., Vanhuysse, S.,
1141
Mboga, N. O., Wolff, E., & Kalogirou, S. (2019). Geographical random forests: A
1142
spatial extension of the random forest algorithm to address spatial heterogeneity
1143
in remote sensing and population modelling.
Geocarto International
,
1
, 1–12.
1144
Goldberg, D. E., & Holland, J. H. (1988).
Genetic algorithms and machine learning
.
1145
Kluwer Academic Publishers-Plenum Publishers. Kluwer Academic Publishers.
1146
Goldstein, A., Kapelner, A., Bleich, J., & Pitkin, E. (2015). Peeking inside the black
1147
box: Visualizing statistical learning with plots of individual conditional expecta-
1148
tion.
Journal of Computational and Graphical Statistics
,
24
, 44–65.
1149
Gomes, L. C., Faria, R. M., de Souza, E., Veloso, G. V., Schaefer, C. E. G. R., &
1150
Fernandes Filho, E. I. (2019). Modelling and mapping soil organic carbon stocks
1151
in Brazil.
Geoderma
,
340
, 337–350.
1152
Grimm, R., Behrens, T., M¨
arker, M., & Elsenbeer, H. (2008). Soil organic carbon
1153
concentrations and stocks on Barro Colorado Island—Digital soil mapping using
1154
Random Forests analysis.
Geoderma
,
146
, 102–113.
1155
41
Guevara, M., Olmedo, G. F., Stell, E., Yigini, Y., Aguilar Duarte, Y., Arel-
1156
lano Hern´
andez, C., Ar´
evalo, G. E., Arroyo-Cruz, C. E., Bolivar, A., Bunning,
1157
S. et al. (2018). No silver bullet for digital soil mapping: country-specific soil
1158
organic carbon estimates across Latin America.
SOIL
,
4
, 173–193.
1159
Hamzehpour, N., Shafizadeh-Moghadam, H., & Valavi, R. (2019). Exploring the
1160
driving forces and digital mapping of soil organic carbon using remote sensing and
1161
soil texture.
CATENA
,
182
, 104141.
1162
Hansen, M. K., Brown, D. J., Dennison, P. E., Graves, S. A., & Bricklemyer, R. S.
1163
(2009). Inductively mapping expert-derived soil-landscape units within dambo
1164
wetland catenae using multispectral and topographic data.
Geoderma
,
150
, 72–
1165
84.
1166
H¨
aring, T., Dietz, E., Osenstetter, S., Koschitzki, T., & Schr¨
oder, B. (2012). Spa-
1167
tial disaggregation of complex soil map units: A decision-tree based approach in
1168
Bavarian forest soils.
Geoderma
,
185
, 37–47.
1169
Hartmann, J., & Moosdorf, N. (2012). The new global lithological map database
1170
GLiM: A representation of rock properties at the Earth surface.
Geochemistry,
1171
Geophysics, Geosystems
,
13
, 1–37.
1172
Hawkins, B. A. (2012). Eight (and a half) deadly sins of spatial analysis.
Journal of
1173
Biogeography
,
39
, 1–9.
1174
He, H., & Garcia, E. A. (2009). Learning from imbalanced data.
IEEE Transactions
1175
on Knowledge and Data Engineering
,
21
, 1263–1284.
1176
Henderson, B. L., Bui, E. N., Moran, C. J., & Simon, D. A. P. (2005). Australia-wide
1177
predictions of soil properties using decision trees.
Geoderma
,
124
, 383–398.
1178
Hengl, T., de Jesus, J. M., Heuvelink, G. B. M., Gonzalez, M. R., Kilibarda, M.,
1179
Blagoti´
c, A., Shangguan, W., Wright, M. N., Geng, X., Bauer-Marschallinger, B.
1180
et al. (2017a). SoilGrids250m: Global gridded soil information based on machine
1181
learning.
PLOS ONE
,
12
, e0169748.
1182
Hengl, T., de Jesus, J. M., MacMillan, R. A., Batjes, N. H., Heuvelink, G. B. M.,
1183
Ribeiro, E., Samuel-Rosa, A., Kempen, B., Leenaars, J. G., Walsh, M. G. et al.
1184
(2014).
SoilGrids1km—global soil information based on automated mapping.
1185
PLOS ONE
,
9
, e105992.
1186
42
Hengl, T., Leenaars, J. G. B., Shepherd, K. D., Walsh, M. G., Heuvelink, G. B. M.,
1187
Mamo, T., Tilahun, H., Berkhout, E., Cooper, M., Fegraus, E. et al. (2017b). Soil
1188
nutrient maps of Sub-Saharan Africa: assessment of soil nutrient content at 250
1189
m spatial resolution using machine learning.
Nutrient Cycling in Agroecosystems
,
1190
109
, 77–102.
1191
Hengl, T., Nussbaum, M., Wright, M. N., Heuvelink, G. B. M., & Gr¨
aler, B. (2018).
1192
Random forest as a generic framework for predictive modeling of spatial and spatio-
1193
temporal variables.
PeerJ
,
6
, e5518.
1194
Heung, B., Ho, H. C., Zhang, J., Knudby, A., Bulmer, C. E., & Schmidt, M. G.
1195
(2016). An overview and comparison of machine-learning techniques for classifica-
1196
tion purposes in digital soil mapping.
Geoderma
,
265
, 62–77.
1197
Heuvelink, G. B. M., & Webster, R. (2001). Modelling soil variation: past, present,
1198
and future.
Geoderma
,
100
, 269–301.
1199
Holmes, K. W., Odgers, N. P., Griffin, E. A., & van Gool, D. (2014). Spatial disag-
1200
gregation of conventional soil mapping across Western Australia using DSMART.
1201
In D. Arrouays, N. McKenzie, J. Hempel, A. Richer de Forges, & A. B. McBratney
1202
(Eds.),
GlobalSoilMap: Basis of the Global Spatial Soil Information System
(pp.
1203
273–279). Taylor & Francis, London, UK.
1204
Hounkpatin, K. O. L., Schmidt, K., Stumpf, F., Forkuor, G., Behrens, T., Scholten,
1205
T., Amelung, W., & Welp, G. (2018). Predicting reference soil groups using legacy
1206
data: A data pruning and Random Forest approach for tropical environment (Dano
1207
catchment, Burkina Faso).
Scientific Reports
,
8
, 9959.
1208
Ill´
es, G., Sutikno, S., Szatm´
ari, G., Sandhyavitri, A., P´
asztor, L., Kristijono, A.,
1209
Moln´
ar, G., Yusa, M., & Sz´
ekely, B. (2019). Facing the peat co
2
threat: digital
1210
mapping of Indonesian peatlands—a proposed methodology and its application.
1211
Journal of Soils and Sediments
, (pp. 1–16).
1212
Jenny, H. (1941).
Factors of Soil Formation: A System of Quantitative Pedology
.
1213
McGrawHill, New York, USA.
1214
Jiang, Z., Li, Y., Shekhar, S., Rampi, L., & Knight, J. (2017). Spatial ensemble
1215
learning for heterogeneous geographic data with class ambiguity: A summary of
1216
results. In
Proceedings of the 25th ACM SIGSPATIAL International Conference
1217
on Advances in Geographic Information Systems
23 (pp. 1–10). ACM.
1218
43
Kalambukattu, J. G., Kumar, S., & Raj, R. A. (2018). Digital soil mapping in
1219
a Himalayan watershed using remote sensing and terrain parameters employing
1220
artificial neural network model.
Environmental Earth Sciences
,
77
, 203.
1221
Karpatne, A., Atluri, G., Faghmous, J. H., Steinbach, M., Banerjee, A., Ganguly, A.,
1222
Shekhar, S., Samatova, N., & Kumar, V. (2017). Theory-guided data science: A
1223
new paradigm for scientific discovery from data.
IEEE Transactions on Knowledge
1224
and Data Engineering
,
29
, 2318–2331.
1225
Keskin, H., Grunwald, S., & Harris, W. G. (2019). Digital mapping of soil carbon
1226
fractions with machine learning.
Geoderma
,
339
, 40–58.
1227
Kheir, R. B., Greve, M. H., Abdallah, C., & Dalgaard, T. (2010a). Spatial soil zinc
1228
content distribution from terrain parameters: A GIS-based decision-tree model in
1229
Lebanon.
Environmental Pollution
,
158
, 520–528.
1230
Kheir, R. B., Greve, M. H., Bøcher, P. K., Greve, M. B., Larsen, R., & McCloy,
1231
K. (2010b). Predictive mapping of soil organic carbon in wet cultivated lands
1232
using classification-tree based models: The case study of Denmark.
Journal of
1233
Environmental Management
,
91
, 1150–1160.
1234
Kirkwood, C., Cave, M., Beamish, D., Grebby, S., & Ferreira, A. (2016). A machine
1235
learning approach to geochemical mapping.
Journal of Geochemical Exploration
,
1236
167
, 49–61.
1237
Koch, J., Stisen, S., Refsgaard, J. C., Ernstsen, V., Jakobsen, P. R., & Højberg, A. L.
1238
(2019). Modeling Depth of the Redox Interface at High Resolution at National
1239
Scale Using Random Forest and Residual Gaussian Simulation.
Water Resources
1240
Research
,
55
, 1451–1469.
1241
Kovaˇ
cevi´
c, M., Bajat, B., & Gaji´
c, B. (2010). Soil type classification and estimation
1242
of soil properties using support vector machines.
Geoderma
,
154
, 340–347.
1243
K¨
uhn, I., & Dormann, C. F. (2012). Less than eight (and a half) misconceptions of
1244
spatial analysis.
Journal of Biogeography
,
39
, 995–998.
1245
K¨
uhn, I., Nobis, M. P., & Durka, W. (2009). Combining spatial and phylogenetic
1246
eigenvector filtering in trait analysis.
Global Ecology and Biogeography
,
18
, 745–
1247
758.
1248
44
Lacoste, M., Minasny, B., McBratney, A., Michot, D., Viaud, V., & Walter, C.
1249
(2014). High resolution 3d mapping of soil organic carbon in a heterogeneous
1250
agricultural landscape.
Geoderma
,
213
, 296–311.
1251
Lagacherie, P. (2008). Digital soil mapping: a state of the art. In A. E. Hartemink,
1252
A. McBratney, & M. de Lourdes Mendon¸ca-Santos (Eds.),
Digital Soil Mapping
1253
with Limited Data
(pp. 3–14). Springer, Dordrecht, Netherlands.
1254
Lagacherie, P., & Holmes, S. (1997). Addressing geographical data errors in a classi-
1255
fication tree for soil unit prediction.
International Journal of Geographical Infor-
1256
mation Science
,
11
, 183–198.
1257
Lagacherie, P., & McBratney, A. (2006). Spatial soil information systems and spatial
1258
soil inference systems: perspectives for digital soil mapping.
Developments in Soil
1259
Science
,
31
, 3–22.
1260
Lamichhane, S., Kumar, L., & Wilson, B. (2019). Digital soil mapping algorithms
1261
and covariates for soil organic carbon mapping and their implications: A review.
1262
Geoderma
, (pp. 395–413).
1263
Lapuschkin, S., W¨
aldchen, S., Binder, A., Montavon, G., Samek, W., & M¨
uller, K.-
1264
R. (2019). Unmasking Clever Hans predictors and assessing what machines really
1265
learn.
Nature Communications
,
10
, 1096.
1266
Le Rest, K., Pinaud, D., Monestiez, P., Chadoeuf, J., & Bretagnolle, V. (2014). Spa-
1267
tial leave-one-out cross-validation for variable selection in the presence of spatial
1268
autocorrelation.
Global Ecology and Biogeography
,
23
, 811–820.
1269
Ließ, M., Glaser, B., & Huwe, B. (2012). Uncertainty in the spatial prediction of
1270
soil texture: comparison of regression tree and random forest models.
Geoderma
,
1271
170
, 70–79.
1272
Ließ, M., Schmidt, J., & Glaser, B. (2016). Improving the spatial prediction of soil
1273
organic carbon stocks in a complex tropical mountain landscape by methodological
1274
specifications in machine learning approaches.
PLOS ONE
,
11
, e0153673.
1275
Lipton, Z. C. (2018). The mythos of model interpretability.
Queue
,
16
, 31–57.
1276
Liu, F., Zhang, G.-L., Song, X., Li, D., Zhao, Y., Yang, J., Wu, H., & Yang, F.
1277
(2019). High-resolution and three-dimensional mapping of soil texture of China.
1278
Geoderma
,
361
, 114061.
1279
45
Lorenzetti, R., Barbetti, R., Fantappi`
e, M., L’Abate, G., & Costantini, E. A. (2015).
1280
Comparing data mining and deterministic pedology to assess the frequency of wrb
1281
reference soil groups in the legend of small scale maps.
Geoderma
,
237
, 237–245.
1282
Ma, Y., Minasny, B., Malone, B. P., & Mcbratney, A. B. (2019). Pedology and digital
1283
soil mapping (DSM).
European Journal of Soil Science
,
70
, 216–235.
1284
MacKay, D. J. C. (1992). Information-based objective functions for active data
1285
selection.
Neural Computation
,
4
, 590–604.
1286
Mahmoudabadi, E., Karimi, A., Haghnia, G. H., & Sepehr, A. (2017). Digital soil
1287
mapping using remote sensing indices, terrain attributes, and vegetation features in
1288
the rangelands of northeastern Iran.
Environmental Monitoring and Assessment
,
1289
189
, 500.
1290
Malone, B. P., McBratney, A. B., Minasny, B., & Laslett, G. M. (2009). Mapping
1291
continuous depth functions of soil carbon storage and available water capacity.
1292
Geoderma
,
154
, 138–152.
1293
Mansuy, N., Thiffault, E., Par´
e, D., Bernier, P., Guindon, L., Villemaire, P., Poirier,
1294
V., & Beaudoin, A. (2014). Digital mapping of soil properties in Canadian managed
1295
forests at 250 m of resolution using the k-nearest neighbor method.
Geoderma
,
1296
235
, 59–73.
1297
Marchant, B. P., Viscarra Rossel, R. A., & Webster, R. (2013). Fluctuations in
1298
method-of-moments variograms caused by clustered sampling and their elimination
1299
by declustering and residual maximum likelihood estimation.
European Journal of
1300
Soil Science
,
64
, 401–409.
1301
Massawe, B. H. J., Subburayalu, S. K., Kaaya, A. K., Winowiecki, L., & Slater, B. K.
1302
(2018). Mapping numerically classified soil taxa in Kilombero valley, Tanzania
1303
using machine learning.
Geoderma
,
311
, 143–148.
1304
Mayr, T., Rivas-Casado, M., Bellamy, P., Palmer, R., Zawadzka, J., & Corstanje, R.
1305
(2010). Two methods for using legacy data in digital soil mapping. In M. A. H. A.
1306
K.-B. S. Boettinger J.L., Howell D.W. (Ed.),
Digital Soil Mapping
(pp. 191–202).
1307
Springer, Dordrecht, NL.
1308
McBratney, A. B., Santos, M. M., & Minasny, B. (2003). On digital soil mapping.
1309
Geoderma
,
117
, 3–52.
1310
46
McNicol, G., Bulmer, C., D’Amore, D., Sanborn, P., Saunders, S., Giesbrecht, I.,
1311
Arriola, S. G., Bidlack, A., Butman, D., & Buma, B. (2019). Large, climate-
1312
sensitive soil carbon stocks mapped with pedology-informed machine learning in
1313
the North Pacific coastal temperate rainforest.
Environmental Research Letters
,
1314
14
, 014004.
1315
Meyer, H., Reudenbach, C., Hengl, T., Katurji, M., & Nauss, T. (2018). Improving
1316
performance of spatio-temporal machine learning models using forward feature
1317
selection and target-oriented validation.
Environmental Modelling & Software
,
1318
101
, 1–9.
1319
Meyer, H., Reudenbach, C., W¨
ollauer, S., & Nauss, T. (2019). Importance of spatial
1320
predictor variable selection in machine learning applications–moving from data
1321
reproduction to spatial prediction.
Ecological Modelling
,
411
, 108815.
1322
Micheletti, N., Foresti, L., Robert, S., Leuenberger, M., Pedrazzini, A., Jaboyed-
1323
off, M., & Kanevski, M. (2014). Machine learning feature selection methods for
1324
landslide susceptibility mapping.
Mathematical Geosciences
,
46
, 33–57.
1325
Miller, B. A., Koszinski, S., Wehrhan, M., & Sommer, M. (2015a). Comparison of
1326
spatial association approaches for landscape mapping of soil organic carbon stocks.
1327
Soil
,
1
, 217––233.
1328
Miller, B. A., Koszinski, S., Wehrhan, M., & Sommer, M. (2015b). Impact of multi-
1329
scale predictor selection for modeling soil properties.
Geoderma
,
239
, 97–106.
1330
Miller, T. (2019). Explanation in artificial intelligence: Insights from the social
1331
sciences.
Artificial Intelligence
,
267
, 1–38.
1332
Mira, M., Weiss, M., Baret, F., Courault, D., Hagolle, O., Gallego-Elvira, B., &
1333
Olioso, A. (2015). The MODIS (collection V006) BRDF/albedo product MCD43D:
1334
Temporal course evaluated over agricultural landscape.
Remote Sensing of Envi-
1335
ronment
,
170
, 216–228.
1336
Mirza, M., & Osindero, S. (2014). Conditional generative adversarial nets.
1337
1338
Molnar, C. (2019).
Interpretable machine learning
. Lulu, Morrisville, USA.
1339
Montavon, G., Samek, W., & M¨
uller, K.-R. (2018). Methods for interpreting and
1340
understanding deep neural networks.
Digital Signal Processing
,
73
, 1–15.
1341
47
Moran, C. J., & Bui, E. N. (2002). Spatial data mining for enhanced soil map
1342
modelling.
International Journal of Geographical Information Science
,
16
, 533–
1343
549.
1344
Mosleh, Z., Salehi, M. H., Jafari, A., Borujeni, I. E., & Mehnatkesh, A. (2016). The
1345
effectiveness of digital soil mapping to predict soil properties over low-relief areas.
1346
Environmental Monitoring and Assessment
,
188
, 195.
1347
Mulder, V. L., Lacoste, M., Richer-de Forges, A. C., Martin, M. P., & Arrouays, D.
1348
(2016). National versus global modelling the 3D distribution of soil organic carbon
1349
in mainland France.
Geoderma
,
263
, 16–34.
1350
Nauman, T. W., & Duniway, M. C. (2019). Relative prediction intervals reveal larger
1351
uncertainty in 3D approaches to predictive digital soil mapping of soil properties
1352
with legacy data.
Geoderma
,
347
, 170–184.
1353
Nussbaum, M., Spiess, K., Baltensweiler, A., Grob, U., Keller, A., Greiner, L.,
1354
Schaepman, M. E., & Papritz, A. (2018). Evaluation of digital soil mapping ap-
1355
proaches with large sets of environmental covariates.
Soil
,
4
, 1–22.
1356
Odgers, N. P., Sun, W., McBratney, A. B., Minasny, B., & Clifford, D. (2014).
1357
Disaggregating and harmonising soil map units through resampled classification
1358
trees.
Geoderma
,
214
, 91–100.
1359
Oldeman, L. R., & Van Engelen, V. W. P. (1993). A world soils and terrain digital
1360
database (SOTER)—An improved assessment of land resources.
Geoderma
,
60
,
1361
309–325.
1362
Oliver, M. A. (1987). Geostatistics and its application to soil science.
Soil Use and
1363
Management
,
3
, 8–20.
1364
Ottoy, S., De Vos, B., Sindayihebura, A., Hermy, M., & Van Orshoven, J. (2017).
1365
Assessing soil organic carbon stocks under current and potential forest cover using
1366
digital soil mapping and spatial generalisation.
Ecological Indicators
,
77
, 139–150.
1367
Padarian, J., Minasny, B., & McBratney, A. B. (2019). Using deep learning for
1368
digital soil mapping.
Soil
,
5
, 79–89.
1369
Pahlavan-Rad, M. R., & Akbarimoghaddam, A. (2018). Spatial variability of soil
1370
texture fractions and pH in a flood plain (case study from eastern Iran).
CATENA
,
1371
160
, 275–281.
1372
48
Peres-Neto, P. R., Legendre, P., Dray, S., & Borcard, D. (2006). Variation partition-
1373
ing of species data matrices: estimation and comparison of fractions.
Ecology
,
87
,
1374
2614–2625.
1375
Poggio, L., Lassauce, A., & Gimona, A. (2019). Modelling the extent of northern
1376
peat soil and its uncertainty with sentinel: Scotland as example of highly cloudy
1377
region.
Geoderma
,
346
, 63–74.
1378
Pohjankukka, J., Pahikkala, T., Nevalainen, P., & Heikkonen, J. (2017). Estimating
1379
the prediction performance of spatial models via spatial k-fold cross validation.
1380
International Journal of Geographical Information Science
,
31
, 2001–2019.
1381
Pouladi, N., Møller, A. B., Tabatabai, S., & Greve, M. H. (2019). Mapping soil
1382
organic matter contents at field level with cubist, random forest and kriging.
Geo-
1383
derma
,
342
, 85–92.
1384
Pozdnoukhov, A., & Kanevski, M. (2006). Monitoring network optimisation for
1385
spatial data classification using support vector machines.
International Journal of
1386
Environment and Pollution
,
28
, 465–484.
1387
Rahman, R., Otridge, J., & Pal, R. (2017). IntegratedMRF: random forest-based
1388
framework for integrating prediction from different data types.
Bioinformatics
,
1389
33
, 1407–1410.
1390
Ramcharan, A., Hengl, T., Nauman, T., Brungard, C., Waltman, S., Wills, S., &
1391
Thompson, J. (2018). Soil property and class maps of the conterminous United
1392
States at 100-meter spatial resolution.
Soil Science Society of America Journal
,
1393
82
, 186–201.
1394
Reichstein, M., Camps-Valls, G., Stevens, B., Jung, M., Denzler, J., Carvalhais,
1395
N. et al. (2019). Deep learning and process understanding for data-driven earth
1396
system science.
Nature
,
566
, 195–204.
1397
Ribeiro, M. T., Singh, S., & Guestrin, C. (2016). Why should I trust you?: Explain-
1398
ing the predictions of any classifier. In
Proceedings of the 22nd ACM SIGKDD
1399
International Conference on Knowledge Discovery and Data Mining
(pp. 1135–
1400
1144). ACM.
1401
Roberts, D. R., Bahn, V., Ciuti, S., Boyce, M. S., Elith, J., Guillera-Arroita, G.,
1402
Hauenstein, S., Lahoz-Monfort, J. J., Schr¨
oder, B., Thuiller, W. et al. (2017).
1403
Cross-validation strategies for data with temporal, spatial, hierarchical, or phylo-
1404
genetic structure.
Ecography
,
40
, 913–929.
1405
49
Roscher, R., Bohn, B., Duarte, M. F., & Garcke, J. (2019). Explainable machine
1406
learning for scientific insights and discoveries.
1407
Rudiyanto, B., Minasny, Setiawan, B. I., Saptomo, S. K., McBratney, A. B. et al.
1408
(2018). Open digital mapping as a cost-effective method for mapping peat thickness
1409
and assessing the carbon stock of tropical peatlands.
Geoderma
,
313
, 25–40.
1410
Ruß, G., & Brenning, A. (2010). Data mining in precision agriculture: management
1411
of spatial information. In
International Conference on Information Processing and
1412
Management of Uncertainty in Knowledge-Based Systems
(pp. 350–359). Springer.
1413
Schratz, P., Muenchow, J., Iturritxa, E., Richter, J., & Brenning, A. (2019). Hyper-
1414
parameter tuning and performance assessment of statistical and machine-learning
1415
algorithms using spatial data.
Ecological Modelling
,
406
, 109–120.
1416
Scull, P., Franklin, J., & Chadwick, O. A. (2005). The application of classification
1417
tree analysis to soil type prediction in a desert landscape.
Ecological Modelling
,
1418
181
, 1–15.
1419
Segal, M., & Xiao, Y. (2011). Multivariate random forests.
Wiley Interdisciplinary
1420
Reviews: Data Mining and Knowledge Discovery
,
1
, 80–87.
1421
Sergeev, A., Buevich, A., Baglaeva, E., & Shichkin, A. (2019). Combining spatial
1422
autocorrelation with machine learning increases prediction accuracy of soil heavy
1423
metals.
CATENA
,
174
, 425–435.
1424
Sharififar, A., Sarmadian, F., Malone, B. P., & Minasny, B. (2019). Addressing
1425
the issue of digital mapping of soil classes with imbalanced class observations.
1426
Geoderma
,
350
, 84–92.
1427
Shi, J., Yang, L., Zhu, A., Qin, C., Liang, P., Zeng, C., Pei, T. et al. (2018). Machine-
1428
learning variables at different scales vs. knowledge-based variables for mapping
1429
multiple soil properties.
Soil Science Society of America Journal
,
82
, 645–656.
1430
Siewert, M. B. (2018). High-resolution digital mapping of soil organic carbon in
1431
permafrost terrain using machine learning: a case study in a sub-Arctic peatland
1432
environment.
Biogeosciences
,
15
, 1663–1682.
1433
da Silva Chagas, C., de Carvalho Junior, W., Bhering, S. B., & Calderano Filho, B.
1434
(2016). Spatial prediction of soil surface texture in a semiarid region using random
1435
forest and multiple linear regressions.
CATENA
,
139
, 232–240.
1436
50
Sinha, P., Gaughan, A. E., Stevens, F. R., Nieves, J. J., Sorichetta, A., & Tatem, A. J.
1437
(2019). Assessing the spatial sensitivity of a random forest model: Application in
1438
gridded population modeling.
Computers, Environment and Urban Systems
,
75
,
1439
132–145.
1440
Somarathna, P. D. S. N., Malone, B. P., & Minasny, B. (2016). Mapping soil organic
1441
carbon content over New South Wales, Australia using local regression kriging.
1442
Geoderma regional
,
7
, 38–48.
1443
Somarathna, P. D. S. N., Minasny, B., & Malone, B. P. (2017). More data or a better
1444
model? figuring out what matters most for the spatial prediction of soil carbon.
1445
Soil Science Society of America Journal
,
81
, 1413–1426.
1446
Song, X.-D., Wu, H.-Y., Ju, B., Liu, F., Yang, F., Li, D.-C., Zhao, Y.-G., Yang, J.-
1447
L., & Zhang, G.-L. (2020). Pedoclimatic zone-based three-dimensional soil organic
1448
carbon mapping in China.
Geoderma
,
363
, 114145.
1449
Song, Y.-Q., Zhao, X., Su, H.-Y., Li, B., Hu, Y.-M., & Cui, X.-S. (2018). Predicting
1450
spatial variations in soil nutrients with hyperspectral remote sensing at regional
1451
scale.
Sensors
,
18
, 3086.
1452
Stojanova, D., Ceci, M., Appice, A., Malerba, D., & Dˇ
zeroski, S. (2013). Dealing
1453
with spatial autocorrelation when learning predictive clustering trees.
Ecological
1454
Informatics
,
13
, 22–39.
1455
Stumpf, F., Schmidt, K., Behrens, T., Sch¨
onbrodt-Stitt, S., Buzzo, G., Dumperth,
1456
C., Wadoux, A., Xiang, W., & Scholten, T. (2016). Incorporating limited field
1457
operability and legacy soil samples in a hypercube sampling design for digital soil
1458
mapping.
Journal of Plant Nutrition and Soil Science
,
179
, 499–509.
1459
Subburayalu, S. K., & Slater, B. K. (2013). Soil series mapping by knowledge dis-
1460
covery from an Ohio county soil map.
Soil Science Society of America Journal
,
1461
77
, 1254–1268.
1462
Szatm´
ari, G., & P´
asztor, L. (2019). Comparison of various uncertainty modelling
1463
approaches based on geostatistics and machine learning algorithms.
Geoderma
,
1464
337
, 1329–1340.
1465
Szatm´
ari, G., Pirk´
o, B., Ko´
os, S., Laborczi, A., Bakacsi, Z., Szab´
o, J., & P´
asztor,
1466
L. (2019). Spatio-temporal assessment of topsoil organic carbon stock change in
1467
Hungary.
Soil and Tillage Research
,
195
, 104410.
1468
51
Taalab, K., Corstanje, R., Mayr, T., Whelan, M., & Creamer, R. (2015).
The
1469
application of expert knowledge in bayesian networks to predict soil bulk density
1470
at the landscape scale.
European Journal of Soil Science
,
66
, 930–941.
1471
Taghizadeh-Mehrjardi, R., Minasny, B., Sarmadian, F., & Malone, B. P. (2014).
1472
Digital mapping of soil salinity in Ardakan region, central Iran.
Geoderma
,
213
,
1473
15–28.
1474
Taghizadeh-Mehrjardi, R., Minasny, B., Toomanian, N., Zeraatpisheh, M., Amirian-
1475
Chakan, A., & Triantafilis, J. (2019a). Digital mapping of soil classes using en-
1476
semble of models in Isfahan region, Iran.
Soil Systems
,
3
, 37.
1477
Taghizadeh-Mehrjardi, R., Nabiollahi, K., & Kerry, R. (2016). Digital mapping of
1478
soil organic carbon at multiple depths using different data mining techniques in
1479
Baneh region, Iran.
Geoderma
,
266
, 98–110.
1480
Taghizadeh-Mehrjardi, R., Schmidt, K., Eftekhari, K., Behrens, T., Jamshidi, M.,
1481
Davatgaar, N., Toomanian, N., & Scholten, T. (2019b). Synthetic resampling
1482
strategies and machine learning for digital soil mapping in Iran.
European Journal
1483
of Soil Science
, .
1484
Taghizadeh-mehrjardi, R., Toomanian, N., Khavaninzadeh, A. R., Jafari, A., & Tri-
1485
antafilis, J. (2016). Predicting and mapping of soil particle-size fractions with
1486
adaptive neuro-fuzzy inference and ant colony optimization in central Iran.
Euro-
1487
pean Journal of Soil Science
,
67
, 707–725.
1488
Tajik, S., Ayoubi, S., Shirani, H., & Zeraatpisheh, M. (2019). Digital mapping of
1489
soil invertebrates using environmental attributes in a deciduous forest ecosystem.
1490
Geoderma
,
353
, 252–263.
1491
Tuia, D., Pozdnoukhov, A., Foresti, L., & Kanevski, M. (2013). Active learning for
1492
monitoring network optimization. In J. Mateu, & W. G. M¨
uller (Eds.),
Spatio-
1493
Temporal Design: Advances in Efficient Data Acquisition
(pp. 285–318). Wiley
1494
Online Library, Chichester, UK.
1495
Tziachris, P., Aschonitis, V., Chatzistathis, T., & Papadopoulou, M. (2019). As-
1496
sessment of spatial hybrid methods for predicting soil organic matter using DEM
1497
derivatives and soil parameters.
CATENA
,
174
, 206–216.
1498
Uria, B., Cˆ
ot´
e, M.-A., Gregor, K., Murray, I., & Larochelle, H. (2016). Neural
1499
autoregressive distribution estimation.
The Journal of Machine Learning Research
,
1500
17
, 7184–7220.
1501
52
Vaysse, K., & Lagacherie, P. (2015). Evaluating digital soil mapping approaches for
1502
mapping GlobalSoilMap soil properties from legacy data in Languedoc-Roussillon
1503
(France).
Geoderma Regional
,
4
, 20–30.
1504
Vaysse, K., & Lagacherie, P. (2017). Using quantile regression forest to estimate
1505
uncertainty of digital soil mapping products.
Geoderma
,
291
, 55–64.
1506
Vermeulen, D., & Van Niekerk, A. (2017). Machine learning performance for pre-
1507
dicting soil salinity using different combinations of geomorphometric covariates.
1508
Geoderma
,
299
, 1–12.
1509
Vincent, S., Lemercier, B., Berthier, L., & Walter, C. (2018). Spatial disaggregation
1510
of complex soil map units at the regional scale based on soil-landscape relation-
1511
ships.
Geoderma
,
311
, 130–142.
1512
Viscarra-Rossel, R. A., & Chen, C. (2011). Digitally mapping the information content
1513
of visible–near infrared spectra of surficial Australian soils.
Remote Sensing of
1514
Environment
,
115
, 1443–1455.
1515
Viscarra-Rossel, R. A., Chen, C., Grundy, M. J., Searle, R., Clifford, D., & Campbell,
1516
P. H. (2015). The Australian three-dimensional soil grid: Australia’s contribution
1517
to the GlobalSoilMap project.
Soil Research
,
53
, 845–864.
1518
Wadoux, A. M. J.-C. (2019a).
Sampling design optimization for geostatistical mod-
1519
elling and prediction
. Ph.D. thesis Wageningen University & Research.
1520
Wadoux, A. M. J.-C. (2019b). Using deep learning for multivariate mapping of soil
1521
with quantified uncertainty.
Geoderma
,
351
, 59–70.
1522
Wadoux, A. M. J.-C., Brus, D. J., & Heuvelink, G. B. M. (2019a). Sampling design
1523
optimization for soil mapping with random forest.
Geoderma
,
355
, 113913.
1524
Wadoux, A. M. J.-C., Padarian, J., & Minasny, B. (2019b).
Multi-source data
1525
integration for soil mapping using deep learning.
Soil
,
5
, 107–119.
1526
Wadoux, A. M. J.-C., Samuel-Rosa, A., Poggio, L., & Mulder, V. L. (2019c). A note
1527
on knowledge discovery and machine learning in digital soil mapping.
European
1528
Journal of Soil Science
, .
1529
Wang, B., Waters, C., Orgill, S., Gray, J., Cowie, A., Clark, A., & Li Liu, D.
1530
(2018). High resolution mapping of soil organic carbon stocks using remote sensing
1531
variables in the semi-arid rangelands of eastern Australia.
Science of The Total
1532
Environment
,
630
, 367–378.
1533
53
Wang, S., Zhuang, Q., Wang, Q., Jin, X., & Han, C. (2017). Mapping stocks of soil
1534
organic carbon and soil total nitrogen in Liaoning Province of China.
Geoderma
,
1535
305
, 250–263.
1536
Were, K., Bui, D. T., Dick, Ø. B., & Singh, B. R. (2015). A comparative assessment
1537
of support vector regression, artificial neural networks, and random forests for pre-
1538
dicting and mapping soil organic carbon stocks across an Afromontane landscape.
1539
Ecological Indicators
,
52
, 394–403.
1540
Wiesmeier, M., Barthold, F., Blank, B., & K¨
ogel-Knabner, I. (2011). Digital mapping
1541
of soil organic matter stocks using random forest modeling in a semi-arid steppe
1542
ecosystem.
Plant and Soil
,
340
, 7–24.
1543
Wu, J., Teng, Y., Chen, H., & Li, J. (2016). Machine-learning models for on-site
1544
estimation of background concentrations of arsenic in soils using soil formation
1545
factors.
Journal of Soils and Sediments
,
16
, 1787–1797.
1546
Xu, S., An, X., Qiao, X., Zhu, L., & Li, L. (2013). Multi-output least-squares support
1547
vector regression machines.
Pattern Recognition Letters
,
34
, 1078–1084.
1548
Yamazaki, D., Ikeshima, D., Tawatari, R., Yamaguchi, T., O’Loughlin, F., Neal,
1549
J. C., Sampson, C. C., Kanae, S., & Bates, P. D. (2017). A high-accuracy map of
1550
global terrain elevations.
Geophysical Research Letters
,
44
, 5844–5853.
1551
Yang, R.-M., Zhang, G.-L., Liu, F., Lu, Y.-Y., Yang, F., Yang, F., Yang, M., Zhao,
1552
Y.-G., & Li, D.-C. (2016). Comparison of boosted regression tree and random forest
1553
models for mapping topsoil organic carbon concentration in an alpine ecosystem.
1554
Ecological Indicators
,
60
, 870–878.
1555
Zeraatpisheh, M., Ayoubi, S., Jafari, A., & Finke, P. (2017). Comparing the efficiency
1556
of digital and conventional soil mapping to predict soil types in a semi-arid region
1557
in Iran.
Geomorphology
,
285
, 186–204.
1558
Zeraatpisheh, M., Ayoubi, S., Jafari, A., Tajik, S., & Finke, P. (2019). Digital
1559
mapping of soil properties using multiple machine learning in a semi-arid region,
1560
central Iran.
Geoderma
,
338
, 445–452.
1561
Zhu, M., Feng, Q., Zhang, M., Liu, W., Deo, R. C., Zhang, C., & Yang, L. (2019).
1562
Soil organic carbon in semiarid alpine regions: the spatial distribution, stock esti-
1563
mation, and environmental controls.
Journal of Soils and Sediments
,
19
, 1–15.
1564
54