SOIL, 7, 217–240, 2021
https://doi.org/10.5194/soil-7-217-2021
© Author(s) 2021. This work is distributed under
the Creative Commons Attribution 4.0 License.
SOIL
SoilGrids 2.0: producing soil information for the globe
with quantified spatial uncertainty
Laura Poggio, Luis M. de Sousa, Niels H. Batjes, Gerard B. M. Heuvelink, Bas Kempen, Eloi Ribeiro,
and David Rossiter
ISRIC – World Soil Information, Wageningen, the Netherlands
Correspondence:
Laura Poggio (laura.poggio@wur.nl)
Received: 14 October 2020 – Discussion started: 9 November 2020
Revised: 9 April 2021 – Accepted: 18 April 2021 – Published: 14 June 2021
Abstract.
SoilGrids produces maps of soil properties for the entire globe at medium spatial resolution (250 m
cell size) using state-of-the-art machine learning methods to generate the necessary models. It takes as inputs soil
observations from about 240 000 locations worldwide and over 400 global environmental covariates describing
vegetation, terrain morphology, climate, geology and hydrology. The aim of this work was the production of
global maps of soil properties, with cross-validation, hyper-parameter selection and quantification of spatially
explicit uncertainty, as implemented in the SoilGrids version 2.0 product incorporating state-of-the-art practices
and adapting them for global digital soil mapping with legacy data. The paper presents the evaluation of the
global predictions produced for soil organic carbon content, total nitrogen, coarse fragments, pH (water), cation
exchange capacity, bulk density and texture fractions at six standard depths (up to 200 cm). The quantitative
evaluation showed metrics in line with previous global, continental and large-region studies. The qualitative
evaluation showed that coarse-scale patterns are well reproduced. The spatial uncertainty at global scale high-
lighted the need for more soil observations, especially in high-latitude regions.
1
Introduction
Healthy soils provide important ecosystem services at the lo-
cal, landscape and global level and are important for the func-
tioning of terrestrial ecosystems (Banwart et al., 2014; FAO
and ITPS, 2015; UNEP, 2012). Information on world soil re-
sources, based on the currently “best available” (shared) soil
profile data, at a scale level commensurate with user needs, is
required to address a range of pressing global issues. These
include avoiding and reducing soil erosion through land re-
habilitation and development (Borrelli et al., 2017; WOCAT,
2007), mitigating and adapting to climate change (Batjes,
2019; Harden et al., 2017; Sanderman et al., 2017; Yigini
and Panagos, 2016; Smith et al., 2019) and ensuring water
security (Rockstroem et al., 2012), food production and food
security (FAO et al., 2018; Soussana et al., 2017; Spring-
mann et al., 2018), as well as preserving biodiversity (Barnes,
2015; IPBES, 2019; van der Esch et al., 2017) and human
livelihood (Bouma, 2015).
The best available soil data are required to support the
Land Degradation Neutrality (LDN) (Cowie et al., 2018)
initiative, achieve several of the Sustainable Development
Goals and provide input for, for example, Earth system mod-
elling by the IPCC (Dai et al., 2019; Luo et al., 2016;
Todd-Brown et al., 2013) and crop modelling (Han et al.,
2019; van Bussel et al., 2015; van Ittersum et al., 2013),
among many other applications. Such information can in
turn help inform international conventions such as the United
Nations Framework Convention on Climate Change (UN-
FCCC), the United Nation Convention to Combat Desertifi-
cation (UNCCD) and the United Nations Convention on Bi-
ological Diversity (UNCBD).
Until the last decade, most global scale assessments re-
quiring soil data used the Digital Soil Map of the World
(DSMW) FAO (1995), an updated version of the original
printed 1
:
5
×
10
6
scale Soil Map of the World (SMW)
(FAO-Unesco, 1971–1981). The soil geographic data from
the DSMW provided the basis for generating a range of de-
Published by Copernicus Publications on behalf of the European Geosciences Union.
218
L. Poggio et al.: Soil information for the globe
rived soil property databases that drew on a larger selection
of soil profile data held in the WISE database (Batjes, 2012)
and more sophisticated (taxotransfer) procedures for deriving
various soil properties (Batjes et al., 2007). Subsequently, in
a joint effort coordinated by the Food and Agriculture Or-
ganization of the United Nations (FAO), the best available
(newer) soil information collated for central and southern
Africa, China, Europe, northern Eurasia and Latin America
was combined into a new product known as the Harmonised
World Soil Database (HWSD) (FAO et al., 2012).
Until recently, the HWSD was the only digital map an-
nex database available for global analyses. However, it has
several limitations (GSP and FAO, 2016; Hengl et al., 2014;
Ivushkin et al., 2019; Omuto et al., 2012). Some of these re-
late to the partly outdated soil geographic data, as well as the
use of a two-layer model (0–30 and 30–100 cm) for deriv-
ing soil properties. Others concern the derived attribute data
themselves, in particular their unquantified uncertainty, and
the use of three different versions of the FAO legend (i.e.
FAO74, FAO85 and FAO90). These issues have been ad-
dressed to varying degrees in various new global soil datasets
(Batjes, 2016; Shangguan et al., 2014; Stoorvogel et al.,
2017) that still largely draw on a traditional soil mapping ap-
proach (Dai et al., 2019).
In the last decade, digital soil mapping (DSM) has be-
come a widely used approach to obtain maps of soil infor-
mation (Minasny and McBratney, 2016). DSM consists pri-
marily in building a quantitative numerical model between
soil observations and environmental information acting as
proxies for the soil forming factors (McBratney et al., 2003;
Minasny and McBratney, 2016). DSM can also integrate di-
rect information as proxies for soil properties, for example
proximal sensing measurements. The number of studies us-
ing DSM to produce maps of soil properties is ever grow-
ing. Numerous modelling approaches are considered, from
linear models to geostatistics, machine learning and artifi-
cial intelligence (e.g. deep learning). Keskin and Grunwald
(2018) provide a recent review of methods and applications
in the field of DSM. DSM techniques have been applied at
various spatial resolutions (e.g. 30 to 1000 m) to support pre-
cision farming (e.g. Piikki et al., 2017) as well as applications
at landscape (e.g. Ellili et al., 2019; Kempen et al., 2015),
country (e.g. Mora-Vallejo et al., 2008; Nijbroek et al., 2018;
Vitharana et al., 2019; Poggio and Gimona, 2017b; Kempen
et al., 2019), regional (e.g. Dorji et al., 2014; Moulatlet et al.,
2017), continental (e.g. Grunwald et al., 2011; Guevara et al.,
2018; Hengl et al., 2017a) and global levels (e.g. Hengl et al.,
2014, 2017b; GSP and ITPS, 2018; Stockmann et al., 2015).
The aim of this paper is to present the development of new
soil property maps for the world at 250 m grid resolution with
a process incorporating state-of-the-art practices and adapt-
ing them to the challenges of global digital soil mapping with
legacy data. It builds on previous global soil properties maps
(SoilGrids250m) (Hengl et al., 2017b), integrating up-to-date
machine learning methods, the increased availability of stan-
dardised soil profile data for the world (Batjes et al., 2020)
and environmental covariates (Nussbaum et al., 2018; Pog-
gio et al., 2013; Reuter and Hengl, 2012). In particular, this
paper addresses the following elements at global scale:
1. incorporation of soil profile data derived from IS-
RIC’s World Soil Information Service (WoSIS), with
expanded number and spatial distribution of observa-
tions (Batjes et al., 2020);
2. a reproducible covariate selection procedure, relying on
recursive feature elimination (Guyon et al., 2002);
3. improved cross-validation procedure, based on spatial
stratification; and
4. quantification of prediction uncertainty using quantile
regression forests (Meinshausen, 2006).
2
Materials and methods
This study uses quantile regression forests (Meinshausen,
2006), a method with a limited number of parameters to be
tuned and that has proven to be an effective compromise be-
tween accuracy and feasibility for large datasets. Selected
primary soil properties as defined and described in the Glob-
alSoilMap specifications (Arrouays et al., 2014) were mod-
elled. The following sections describe each step of the work-
flow (Fig. 1) in detail. These include the following:
1. input soil data preparation
2. covariates’ selection
3. model tuning and cross-validation
4. final model fitting for prediction
5. predictions with uncertainty estimation.
2.1
Soil observation data
Soil property data for this study were derived from the IS-
RIC World Soil Information Service (WoSIS), which pro-
vides consistent, standardised soil profile data for the world
(Batjes et al., 2020). All soil data shared with ISRIC to sup-
port global mapping activities are first stored in the ISRIC
Data Repository, together with their metadata (including the
name of the data owner and licence defining access rights).
Subsequently, the source data are imported “as is” into Post-
greSQL, after which they are ingested into the WoSIS data
model itself. Following data quality assessment and con-
trol (including consistency checks on latitude–longitude and
depth of horizon/layer; flagging of duplicate profiles; and
providing measures for geographic and attribute accuracy,
as well as time stamps), the descriptions for the soil ana-
lytical methods and the units of measurement are standard-
ised using consistent procedures, with additional checks for
SOIL, 7, 217–240, 2021
https://doi.org/10.5194/soil-7-217-2021
L. Poggio et al.: Soil information for the globe
219
Figure 1.
Workflow of the methodological approach.
Table 1.
Soil properties description and units.
Soil property
Acronym Units
Mapped units
Description
Bulk density
BDOD
kg/dm
3
cg/cm
3
Bulk density of the fine earth fraction oven dry
Cation exchange
capacity
CEC
cmol(c)/kg
mmol(c)/kg
Capacity of the fine earth fraction to hold exchangeable
cations
Coarse fragments
CFVO
cm
3
/100 cm
3
(volume %)
cm
3
/dm
3
Volumetric content of fragments larger than 2 mm in the
whole soil
Nitrogen
N
g/kg
cg/kg
Sum of total nitrogen (ammonia, organic and re-
duced nitrogen) as measured by Kjeldahl digestion plus
nitrate–nitrite
pH (water)
pH
–
10
∗
Negative common logarithm of the activity of hydro-
nium ions (H
+
) in water
Organic carbon
concentration
SOC
g/kg
dg/kg
Gravimetric content of organic carbon in the fine earth
fraction of the soil
Soil texture fraction
STF
%
g/kg
Gravimetric contents of sand, silt and clay in the fine
earth fraction of the soil
∗
unitless.
possible erroneous entries for the soil analytical data them-
selves (Ribeiro et al., 2018). Ultimately, upon final consis-
tency checks, the standardised data are made available via
the ISRIC Soil Data Hub (https://data.isric.org, last access:
20 May 2021) in accord with the licence specified by the data
providers. As a result, not all data standardised in WoSIS are
freely available to the international community. Hence, this
study considers two “sources” of point data.
The first is the latest publicly available snapshot of WoSIS
(Batjes et al., 2020). It contains, among others, data for
chemical (organic carbon, total nitrogen, soil pH, cation ex-
change capacity) and physical properties (soil texture (sand,
silt and clay), coarse fragments). The snapshot comprises
196 498 georeferenced profiles originating from 173 coun-
tries, representing over 832 000 soil layers (or horizons), in
total over 5.8 million records. Generally, there are more ob-
servations for the superficial than the deeper layers. About
5 % of the profiles were sampled before 1960, 14 % between
1961–1980, 32 % between 1981–2000 and 16 % between
2001–2020; the date of sampling is unknown for 34 % of the
shared profiles (Batjes et al., 2020).
https://doi.org/10.5194/soil-7-217-2021
SOIL, 7, 217–240, 2021
220
L. Poggio et al.: Soil information for the globe
Second, in addition to the freely shareable data, several
soil observation databases in our repository have licences
stipulating that ISRIC may only use them for SoilGrids ap-
plications or visualisations, for example EU-LUCAS (Tóth
et al., 2013) and soil data for the state of Victoria (Australia).
The corresponding source datasets were screened and pro-
cessed using the same procedures as used for the regular
WoSIS workflow (some 42 000 profiles). As a result, some
240 000 profiles in total were used as the data source for the
present 2020 SoilGrids run, comprising more than 920 000
observed soil layers. During data processing some minor cor-
rections were made to the merged input dataset, for example
further depth congruence checks.
2.1.1
Soil properties
For the purposes of SoilGrids, “soil” is up to 2 m thick uncon-
solidated material at the Earth’s epidermis in direct contact
with the atmosphere; thus subaqueous and tidally exposed
soils are not considered here. Neither are materials deeper
than 2 m. This decision has consequences for computations
of total stocks, in particular soil organic carbon.
Table 1 describes the soil properties that are considered in
this version of SoilGrids: organic carbon content, total nitro-
gen content, soil pH (measured in water), cation exchange
capacity, soil texture fractions and proportion of coarse frag-
ments. These properties were modelled for the six standard
depths intervals as defined in the GlobalSoilMap specifica-
tions (Arrouays et al., 2014): 0–5, 5–15, 15–30, 30–60, 60–
100 and 100–200 cm.
“Litter layers” on top of minerals soils were excluded from
further modelling using the following assumptions. Consis-
tency in layer depth (e.g. sequential increase in the upper
and lower depth reported for each layer down the profile)
in WoSIS was checked using automated procedures. In ac-
cord with current internationally accepted conventions, such
depth increments are given as “measured from the surface,
including organic layers and mineral covers” (FAO, 2006;
Schoeneberger et al., 2012). Prior to 1993, however, the start
(zero depth) of the profile was set at the top of the mineral
surface (the solum proper), except when “thick” organic lay-
ers as defined for peat soils (FAO-ISRIC, 1986) were present
at the surface. Then the top of the peat layer was taken as
the soil surface. Organic horizons were recorded as above
and mineral horizons recorded as below, relative to the min-
eral surface (Schoeneberger et al., 2012) (p. 2–6). Insofar as
is possible, “superficial litter” on top of mineral layers was
flagged as an auxiliary (Boolean) variable, also with refer-
ence to the original soil horizon designation when provided,
so it can be filtered out during auxiliary computations of soil
properties.
2.1.2
Transformation of texture data
A transformation was applied to the texture fractions, as fol-
lows. The relative percentage of sand, silt and clay can be
treated as compositional variables, as the sum of the com-
ponents always equals 100 %. Therefore, these components
were transformed using the addictive log ratio (ALR) trans-
formation with the Gauss–Hermite quadrature (Aitchison,
1986). ALR has previously been applied to soil texture data
(Lark and Bishop, 2007; Akpa et al., 2014; Ballabio et al.,
2016; Poggio and Gimona, 2017a), and it has been shown
(Lark and Bishop, 2007) that ALR-transformed variables
preserve information on the spatial correlation and maintain
the compositional integrity of the original components. In
this study, clay was used as the denominator variable. There-
fore the two ALR components that were interpolated can be
defined as
ALR1
=
log
sand
clay
ALR2
=
log
silt
clay
.
(1)
2.1.3
Spatial stratification of observations
Random splitting of profile observations into
n
cross-
validation folds is not suitable in this context, considering
the high spatial variation in observation density as it would
provide biased results (Brus, 2014). For regions like Europe
and North America there are over four profiles per 10 km
2
,
whereas for large countries in Asia, such as Kazakhstan, In-
dia or Mongolia, the number of available profiles is still quite
limited (
<
one profile per 100 km
2
) (see Batjes et al., 2020
for further details).
Therefore, soil observations were spatially stratified in the
geodetic domain to guarantee a balanced spatial distribution
within each cross-validation fold. Spatial strata, in the form
of hexagons, were created with an Icosahedral Snyder Equal-
Area Grid (ISEAG) of aperture 3 and resolution 6, resulting
in 7292 strata (i.e. hexagonal cells), each with an area around
70 000 km
2
. This ISEAG was generated with the
dggridR
package for the R language (Barnes et al., 2016).
The profiles were assigned to 1 of 10 folds, each equally
represented in each stratum, i.e. each cell of the grid pre-
viously described. All observations (layers or horizons) be-
longing to a profile were always in the same fold for both
model calibration and evaluation. The
caret
R package was
used to subdivide the locations in the folds while maintaining
the spatial distribution.
2.2
Environmental covariates
Over 400 geographic layers were available as environmental
covariates for this work. These were chosen for their pre-
sumed relation to the major soil forming factors, including
long-term soil conditions, i.e. the “time” factor. Appendix A
SOIL, 7, 217–240, 2021
https://doi.org/10.5194/soil-7-217-2021
L. Poggio et al.: Soil information for the globe
221
provides a list of the products used as covariates and their
sources. The layers considered can be grouped as follows.
–
Climate: temperature, precipitation, snowfall, cloud
cover, solar radiation, wind speed;
–
ecology: bioclimatic zones and ecophysiographic re-
gions;
–
geology: soil and sedimentary thickness, rock types;
–
land use and cover: from sources such as the Euro-
pean Space Agency (ESA) and U.S. Geological Survey
(USGS);
–
elevation and terrain morphology: including numerous
morphology indexes and landform classes;
–
vegetation indexes: such as the normalised difference
vegetation index (NDVI), enhanced vegetation index
(EVI) and net primary production (NPP);
–
raw bands from Landsat and MODIS products;
–
hydrography: global water table, inundation and glacier
extent, and surface water change.
The average and standard deviation of climatic variables
and vegetation indices over 15 years (2001–2015) were com-
puted from monthly data to capture their seasonal dynamics.
All covariates were projected to a common coordinate ref-
erence system (CRS), i.e. Goode’s homolosine projection for
land masses applied to the WGS84 datum. This projection
was selected since among the equal-area projections sup-
ported by open-source software it is the most effective min-
imising distortions over land (de Sousa et al., 2019). The
projected covariates were imported to GRASS GIS in a nor-
malised raster structure with cells of 250 m by 250 m. Co-
variates, and hence mapped areas, were restricted to land ar-
eas without built-up, water and glacier areas using a mask
created from the ESA Land Cover layer for 2015 (Buchhorn
et al., 2020). Thus properties of urban and subaqueous soils
are not considered.
2.3
Covariates’ selection
Considering the large number of available environmental lay-
ers, a standardised and reproducible procedure to select co-
variates used for modelling was implemented to (i) reduce
redundancy between covariates, (ii) obtain a more parsimo-
nious and computationally efficient model, (iii) decrease the
risk of over-fitting (Gomes et al., 2019) and (iv) avoid a bi-
ased assessment of variable importance (Strobl et al., 2008).
The covariates’ selection procedure consisted of two steps,
de-correlation and recursive feature elimination.
2.3.1
De-correlation analysis
De-correlation analysis was carried out as an initial step to
reduce the redundancy of information from more than 400
environmental layers. Only covariate layers that had a pair-
wise correlation coefficient
<
=
0
.
85 with all other covariates
were included in the subsequent analyses. For each pair of
covariates correlated above this threshold, only the first one
in alphabetical order was selected for inclusion in the mod-
elling phase. This step reduced the number of initial covari-
ates to approximately 150 layers.
2.3.2
Recursive feature elimination
Recursive feature elimination (RFE) (Guyon et al., 2002) is
a methodology that has proven effective to select an optimal
set of covariates for regression trees models (Gomes et al.,
2019; Hounkpatin et al., 2018). In this study, the RFE pro-
cedure implemented in the
caret
package for the R lan-
guage (Kuhn, 2015) was used, as it offers a good compro-
mise between accuracy and computation time. The algorithm
starts by fitting a model using all covariates, assessing its per-
formance and ranking covariate importance. The least impor-
tant covariates are then removed from the pool, and again
the model is fitted and assessed and the least important co-
variates removed. The procedure is repeated down to a pool
between 0 and
n
covariates. This procedure is based on out-
of-bag (OOB) cross-validation and does not test all covari-
ates’ combinations, but it is considered one of the most ro-
bust covariates’ selection approaches for models like random
forests (Nussbaum et al., 2018).
The RFE procedure on the full set of observations and
covariates would prove computationally prohibitive. To im-
prove computational feasibility for large datasets, additional
steps were developed. Four sets of observations were used
for RFE, each obtained using three cross-validation folds (see
Sect. 2.1.3 for further details): set 1 contained folds 1 to 3,
set 2 folds 4 to 6, set 3 folds 7 to 9 and set 4 contained fold 10
and two other random selected folds. In a first step, the RFE
procedure from
caret
was run independently on each set
with default model hyper-parameters for the random forests
algorithm as implemented in the
ranger
package (i.e.
ntree
as 500 and
mtry
as the rounded square root of the number
variables). In each set, the optimal number and combination
of covariates were automatically selected when the model
performances stopped increasing, i.e. when the loss function
reached its minimum. In this study, the loss function was the
OOB root-mean-square error (RMSE).
In the second step, the RFE procedure was applied with
all observations and all covariates selected in at least one of
the four sets used in the previous step. The final covariate set
was the set minimising the loss function.
https://doi.org/10.5194/soil-7-217-2021
SOIL, 7, 217–240, 2021
222
L. Poggio et al.: Soil information for the globe
2.4
Hyper-parameter selection and cross-validation
Figure 2 summarises the approach used for the selection of
the model hyper-parameters and the cross-validation. Further
details are provided in the following sections.
2.4.1
Model tuning and numeric evaluation
Model tuning was performed with a 10-fold cross-validation
procedure applied to multiple combinations of hyper-
parameters.
Different numbers of decision trees (
ntree
parameter) were
combined with different numbers of covariates used in tree
splits (
mtry
parameter). The number of trees was progres-
sively increased with the following values: 100, 150, 200,
250, 500, 750 and 1000. The different
mtry
values were mul-
tiples of the square root of the number of covariates. Four
multipliers were tested, 1 (default in
ranger
), 1.5, 2 and
3. For example, if the RFE procedure identified a set of 50
covariates, the
mtry
values assessed were 7, 11, 14 and 21.
Each of the resulting combinations of
ntree
and
mtry
pa-
rameters was used to train a different model with observa-
tions from nine folds. Predictions were then assessed on the
remaining fold with classical performance measures, i.e. root
mean squared error (RMSE) and model efficiency coefficient
(MEC; Janssen and Heuberger, 1995). MEC is equal to the
fraction of the explained variance based on the 1
:
1 line of
predicted versus observed that is defined as 1 minus the ratio
between residual sum of squares and total sum of squares.
The final hyper-parameter selection was based on an optimi-
sation of model performance and computational constraints,
in this case memory consumption. For example an increase
of the
ntree
parameter above 200 provided a minor increment
in the metrics (usually less than 0.1 %, not reported here)
while requiring considerably more memory and computation
time.
The model evaluation was based on the performance met-
rics of the selected hyper-parameters’ combination. Predic-
tions at the centre of the six standard depth intervals were
compared with observations having the midpoint included
within the considered interval.
2.5
Prediction and uncertainty quantification
2.5.1
Model fit
The final model for each soil property was fitted with all
available observations, the covariates and the hyperparam-
eters selected in the previous steps. Observation depth was
included in the model as a covariate. It was calculated at the
midpoint of the sampled layer or horizon.
Models were obtained with the
ranger
package (Wright
and Ziegler, 2017), with the option
quantreg
to build
quantile random forests (QRF; Meinshausen, 2006). With
this option, the prediction is not a single value, e.g. the av-
erage of predictions from the group of decision trees in the
random forest, but rather a cumulative probability distribu-
tion of the soil property at each location and depth.
For each property (see Table 1) and standard depth from
the GlobalSoilMap specification (0–5, 5–15, 15–30, 30–60,
60–100 and 100–200 cm), four different values were com-
puted to characterise this distribution: median (0
.
50 quantile,
q
0
.
50
), mean, 0
.
05 quantile (
q
0
.
05
) and 0
.
95 quantile (
q
0
.
95
),
i.e. the lower and upper limits of a 90 % prediction inter-
val. This uncertainty interval is as described in the Global-
SoilMap specifications (Arrouays et al., 2014). The predic-
tions were computed for the mid-point of the depth interval
and considered constant for the whole depth interval.
In order to compute the prediction uncertainty for soil tex-
ture, the back-transformation was applied at the level of indi-
vidual tree predictions and the quantiles of the tree prediction
distributions obtained from the resulting values.
2.5.2
Uncertainty
The percentage of cross-validation observations contained
in the 0.9 prediction interval was calculated (prediction in-
terval coverage probability, PICP) (Shrestha and Soloma-
tine, 2006). Ideally the PICP is close to 0.9, indicating that
the uncertainty was correctly assessed. A PICP substantially
greater than 0.9 suggests that the uncertainty was underes-
timated; a substantially smaller PICP indicates that it was
overestimated.
Furthermore, to visualise the uncertainty as a map, the fol-
lowing indicators were calculated:
1. 90th prediction interval (PI90)
PI90
=
q
0
.
95
−
q
0
.
05
;
(2)
2. ratio of the interquartile range over the median (predic-
tion interval ratio, PIR):
PIR
=
q
0
.
95
−
q
0
.
05
q
0
.
50
.
(3)
2.6
Qualitative evaluation of spatial patterns
Expert judgement was used to evaluate the reasonableness
of the maps, by comparing well-known spatial patterns at
global, regional and local scales with SoilGrids predictions
(see Sect. 3.4). Obviously these are not definitive evaluations,
only indicative.
2.7
Software and computational framework
SoilGrids requires an intensive computational workflow,
with numerous steps integrating different software. Soil-
Grids is entirely based on open-source software, in particular
SLURM (Yoo et al., 2003) for job management, GRASS GIS
(GRASS Development Team, 2020) for data and tiles’ man-
agement and R statistical software (R Core Team, 2020) for
model fitting and statistical analysis.
SOIL, 7, 217–240, 2021
https://doi.org/10.5194/soil-7-217-2021
L. Poggio et al.: Soil information for the globe
223
Figure 2.
Detailed workflow for the Hyper-parameter selection and cross-validation.
Predictions were computed in a high-performance com-
puting cluster. A dynamic geographic tiling system was de-
veloped with GRASS GIS to maximise the use of memory
for each job. Technical details on this parallelisation scheme
are given in de Sousa et al. (2020).
The predictions were multiplied by a conversion factor of
10 or 100 to maintain the required precision while using in-
teger type in the file geotiff to reduce space occupied on disk.
Application of the conversion factor resulted in mapped lay-
ers with units differing from those of the input observations
(see Table 1).
The total computation time with the selected covariates
and hyper-parameters differed per property. On average, the
complete computation of the 24 maps (mean and three quan-
tiles for each of the six standard depths) for a single prop-
erty, including (i) RFE, (ii) model training and (iii) predic-
tion, took approximately 1500 CPU hours. The prediction ac-
counted for about two-thirds of the total time.
3
Results and discussion
3.1
Input soil observations
Table 2 breaks down the distribution of the legacy soil ob-
servations for each soil property by depth interval. Table B1,
in Appendix B, shows the number of observations by biocli-
matic region.
Figures 3 and 4 show examples of observation density of
the soil calibration data for two soil properties, pH
water
and
proportion of coarse fragments, that show a large difference
in density.
As indicated, the number of observations for each property
varies greatly with depth and bioclimatic region, with higher
densities observed for North America and Europe. Generally,
there are more observations for agricultural areas. Further,
the available profiles have been collated over several decades,
some 62 % of the data being from 1960–2020; the time of
sampling is unknown for around 34 % of the profiles. As in-
dicated by Batjes et al. (2020), in principle, the age of the
observations should be taken into account during the map-
ping process via covariate layers for time periods commen-
surate with the sampling dates, especially for soil properties
that are readily affected by changes in land use or manage-
ment practices. However, for these so-called “dynamic” soil
properties, such as pH and soil organic matter content, we
consider that the spatial variation will be much greater that
the temporal variation, so that not taking the age of observa-
tions into account will not greatly affect the map. In addition,
it is difficult or impossible to find comparable covariates, in
particular remote-sensing-derived covariates, for each time
period. Space–time relations should be considered in future
assessments (Heuvelink et al., 2020).
This study considers standardised data for some 240 000
profiles, derived from WoSIS. This is over 60 000 more pro-
files than considered in the data compilation underpinning
the preceding SoilGrids runs (Hengl et al., 2017b), thus pro-
viding substantial new information for calibration of the new
global models. However, as indicated, there are still signif-
icant geographic gaps (e.g. arid regions, boreal regions and
“forest” soils). Some of these are related to the physical re-
moteness or inaccessibility of some regions, while others are
related to the fact that many soil datasets still are not or can
https://doi.org/10.5194/soil-7-217-2021
SOIL, 7, 217–240, 2021
224
L. Poggio et al.: Soil information for the globe
not be shared for various reasons as described by Arrouays
et al. (2017).
In the previous version of SoilGrids (Hengl et al., 2017b),
synthetic observations were randomly placed in regions with
few or no observations, e.g. the Sahara and the Arabian
Peninsula. This approach is worth further exploring, includ-
ing information derived from other regional datasets, expert
opinion and transfer learning from similar areas according to
the
Homosoil
concept (Mallavan et al., 2010), which assumes
similarity of soil-forming factors across regions. However,
SoilGrids already implicitly incorporates the Homosoil con-
cept, as long as there are sufficient observations in a given
soil-forming environment anywhere in the world. Therefore,
no synthetic observations (“pseudo-points”) were included in
this version of SoilGrids, also by a lack of confidence about
the accuracy of the synthetic data.
In future studies, it will be relevant to identify beforehand
areas of the world with a low observation density that are not
yet represented by a high density of observations in other
areas with similar soil-forming factors. A set of synthetic
profiles could then be generated to describe these areas, by
consulting soil scientists knowledgeable on the soils and soil
properties of these areas.
3.2
Model tuning and hyper-parameter selection
Model hyper-parameters selected for each property are pre-
sented in Table 3.
The numbers of covariates selected using the two-step ap-
proach for covariates’ selection was fairly small in compari-
son with the full set (Table 3), resulting in more parsimonious
models. Figure 5 shows two examples of the loss function for
RFE for two soil properties with different numbers and distri-
butions of input observations. In both cases, there is a clear
improvement of performances while using 15 to 20 covari-
ates. The curve reaches a minimum of the loss function and
then stays on a plateau with a slight decline after the identi-
fied minimum.
All final models were trained with a maximum of 200 de-
cision trees, a number beyond which performance gains did
not noticeably increase.
The
mtry
parameter mainly depended on the number of co-
variates and was always between 1.5 and 2 times the square
root of the number of covariates, which is the default pro-
vided by common random forest packages such as
ranger
(Wright and Ziegler, 2017). This confirms the need to de-
termine optimum model hyper-parameters, especially when
dealing with large numbers of input data (Nussbaum et al.,
2018) as is the case here.
3.3
Quantitative evaluation
Cross-validation results are summarised in Table 4, present-
ing the root mean squared error (RMSE) and model effi-
ciency coefficient (MEC). The MEC varies from a mini-
mum of 0.31 for coarse fragments to a maximum of 0.74
for BDOD. Clay is less well modelled than the other two
particle-size classes. This may be an effect of the chosen
ALR transformation that had clay as denominator (Lark and
Bishop, 2007). Metrics of the mean were always better than
or equal to those for the median for all properties.
Overall, these metrics are in line with continental or large-
region DSM studies (Keskin and Grunwald, 2018). However,
they are slightly lower than those presented by Hengl et al.
(2017b). The latter difference can be explained by the more
prudent cross-validation approach now taken, with spatially
balanced folds and all observations belonging to the same
profile in the same fold. This prevents the use of data from the
same profile both for calibration and numerical evaluation.
Table 4 shows that the models with a higher number of
retained covariates (Table 3) have better predictive perfor-
mances. However, these models are also the models with the
largest number of observations (Table 2). The considered soil
properties are also different. Therefore, no general conclu-
sion can be drawn from this observation.
Table 5 shows the MEC for mean predictions by depth in-
terval. Performances decreased with depth, in line with many
other DSM studies (Keskin and Grunwald, 2018). This pat-
tern can be explained mainly by weakened relationships be-
tween environmental layers and soil properties of the deeper
layers.
In this study, the vertical dimension of soil variability was
only taken into account by using the depth of the observa-
tion as a covariate. Recent publications (Ma et al., 2021;
Nauman and Duniway, 2019) indicate that such an approach
can be too simplistic or lead to problems with consistency
over the predicted depth sequence. This may be true for lo-
cal datasets, in which the short-range spatial variability is
of a similar magnitude as the vertical variability. Further re-
search is necessary to assess the effects of using depth as a
covariate on global datasets and models. Alternatives such
as 3D smoothers (Poggio and Gimona, 2017b) or geostatis-
tical models exploiting 3D spatial auto-correlation are worth
exploring in further studies.
Table 6 summarises the PICPs, globally and by predicted
depth interval. Most of the values are between 0.88 and 0.92,
indicating that the prediction intervals obtained with QRF are
a realistic representation of the prediction uncertainty, as the
expected value for a 90 % prediction interval is 0.90. Excep-
tions are the models for coarse fragments with higher values
around 0.95, indicating an overestimation of prediction un-
certainty. The texture components have values with a larger
spread, around 0.78 to 0.80 for sand and closer to 0.96 for silt
and clay. These indicate a potential underestimation of pre-
diction intervals for sand and overestimation for silt and clay.
These results may be related with the range of these prop-
erties in the input observations. The transformation method
used to derive the prediction intervals for the texture compo-
nents could also be a contributing factor. Further exploration
of the causes is worthwhile.
SOIL, 7, 217–240, 2021
https://doi.org/10.5194/soil-7-217-2021
L. Poggio et al.: Soil information for the globe
225
Table 2.
Number of observations per standard depth interval for each soil property. See Table 1 for abbreviations and units of the soil
properties considered.
Depth interval
BDOD
CEC
CFVO
N
pH
SOC
STF
0–5 cm
8122
20 576
15 541
27192
44 049
48 616
42 983
5–15 cm
19 817
49 463
66 833
82 856
146 677
148 918
155 302
15–30 cm
17 819
40 673
35 254
39 568
91 326
91 682
98 659
30–60 cm
27 146
63 444
56 755
48 804
141 812
122 338
140 353
60–100 cm
23 130
58 038
50 912
36 946
131 172
102 687
12 7073
100–200 cm
23 396
66 236
49 995
28 135
129 373
92 327
116 847
Figure 3.
Number of observations per grid cell (70 000 km
2
) for soil pH
water
.
Table 3.
Hyper-parameters for each considered soil property. See
Table 1 for abbreviations and units of the soil properties considered.
Number of
Number of
mtry
covariates
trees
BDOD
40
200
12
CEC
25
200
10
CFVO
20
200
6
N
30
200
10
pH
32
200
9
SOC
40
200
12
Texture ALR I
25
150
10
Texture ALR II
27
150
10
A key issue for DSM applications using legacy soil data is
the evaluation of the results. The aspect of evaluation which
compares actual with predicted values is numerical (or sta-
tistical) evaluation, often termed “validation” in the DSM lit-
erature; Oreskes (1998) and Rossiter (2017) explain why the
term “evaluation” is preferred for the overall process of as-
sessing the success of models, including DSM models. The
best approach to numerical evaluation is to have an inde-
pendent dataset obtained with probability sampling using a
known sampling design (Brus et al., 2011; Brus, 2014). How-
ever, this is not feasible when only legacy data are available.
In this case, a so-called “cross-validation” approach is often
used. This needs to be tuned to avoid over- or underestima-
tion of the numeric evaluation metrics, especially in case of
large differences in observation density, i.e. clustered spatial
observations. This is especially important at global scale, as
the distribution of the soil observations is not uniform across
the globe. It can not be guaranteed that the numeric evalua-
tion metrics derived from cross-validation are unbiased esti-
mates of the true numeric evaluation metrics, i.e. those that
would have been computed on a probability sample of the
whole population. It is also not possible to quantify how close
the cross-validation metrics estimates are to the true metrics,
as it is not possible to obtain confidence intervals (Brus et al.,
2011). When using cross-validation it is important to pre-
vent over- or under-optimistic estimates. For example, it is
likely that prediction errors are smaller in areas where the
sampling density is higher. Because of their high sampling
density, such areas will be over-represented in the sample
as the percentage of cross-validation points in clustered ar-
eas will be higher than the percentage of the total land area
covered by those areas. Results of standard cross-validation
will be strongly influenced by the performances in clustered
https://doi.org/10.5194/soil-7-217-2021
SOIL, 7, 217–240, 2021
226
L. Poggio et al.: Soil information for the globe
Figure 4.
Number of observations per grid cell (70 000 km
2
) for coarse fragments.
Figure 5.
Example of loss function (RMSE) used in the RFE step of covariates’ selection.
areas. Using spatial cross-validation as suggested by Meyer
et al. (e.g. 2018), where it is ensured that calibration data are
never too close to a cross-validation point, on the other hand
could produce over-pessimistic results. In order to address
some of these concerns, this study adopted a practical solu-
tion in which the folds were created to guarantee a spatially
balanced distribution between cross-validation folds, main-
taining the same densities of the input data in each fold so
that they represent approximately the same population.
Although the numerical evaluation procedure used in this
work takes into account the spatial distribution of the ob-
servations and their density, further improvement is possi-
ble in both model training and evaluation. For example, the
weights assigned to observations in heavily sampled areas
could be reduced. The United States and large regions of Eu-
rope and Australia have high numbers of observations that
could be downweighted to strengthen the spatial robustness
of the evaluation procedure. Declustering or debiasing tech-
niques (Goovaerts, 1997; Deutsch and Journel, 1998) have
been applied with success in other geostatistics exercises and
could be adapted to global soil mapping. The creation of the
folds could also be modified to take into account the density
of the observations.
3.4
Qualitative evaluation of spatial patterns
At global scale, well-known patterns are reproduced, and
typical properties associated with many World Reference
Base for Soil Resources (WRB) (IUSS Working Group
WRB, 2015) Reference Soil Groups can be recognised.
For example, the pH map identifies the large regions of
alkaline soils (Solonetz, Solochak), highly weathered soils
(e.g. Acrisols, Alisols, Plithosols), acid forest soils (e.g. Pod-
zols) and young soils from calcareous glacial deposits (e.g.
Luvisols). The low pH of Andosols (e.g. Pacific North-
west United States, Japan, New Zealand) is also correctly
represented. The texture components (particle size class –
PSC) maps correctly identify the siltier deltas (e.g. Yellow–
SOIL, 7, 217–240, 2021
https://doi.org/10.5194/soil-7-217-2021
L. Poggio et al.: Soil information for the globe
227
Table 4.
Global cross-validation results for both mean and median predictions. See Table 1 for abbreviations and units of the soil properties
considered.
Property
RMSE (median)
RMSE (mean)
MEC (median)
MEC (mean)
BDOD
0.19
0.19
0.73
0.74
CEC
11.01
10.69
0.40
0.43
CFVO
13.46
12.69
0.22
0.31
N
2.62
2.50
0.47
0.52
pH
0.78
0.77
0.67
0.68
SOC
39.67
36.48
0.37
0.47
Sand
0.19
0.18
0.51
0.54
Silt
0.13
0.13
0.60
0.62
Clay
0.13
0.13
0.42
0.43
Table 5.
MEC per depth layer for mean predictions. See Table 1 for abbreviations and units of the soil properties considered.
Depth layer
BDOD
CEC
CFVO
N
pH
SOC
Sand
Silt
Clay
0–5 cm
0.78
0.46
0.33
0.65
0.69
0.55
0.59
0.71
0.45
5–15 cm
0.74
0.42
0.35
0.41
0.66
0.39
0.58
0.64
0.42
15–30 cm
0.72
0.39
0.33
0.44
0.68
0.38
0.57
0.68
0.42
30–60 cm
0.70
0.42
0.31
0.46
0.68
0.38
0.54
0.62
0.41
60–100 cm
0.61
0.41
0.29
0.48
0.68
0.42
0.50
0.57
0.40
100–200 cm
0.59
0.45
0.29
0.49
0.67
0.59
0.48
0.54
0.40
Yangtze, Ganges–Brahmaputra), broad river plains (e.g. Po,
Danube, Mississippi, Rio Plate, upper Amazon), the loess re-
gions (e.g. Midwestern United States, NW Europe, Ukraine)
and the sandy North German–Polish plain. The cation ex-
change capacity (CEC) map clearly identifies large regions of
highly weathered clays (e.g. Southeastern United States and
China, central Brazil) and high-CEC 2
:
1 clays (e.g. “black
cotton” Vertisols in the Deccan plateau and the Sudan). This
map, together with the soil organic carbon (SOC) concentra-
tion maps, identifies large regions of Histosols (e.g. north-
ern Canada, Scotland, Siberia, Borneo). The SOC stock map
identifies deep Histosols and cool, wet regions (e.g. Pacific
Northwest North American coast, Ireland, southern Chile).
The coarse fragment map identifies large areas of the Tibetan
Plateau and the principal mountain chains, as well as recently
glaciated soils on igneous bedrock (e.g. Scandinavia, north-
ern Quebec and Ontario).
Many regional patterns are also clear, for example the pH
transition from Sahara through Sahel to the West African
coast and the PSC transitions from the Des Moines glacial
lobe to the proglacial loess deposits in Iowa (United States)
as well as the PSC transition from clayey marine sediments
along the North and Baltic Sea coasts through the sandy
plains to the central German loess belt. The CEC map iden-
tifies contrasting areas of Vertisols (e.g, “black belts” in
Alabama–Mississippi and Texas, United States). The coarse
fragment map shows the detailed pattern of the basin-and-
range region of the western United States and the ridge-and-
valley region of Appalachia.
However, at the local scale, a preliminary assessment of
SoilGrids in the United States, compared with a gridded ver-
sion of the national detailed gSSURGO (NRCS National Soil
Survey Center, 2016) soil geographic database based on a de-
tailed field survey, reveals that SoilGrids may fail to account
for local parent material transitions, e.g. sedimentary facies
of coastal plain marine sediments, as well as glacial features
such as proglacial lacustrine sediments and relic beach lines,
so that the local PSC pattern is not accurate, sometimes on
the order of 20 %–30 % of a particle-size class.
For example, Fig. 6a shows predicted sand concentration
of the 0–5 cm layer in an approximately 50
×
50 km area in
central Pennsylvania (United States), ranging from approxi-
mately 25 % (darkest red) to 80 % (darkest blue). Important
local differences are clear: the low sand concentrations of the
clayey soils in the limestone valleys trending SW–NE and
the high concentrations in soils from glacial till developed
on sandstones in the north, as well as the residual soils on the
resistant sandstone ridges of the Ridge and Valley province
of Appalachia in the south. These do generally agree with the
detailed soil survey.
At the detailed scale (250 m pixel), SoilGrids typically
shows fine details that do not always appear to be related
to obvious landscape or land use differences, when the map
is viewed as a ground overlay in Google Earth. For example,
Fig. 6b shows detail of predicted sand concentration of the 0–
5 cm layer in an approximately 3
×
3 km area of the previous
figure. The effect of some covariates being at 1 km resolu-
tion and others at 250 m is apparent, but the reason for the
https://doi.org/10.5194/soil-7-217-2021
SOIL, 7, 217–240, 2021
228
L. Poggio et al.: Soil information for the globe
Table 6.
Prediction interval coverage probability, global and by predicted depth interval. See Table 1 for abbreviations and units of the soil
properties considered.
Property
Global
[0, 5]
[5, 15]
[15, 30]
[30, 60]
[60, 100]
[100, 200]
BDOD
0.90
0.89
0.91
0.91
0.91
0.90
0.88
CEC
0.88
0.89
0.90
0.88
0.88
0.88
0.87
CFVO
0.95
0.96
0.95
0.95
0.95
0.94
0.94
N
0.92
0.91
0.92
0.93
0.92
0.92
0.92
pH
0.90
0.91
0.91
0.90
0.91
0.90
0.89
SOC
0.92
0.91
0.92
0.92
0.92
0.92
0.92
Sand
0.79
0.82
0.82
0.80
0.78
0.78
0.78
Silt
0.96
0.95
0.96
0.96
0.96
0.96
0.96
Clay
0.96
0.96
0.96
0.96
0.95
0.95
0.96
Figure 6.
Predicted sand concentration, %, 0–5 cm, ground overlay in © Google Earth.
(a)
Overview; centre
≈ −
77
◦
14
0
E, 41
◦
14
0
N, near
Jersey Shore, PA.
(b)
Detail; centre
≈ −
76
◦
56
0
E, 41
◦
33
0
N.
fine-scale differences is not. This area is of similar lithology,
relief and land cover (second-growth dense forest) except the
narrow valley at the north-west edge, yet the predictions are
quite different.
In this context, it should be realised that SoilGrids250m
predictions are not meant for use at a detailed scale, i.e. at
the subnational or local level, as national data providers of-
ten have access to more detailed point datasets and covariate
layers for their country than have been provided to the point
dataset on which SoilGrids250m is based (Chen et al., 2020;
Roudier et al., 2020; Vitharana et al., 2019; Liu et al., 2020).
3.5
Prediction uncertainty
In general, the least sampled areas present the highest pre-
diction uncertainties as expressed by the PICP. Figures 7
and 8 show an example for two properties and depths (maps
for all properties and depths can be accessed at https://data.
isric.org). Figure 9 and 10 show an example representing the
quantiles for pH
water
for the 60–100 cm layer. The north of
Russia and the centre and north-west of Canada are large re-
gions for which few soil observations are available; therefore
prediction distributions are wider than in more densely sam-
pled areas. However, these patterns are different for different
properties. For example, arid areas actually have the narrow-
est prediction ranges of pH
water
. The uncertainty range is of-
SOIL, 7, 217–240, 2021
https://doi.org/10.5194/soil-7-217-2021
L. Poggio et al.: Soil information for the globe
229
ten wide for properties and regions with a wider range of the
property being modelled. This can be explained by the mod-
elling approach performing more accurately within a limited
range of options. These regions also have larger local spatial
variation with more difficulties for predictions.
The communication of uncertainty is an open challenge
(Arrouays et al., 2020). Uncertainty should provide informa-
tion for policymakers and other stakeholders and not only
scientists and modellers. The maps computed with Eq. (3) are
a first step in this direction, but their limitations must be un-
derstood. For properties that have values at or near zero, e.g.
coarse fragments, they do not provide an entirely accurate
uncertainty estimate. The use of uncertainty classes could be
a further step to help domain stakeholders.
3.6
Limitations and outlook
This study represents a considerable effort to provide a glob-
ally consistent product using the point dataset available to IS-
RIC, a large number of relevant covariates and some optimi-
sation of a well-established machine learning method, within
the limits of practical computation. Yet it is clear that this
product has some limitations, which will be considered in
further work.
First, there is an ever-expanding group of new covari-
ates that can help explain and model the spatial variation
of soil properties. Products derived from Earth observation
are particularly relevant in this regard and have considerably
improved over the last decade. For example, the European
Space Agency Sentinel missions (both optical and radar) pro-
vide high-resolution data that have been shown to improve
DSM model performances.
Second, a fundamental problem is a lack of well dis-
tributed point observations within the soil property geo-
graphic and features space. Additional soil data for so far
under-represented regions, for example the northern boreal
regions as being collated by the International Soil Carbon
Network (Malhotra et al., 2019), will be sought for possible
consideration in the WoSIS workflow that provides the point
data underpinning the SoilGrids mapping effort. This effort
would be aided by the provision by more data providers of at
least a representative part of their point data to WoSIS, under
suitable license. It is also important to consider the distribu-
tion of the observations in the covariate space to minimise the
issues related to predictions into unknown regions of feature
space (Meyer and Pebesma, 2020).
Third, DSM methods are under active development, both
new methods and improvements to established methods. The
use of decision-tree-based models in DSM has become fairly
common in recent years. Models such as random forests,
XGBoost or Cubist tend to provide better results than most
multiple linear regression methods with reasonable computa-
tion costs (Khaledian and Miller, 2020). However, methods
such as artificial neural networks promise further improve-
ments in model performances if the amount and distribution
of the data support these highly complex models. This is the
case in particular with convolutional or recursive neural net-
works (deep learning). However, these methods present com-
putational challenges with the amount of training data nec-
essary for a sufficiently accurate DSM exercise, especially
when working at global scale at medium to fine resolutions.
Fourth, the proper method of cross-validation is an-
other important aspect when considering how to assess and
improve model performances. In particular, spatial cross-
validation and declustering of the data need to be further ex-
plored.
Fifth, this research considered only the modelling of some
primary soil properties, as defined and described in the Glob-
alSoilMap specifications. More work is necessary to obtain
maps for soil thickness (either rooting zone, pedogenetic
solum or regolith), soil properties derived with pedo-transfer
functions, e.g. hydrologic soil properties such as saturated
hydraulic conductivity (Pachepsky and Rawls, 2004), and
complex properties that depend on multiple primary prop-
erties, e.g. carbon stocks. These layers are important inputs
to model and map soil functions in the present and in the fu-
ture as well as to support Earth system modelling (Luo et al.,
2016; Dai et al., 2019).
Sixth, the quantification of uncertainty is recommended
and is becoming more common in DSM studies. This work
introduced it at global scale for the first time to the best of our
knowledge. While the provision of quantiles is mentioned
in the GlobalSoilMap specifications, the representation and
communication of uncertainty to end users and stakehold-
ers remain an important research field to be further explored.
The appropriate uncertainty intervals, both in terms of user
acceptance and modelling feasibility, also need to be investi-
gated.
Finally, the integration of highly automatised workflows
with expert opinion should be further explored. DSM prod-
ucts use statistical models to describe soils, and it is im-
portant to take into account the expertise and experience of
pedologists, at least in an evaluation loop if not as part of the
modelling itself. We made a first attempt at this in the Qual-
itative evaluation section above but do not have a method to
effectively incorporate expert observations into a workflow.
4
Conclusions
This study presents and discusses the production of global
maps of soil properties as implemented in the SoilGrids
2.0 product, with cross-validation, hyper-parameter selection
and quantification of uncertainty, using the best available
(shared) soil profile data for the world. In particular, the study
describes a robust and reproducible DSM workflow address-
ing the challenges of global data modelling:
1. non-homogeneous spatial distribution of input soil ob-
servations;
https://doi.org/10.5194/soil-7-217-2021
SOIL, 7, 217–240, 2021
230
L. Poggio et al.: Soil information for the globe
Figure 7.
Mean soil organic carbon content (dg/kg) prediction and range between 5 % and 95 % quantiles in the 5 to 15 cm depth interval,
(a)
for prediction and
(b)
for interquartile range.
Figure 8.
Median total nitrogen prediction (cg/kg) and associated uncertainty for the 15 to 30 cm depth interval,
(a)
for prediction and
(b)
for
uncertainty.
SOIL, 7, 217–240, 2021
https://doi.org/10.5194/soil-7-217-2021
L. Poggio et al.: Soil information for the globe
231
Figure 9.
Prediction distribution for pH
water
(10 pH) in the 60 to 100 cm depth interval,
(a)
for mean and
(b)
for median.
Figure 10.
Prediction distribution for pH
water
(10 pH) in the 60 to 100 cm depth interval,
(a)
for 5 % quantile and
(b)
for 95 % quantile.
https://doi.org/10.5194/soil-7-217-2021
SOIL, 7, 217–240, 2021
232
L. Poggio et al.: Soil information for the globe
2. robust quantitative evaluation with a cross-validation
procedure balancing accuracy and performances;
3. qualitative evaluation of the spatial patterns of the maps
to include information about matching with well recog-
nised pedo-landscape features;
4. quantification and mapping of the spatial uncertainty to
provide users with a measure for and warning for users
of the products.
As such, it describes a next step into global modelling and
mapping of soil properties, explicitly highlighting the im-
portance of quantitative and qualitative evaluation and un-
certainty communication. The actual use of SoilGrids 2.0 in
global and wide-area regional applications, where soil prop-
erties are important model inputs, will be the real test of its
applicability and usefulness.
SOIL, 7, 217–240, 2021
https://doi.org/10.5194/soil-7-217-2021
L. Poggio et al.: Soil information for the globe
233
Appendix A: Environmental covariates
Over 400 geographic layers were available as environmental
covariates for this work. These are chosen for their presumed
relation to the major soil forming factors.
Table A1.
Covariate sets and sources.
Weather and climate
- Temperature and precipitation from the Climatologies at high resolution for the Earth’s land surface areas (CHELSEA)
dataset (Karger et al., 2016)
- Snowfall from ESA’s CCI Land Cover dataset (Bontemps et al., 2013)
- Cloud cover by EarthEnv (Wilson and Jetz, 2016)
- Temperature and water vapour from NASA’s MODIS products (Wan, 2006)
- Precipitation, solar radiation, temperature, water vapour, wind speed plus various indexes from the WorldClim version 2
climate data series (Fick and Hijmans, 2017)
Ecology and ecosystems
- Bioclimatic zones in the Global Ecophysiography product by the USGS Geosciences and Environmental Change Science
Center (GECSC) (Dinerstein et al., 2017)
Geology
- Average soil and sedimentary-deposit thickness by the Distributed Active Archive Centre (DAAC) (Pelletier et al., 2016).
- Rock types by the USGS Geosciences and Environmental Change Science Center (GECSC), based on the Global Lithological
Map database v1.1 (Hartmann and Moosdorf, 2012)
Land use and land cover
- 2010 land cover classes from ESA’s CCI Land Cover (Bontemps et al., 2013)
- Bare ground and tree cover from the USGU Global Land Cover dataset (Hansen et al., 2013)
- 2010 land cover classes from the NGCC’s GLobeLand30 product (Chen et al., 2015)
Elevation and morphology
- Land surface elevation from the EarthEnv-DEM90 dataset (Robinson et al., 2014)
- Land surface elevation and various morphology indexes from the WorldGrids dataset (Reuter and Hengl, 2012)
- Land form classes in the Global Ecophysiography product by the USGS Geosciences and Environmental Change Science
Center (GECSC) (Sayre et al., 2014)
Core satellite outputs
- Bands 3, 4, 5 and 7 from Landsat (Zanter, 2019)
- Middle- and near-infrared bands from MODIS (Savtchenko et al., 2004)
Vegetation indexes
- NDVI from Landsat (Zanter, 2019)
- EVI and NPP from MODIS (Savtchenko et al., 2004)
Hydrography
- Global Inundation Extent from Multi-Satellites (GIEMS) dataset by Estellus (Fluet-Chouinard et al., 2015)
- Extent of glaciers, surface water change and occurrence probability by the JRC (Pekel et al., 2016)
- Global water table depth (Fan et al., 2013)
https://doi.org/10.5194/soil-7-217-2021
SOIL, 7, 217–240, 2021
234
L. Poggio et al.: Soil information for the globe
Appendix B: Bioclimatic regions
Table B1 summarises the number of observations per prop-
erty for each bioclimatic region. An interactive map of the re-
gions is available at http://ecoregions2017.appspot.com/ (last
access: 20 May 2021).
Table B1.
Number of observations per property for each bioclimatic region. See Table 1 for abbreviations.
Biome
CEC
CFVO
N
pH
SOC
STF
Tropical and subtropical moist broadleaf forests
4185
2117
8378
12 872
11 901
11 651
Tropical and subtropical dry broadleaf forests
558
205
1370
3264
2724
3051
Tropical and subtropical coniferous forests
59
30
54
1336
878
1331
Temperate broadleaf and mixed forests
12 585
29 708
24 711
56 569
49 727
61 822
Temperate conifer forests
6058
6417
5812
7597
9490
9834
Boreal forests/taiga
1443
3210
4834
4140
6819
5358
Tropical and subtropical grasslands, savannas and shrublands
8391
8259
20 181
27 633
24 951
23 135
Temperate grasslands, savannas and shrublands
13 442
9885
9812
23 654
24 421
25 416
Flooded grasslands and savannas
246
124
503
754
818
798
Montane grasslands and shrublands
479
1865
1073
1386
3994
3568
Tundra
312
199
466
548
807
695
Mediterranean forests, woodlands and scrub
1747
5951
8034
9126
12 532
11 428
Deserts and Xeric shrublands
3342
3412
3224
8994
8163
9862
Mangroves
88
26
165
264
437
250
SOIL, 7, 217–240, 2021
https://doi.org/10.5194/soil-7-217-2021
L. Poggio et al.: Soil information for the globe
235
Code and data availability.
The code underpinning the Soil-
Grids 2.0 workflow is available under the GPL3 license at https:
//git.wur.nl/isric/soilgrids/soilgrids (last access: 21 May 2021).
SoilGrids predictions themselves are available to the public un-
der the Creative Commons CC-BY 4.0 licence, facilitating their
widespread use. They may be obtained as world mosaics in the Vir-
tual Raster Tile (VRT) format from https://files.isric.org/soilgrids/
latest/ (last access: 21 May 2021). The Web Coverage Service
(WCS; https://maps.isric.org, last access: 21 May 2021) facilitates
automated access, e.g. from computer programmes or modelling
frameworks. A set of notebooks (https://git.wur.nl/isric/soilgrids/
soilgrids.notebooks, last access: 21 May 2021) was developed with
examples for the use of the WCS. A new web-based portal (https:
//soilgrids.org/, last access: 21 May 2021) was also developed with
this release, providing users with a light and swift means to visualise
and explore the new predictions, making the best of state-of-the-art
technologies for the web. A ReST API in beta stage is also available
at https://rest.soilgrids.org/ (last access: 21 May 2021).
Supplement.
The supplement related to this article is available
online at: https://doi.org/10.5194/soil-7-217-2021-supplement.
Author contributions.
LP and LMdS conceived and executed the
research and wrote the paper. NHB, GBMH and BK gave sugges-
tions about the approach and contributed extensively to the paper.
DR designed and executed the qualitative evaluation and wrote the
corresponding sections in the paper. NHB and ER designed and cre-
ated the database of soil observations. All authors reviewed the pa-
per.
Competing interests.
The authors declare that they have no con-
flict of interest.
Acknowledgements.
We thank Rik van den Bosch (ISRIC Di-
rector) for the internal support to the SoilGrids project. We espe-
cially thank the organisations and experts (https://www.isric.org/
explore/wosis/wosis-contributing-institutions-and-experts, last ac-
cess: 21 May 2021) that provided soil point data for consideration
in WoSIS and SoilGrids. Laura Poggio is a member of a consor-
tium supported by LE STUDIUM Loire Valley Institute for Ad-
vanced Studies through its LE STUDIUM Research Consortium
Programme.
Financial support.
This work was supported by ISRIC core fund-
ing, with additional support from the European Union’s EU H2020
Research and Innovation Programme, grant agreement no. 774378
(https://www.circasa-project.eu, last access: 21 May 2021). ISRIC
– World Soil Information, legally registered as the International Soil
Reference and Information Centre, receives core funding from the
Netherlands Government.
Review statement.
This paper was edited by Nicolas P. A. Saby
and reviewed by Feng Liu and Dominique Arrouays.
References
Aitchison, J.: The statistical analysis of compositional data, Chap-
man & Hall, London, UK, 1986.
Akpa, S. I. C., Odeh, I. O. A., Bishop, T. F. A., and
Hartemink,
A.
E.:
Digital
Mapping
of
Soil
Particle-
Size
Fractions
for
Nigeria,
Soil
Sci.,
78,
1953–1966,
https://doi.org/10.2136/sssaj2014.05.0202, 2014.
Arrouays, D., Grundy, M. G., Hartemink, A. E., Hempel, J. W.,
Heuvelink, G. B., Hong, S. Y., Lagacherie, P., Lelyk, G.,
McBratney, A. B., McKenzie, N. J., de Lourdes Mendonca-
Santos, M., Minasny, B., Montanarella, L., Odeh, I. O.,
Sanchez, P. A., Thompson, J. A., and Zhang, G.-L.: Global-
SoilMap: Toward a Fine-Resolution Global Grid of Soil Prop-
erties, in: Advances in Agronomy, Academic Press, 93–134,
https://doi.org/10.1016/B978-0-12-800137-0.00003-0, 2014.
Arrouays, D., Leenaars, J. G. B., Richer-de Forges, A. C., Ad-
hikari, K., Ballabio, C., Greve, M., Grundy, M., Guerrero, E.,
Hempel, J., Hengl, T., Heuvelink, G., Batjes, N., Carvalho, E.,
Hartemink, A., Hewitt, A., Hong, S.-Y., Krasilnikov, P., La-
gacherie, P., Lelyk, G., Libohova, Z., Lilly, A., McBratney, A.,
McKenzie, N., Vasquez, G. M., Mulder, V. L., Minasny, B., Mon-
tanarella, L., Odeh, I., Padarian, J., Poggio, L., Roudier, P., Saby,
N., Savin, I., Searle, R., Solbovoy, V., Thompson, J., Smith, S.,
Sulaeman, Y., Vintila, R., Rossel, R. V., Wilson, P., Zhang, G.-L.,
Swerts, M., Oorts, K., Karklins, A., Feng, L., Ibelles Navarro,
A. R., Levin, A., Laktionova, T., Dell’Acqua, M., Suvannang,
N., Ruam, W., Prasad, J., Patil, N., Husnjak, S., Pásztor, L.,
Okx, J., Hallet, S., Keay, C., Farewell, T., Lilja, H., Juilleret,
J., Marx, S., Takata, Y., Kazuyuki, Y., Mansuy, N., Panagos,
P., Van Liedekerke, M., Skalsky, R., Sobocka, J., Kobza, J.,
Eftekhari, K., Alavipanah, S. K., Moussadek, R., Badraoui, M.,
Da Silva, M., Paterson, G., Gonçalves, M. D. C., Theocharopou-
los, S., Yemefack, M., Tedou, S., Vrscaj, B., Grob, U., Kozák,
J., Boruvka, L., Dobos, E., Taboada, M., Moretti, L., and Ro-
driguez, D.: Soil legacy data rescue via GlobalSoilMap and
other international and national initiatives, GeoResJ, 14, 1–19,
https://doi.org/10.1016/j.grj.2017.06.001, 2017.
Arrouays, D., McBratney, A., Bouma, J., Libohova, Z., de Forges,
A. C. R., Morgan, C. L., Roudier, P., Poggio, L., and Mulder,
V. L.: Impressions of digital soil maps: The good, the not so
good, and making them ever better, Geoderma Regional, 20,
e00255, https://doi.org/10.1016/j.geodrs.2020.e00255, 2020.
Ballabio, C., Panagos, P., and Monatanarella, L.: Mapping topsoil
physical properties at European scale using the LUCAS database,
Geoderma, 261, 110–123, 2016.
Banwart, S., Black, H., Cai, Z., Gicheru, P., Joosten, H., Victo-
ria, R., Milne, E., Noellemeyer, E., Pascual, U., Nziguheba,
G., Vargas, R., Bationo, A., Buschiazzo, D., de Brogniez, D.,
Melillo, J., Richter, D., Termansen, M., van Noordwijk, M., Go-
verse, T., Ballabio, C., Bhattacharyya, T., Goldhaber, M., Niko-
laidis, N., Zhao, Y., Funk, R., Duffy, C., Pan, G., la Scala, N.,
Gottschalk, P., Batjes, N., Six, J., van Wesemael, B., Stocking,
M., Bampa, F., Bernoux, M., Feller, C., Lemanceau, P., and Mon-
tanarella, L.: Benefits of soil carbon: report on the outcomes of
https://doi.org/10.5194/soil-7-217-2021
SOIL, 7, 217–240, 2021
236
L. Poggio et al.: Soil information for the globe
an international scientific committee on problems of the environ-
ment rapid assessment workshop, Carbon Manag., 5, 185–192,
https://doi.org/10.1080/17583004.2014.913380, 2014.
Barnes, M.: Aichi targets: Protect biodiversity, not just area, Nature,
526, 195–195, https://doi.org/10.1038/526195e, 2015.
Barnes, R., Sahr, K., Evenden, G., Johnson, A., and Warmerdam, F.:
dggridR: Discrete Global Grids for R, R package version 2.0.4,
available at: https://github.com/r-barnes/dggridR/ (last access: 21
May 2021), 2016.
Batjes, N.: ISRIC-WISE derived soil properties on a 5 by 5 arc-
minutes global grid (ver. 1.2), Report 2012/01, ISRIC – World
Soil Information, available at: http://www.isric.org/sites/default/
files/isric_report_2012_01.pdf (last access: 21 May 2021), 2012.
Batjes,
N.:
Harmonised
soil
property
values
for
broad-
scale
modelling
(WISE30sec)
with
estimates
of
global
soil
carbon
stocks,
Geoderma,
269,
61–68,
https://doi.org/10.1016/j.geoderma.2016.01.034 2016.
Batjes, N., Al-Adamat, R., Bhattacharyya, T., Bernoux, M.,
Cerri, C., Gicheru, P., Kamoni, P., Milne, E., Pal, D., and
Rawajfih, Z.: Preparation of consistent soil data sets for
SOC modelling purposes: secondary SOTER data sets for
four case study areas, Agr. Ecosyst. Environ., 122, 26–34,
https://doi.org/10.1016/j.agee.2007.01.005, 2007.
Batjes, N. H.: Technologically achievable soil organic carbon se-
questration in world croplands and grasslands, Land Degrad.
Dev., 30, 25–32, https://doi.org/10.1002/ldr.3209, 2019.
Batjes, N. H., Ribeiro, E., and van Oostrum, A.: Standardised
soil profile data to support global mapping and modelling
(WoSIS snapshot 2019), Earth Syst. Sci. Data, 12, 299–320,
https://doi.org/10.5194/essd-12-299-2020, 2020.
Bontemps, S., Defourny, P., Radoux, J., Van Bogaert, E., Lamarche,
C., Achard, F., Mayaux, P., Boettcher, M., Brockmann, C.,
Kirches, G., Zulkhe, M., Kalogirou, V., Seifert, F. M., and Arino,
O.: Consistent global land cover maps for climate modelling
communities: Current achievements of the ESA’s land cover CCI,
in: Proceedings of the ESA Living Planet Symposium, Edin-
burgh, Scotland, 9–13, 2013.
Borrelli, P., Robinson, D. A., Fleischer, L. R., Lugato, E., Bal-
labio, C., Alewell, C., Meusburger, K., Modugno, S., Schütt,
B., Ferro, V., Bagarello, V., Oost, K. V., Montanarella, L., and
Panagos, P.: An assessment of the global impact of 21st cen-
tury land use change on soil erosion, Nat. Commun., 8, 2013,
https://doi.org/10.1038/s41467-017-02142-7, 2017.
Bouma,
J.:
Engaging
Soil
Science
in
Transdisciplinary
Research
Facing
“Wicked”
Problems
in
the
Infor-
mation
Society,
Soil
Sci.
Soc.
Am.
J.,
79,
454–458,
https://doi.org/10.2136/sssaj2014.11.0470, 2015.
Brus, D. J.: Statistical sampling approaches for soil monitoring,
Eur. J. Soil Sci., 65, 779–791, https://doi.org/10.1111/ejss.12176,
2014.
Brus, D. J., Kempen, B., and Heuvelink, G.: Sampling for val-
idation of digital soil maps, Eur. J. Soil Sci., 62, 394–407,
https://doi.org/10.1111/j.1365-2389.2011.01364.x, 2011.
Buchhorn, M., Lesiv, M., Tsendbazar, N. E., Herold, M., Ber-
tels, L., and Smets, B.: Copernicus Global Land Cover
Layers – Collection 2, Remote Sens.-Basel, 12, 1044,
https://doi.org/10.3390/rs12061044, 2020.
Chen, J., Chen, J., Liao, A., Cao, X., Chen, L., Chen, X., He, C.,
Han, G., Peng, S., Lu, M., Zhang, W., Tong, X, and Mills, J.:
Global land cover mapping at 30 m resolution: A POK-based op-
erational approach, ISPRS J. Photogramm., 103, 7–27, 2015.
Chen,
S.,
Mulder,
V.
L.,
Heuvelink,
G.
B.
M.,
Pog-
gio, L., Caubet, M., Román Dobarco, M., Walter, C.,
and
Arrouays,
D.:
Model
averaging
for
mapping
top-
soil organic carbon in France, Geoderma, 366, 114237,
https://doi.org/10.1016/j.geoderma.2020.114237, 2020.
Cowie, A. L., Orr, B. J., Castillo Sanchez, V. M., Chasek, P., Cross-
man, N. D., Erlewein, A., Louwagie, G., Maron, M., Metter-
nicht, G. I., Minelli, S., Tengberg, A. E., Walter, S., and Wel-
ton, S.: Land in balance: The scientific conceptual framework
for Land Degradation Neutrality, Environ. Sci. Policy, 79, 25–
35, https://doi.org/10.1016/j.envsci.2017.10.011, 2018.
Dai, Y., Shangguan, W., Wei, N., Xin, Q., Yuan, H., Zhang, S.,
Liu, S., Lu, X., Wang, D., and Yan, F.: A review of the global
soil property maps for Earth system models, SOIL, 5, 137–158,
https://doi.org/10.5194/soil-5-137-2019, 2019.
de Sousa, L. M., Poggio, L., and Kempen, B.: Comparison
of FOSS4G Supported Equal-Area Projections Using Dis-
crete Distortion Indicatrices, ISPRS Int. J. Geo-Inf., 8, 351,
https://doi.org/10.3390/ijgi8080351, 2019.
de Sousa, L. M., Poggio, L., Dawes, G., Kempen, B., and
Van Den Bosch, R.: Computational Infrastructure of SoilGrids
2.0, in: International Symposium on Environmental Software
Systems, 24–31, 2020.
Deutsch, C. and Journel, A.: GSLIB: Geostatistical Software Li-
brary and User’s Guide, edn. 2, Oxford University Press, New
York, USA, 1998.
Dinerstein, E., Olson, D., Joshi, A., Vynne, C., Burgess, N. D.,
Wikramanayake, E., Hahn, N., Palminteri, S., Hedao, P., Noss,
R., Hansen, M., Locke, H., Ellis, E. C., Jones, B., Barber, C. V.,
Hayes, R., Kormos, C., Martin, V., Crist, E., Sechrest, W., Price,
L., Baillie, J. E. M., Weeden, D., Suckling, K., Davis, C., Sizer,
N., Moore, R., Thau, D., Birch, T., Potapov, P., Turubanova, S.,
Tyukavina, A., de Souza, N., Pintea, L., Brito, J. C., Llewellyn,
O. A., Miller, A. G., Patzelt, A., Ghazanfar, S. A., Timber-
lake, J., Klöser, H., Shennan-Farpón, Y., Kindt, R., Lillesø, J.-
P. B., van Breugel, P., Graudal, L., Voge, M., Al-Shammari,
K. F., and Saleem, M.: An Ecoregion-Based Approach to Pro-
tecting Half the Terrestrial Realm, BioScience, 67, 534–545,
https://doi.org/10.1093/biosci/bix014, 2017.
Dorji, T., Odeh, I. O. A., Field, D. J., and Baillie, I. C.: Dig-
ital soil mapping of soil organic carbon stocks under dif-
ferent land use and land cover types in montane ecosys-
tems, Eastern Himalayas, Forest Ecol. Manag., 318, 91–102,
https://doi.org/10.1016/j.foreco.2014.01.003, 2014.
Ellili, Y., Walter, C., Michot, D., Pichelin, P., and Lemercier, B.:
Mapping soil organic carbon stock change by soil monitoring
and digital soil mapping at the landscape scale, Geoderma, 351,
1–8, https://doi.org/10.1016/j.geoderma.2019.03.005, 2019.
Fan, Y., Li, H., and Miguez-Macho, G.: Global Patterns
of
Groundwater
Table
Depth,
Science,
339,
940–943,
https://doi.org/10.1126/science.1229881, 2013.
FAO: Digital Soil Map of the World and derived properties (ver.
3.5), Report FAO Land and Water Digital Media Series 1, Food
and Agriculture Organization of the United Nations (FAO), avail-
able at: http://www.fao.org/geonetwork/srv/en/metadata.show?
id=14116 (last access: 21 May 2021), 1995.
SOIL, 7, 217–240, 2021
https://doi.org/10.5194/soil-7-217-2021
L. Poggio et al.: Soil information for the globe
237
FAO: Guidelines for soil description, edn. 4, Report, Food and Agri-
culture Organization of the United Nations (FAO), available at:
http://www.fao.org/docrep/019/a0541e/a0541e.pdf (last access:
21 May 2021), 2006.
FAO and ITPS: Status of the world’s soil resources (SWSR) – Main
report, Report, Food and Agriculture Organization of the United
Nations (FAO) and Intergovernmental Technical Panel on Soils
(ITPS), available at: http://www.fao.org/3/a-i5199e.pdf (last ac-
cess: 21 May 2021), 2015.
FAO, IIASA, ISRIC, ISSCAS, and JRC: Harmonized World Soil
Database (version 1.2), Report, edited by: Nachtergaele, F.
O., van Velthuizen, H., Verelst, L., Wiberg, D., Batjes, N.
H., Dijkshoorn, J. A., van Engelen, V. W. P., Fischer, G.,
Jones, A., Montanarella, L., Petri, M., Prieler, S., Teixeira, E.,
and Xuezheng, S., Food and Agriculture Organization of the
United Nations (FAO), International Institute for Applied Sys-
tems Analysis (IIASA), ISRIC – World Soil Information, In-
stitute of Soil Science – Chinese Academy of Sciences (ISS-
CAS), Joint Research Centre of the European Commission
(JRC), available at: http://webarchive.iiasa.ac.at/Research/LUC/
External-World-soil-database/HWSD_Documentation.pdf (last
access: 21 May 2021), 2012.
FAO, IFAD, UNICEF, WFP, and WHO: The State of Food Se-
curity and Nutrition in the World 2018, Building climate re-
silience for food security and nutrition, Report, Food and Agri-
culture Organization of the United Nations (FAO), available at:
http://www.fao.org/3/I9553EN/i9553en.pdf (last access: 21 May
2021), 2018.
FAO-ISRIC: Guidelines for soil description, edn. 3, Food and Agri-
culture Organization of the United Nations (FAO), Rome, Italy,
1986.
FAO-Unesco: FAO-Unesco Soil Map of the World, 1
:
5 000 000
(Vol. 1 to 10), United Nations Educational, Scientific, and Cul-
tural Organization, Paris, France, 1971–1981.
Fick, S. E. and Hijmans, R. J.: WorldClim 2: new 1 km spatial reso-
lution climate surfaces for global land areas, Int. J. Climatol., 37,
4302–4315, 2017.
Fluet-Chouinard, E., Lehner, B., Rebelo, L.-M., Papa, F., and
Hamilton, S. K.: Development of a global inundation map at
high spatial resolution from topographic downscaling of coarse-
scale remote sensing data, Remote Sens. Environ., 158, 348–361,
2015.
Gomes, L. C., Faria, R. M., de Souza, E., Veloso, G. V., Schaefer,
C. E. G., and Fernandes Filho, E. I.: Modelling and mapping soil
organic carbon stocks in Brazil, Geoderma, 340, 337–350, 2019.
Goovaerts, P.: Geostatistics for Natural Resources Evaluation, Ox-
ford University Press, Oxford, UK, 500 pp., 1997.
GRASS Development Team: Geographic Resources Analysis Sup-
port System (GRASS GIS) Software, version 7.8.0, available at:
http://www.grass.osgeo.org (last access: 21 May 2021), 2020.
Grunwald, S., Thompson, J. A., and Boettinger, J. L.: Digital
soil mapping and modeling at continental scales: Finding solu-
tions for global issues, Soil Sci. Soc. Am. J., 75, 1201–1213,
https://doi.org/10.2136/sssaj2011.0025, 2011.
GSP and FAO: Pillar 4 Implementation Plan – Towards a
Global Soil Information System, Report, Global Soil Partner-
ship, available at: http://www.fao.org/3/a-bl102e.pdf (last access:
21 May 2021), 2016.
GSP and ITPS: Global soil organic carbon map (GSOCmap), Re-
port Technical Report, Global Soil Partnership (GSP) and Inter-
national Panel on Soils (ITPS), available at: http://www.fao.org/
3/i8891en/I8891EN.pdf (last access: 21 May 2021), 2018.
Guevara, M., Olmedo, G. F., Stell, E., Yigini, Y., Aguilar Duarte,
Y., Arellano Hernández, C., Arévalo, G. E., Arroyo-Cruz, C. E.,
Bolivar, A., Bunning, S., Bustamante Cañas, N., Cruz-Gaistardo,
C. O., Davila, F., Dell Acqua, M., Encina, A., Figueredo Tacona,
H., Fontes, F., Hernández Herrera, J. A., Ibelles Navarro, A. R.,
Loayza, V., Manueles, A. M., Mendoza Jara, F., Olivera, C., Os-
orio Hermosilla, R., Pereira, G., Prieto, P., Ramos, I. A., Rey
Brina, J. C., Rivera, R., Rodríguez-Rodríguez, J., Roopnarine,
R., Rosales Ibarra, A., Rosales Riveiro, K. A., Schulz, G. A.,
Spence, A., Vasques, G. M., Vargas, R. R., and Vargas, R.: No
silver bullet for digital soil mapping: country-specific soil or-
ganic carbon estimates across Latin America, SOIL, 4, 173–193,
https://doi.org/10.5194/soil-4-173-2018, 2018.
Guyon, I., Weston, J., Barnhill, S., and Vapnik, V.: Gene selection
for cancer classification using support vector machines, Mach.
Learn., 46, 389–422, 2002.
Han, E., Ines, A. V. M., and Koo, J.: Development of a
10 km resolution global soil profile dataset for crop mod-
eling applications, Environ. Modell. Softw., 119, 70–83,
https://doi.org/10.1016/j.envsoft.2019.05.012, 2019.
Hansen, M. C., Potapov, P. V., Moore, R., Hancher, M., Turubanova,
S., Tyukavina, A., Thau, D., Stehman, S., Goetz, S., Loveland,
T. R., Kommareddy, A., Egorov, A., Chini, L., Justice, C. O.,
and Townshend, J. R.,G.: High-resolution global maps of 21st-
century forest cover change, Science, 342, 850–853, 2013.
Harden, J. W., Hugelius, G., Ahlström, A., Blankinship, J. C., Bond-
Lamberty, B., Lawrence, C. R., Loisel, J., Malhotra, A., Jackson,
R. B., Ogle, S., Phillips, C., Ryals, R., Todd-Brown, K., Var-
gas, R., Vergara, S. E., Cotrufo, M. F., Keiluweit, M., Heckman,
K. A., Crow, S. E., Silver, W. L., DeLonge, M., and Nave, L. E.:
Networking our science to characterize the state, vulnerabilities,
and management opportunities of soil organic matter, Global
Change Biol., 24, 705–718, https://doi.org/10.1111/gcb.13896,
2017.
Hartmann, J. and Moosdorf, N.: The new global lithologi-
cal map database GLiM: A representation of rock proper-
ties at the Earth surface, Geochem. Geophy. Geosy., 13,
https://doi.org/10.1029/2012GC004370, 2012.
Hengl, T., Mendes de Jesus, J., MacMillan, R. A., Bat-
jes, N. H., Heuvelink, G. B., Ribeiro, E. C., Samuel-
Rosa, A., Kempen, B., Leenaars, J. G., Walsh, M. G.,
and Gonzalez, M. R.: SoilGrids1km – global soil informa-
tion based on automated mapping, PLoS ONE, 9, e114788,
https://doi.org/10.1371/journal.pone.0105992, 2014.
Hengl, T., Leenaars, J. G. B., Shepherd, K. D., Walsh, M. G.,
Heuvelink, G. B. M., Mamo, T., Tilahun, H., Berkhout,
E., Cooper, M., Fegraus, E., Wheeler, I., and Kwabena,
N. A.: Soil nutrient maps of Sub-Saharan Africa: assess-
ment of soil nutrient content at 250 m spatial resolution us-
ing machine learning, Nutr. Cycl. Agroecosys., 109, 77–102,
https://doi.org/10.1007/s10705-017-9870-x, 2017a.
Hengl, T., Mendes de Jesus, J., Heuvelink, G. B. M., Ruiperez Gon-
zalez, M., Kilibarda, M., Blagoti´c, A., Shangguan, W.,
Wright, M. N., Geng, X., Bauer-Marschallinger, B., Gue-
vara, M. A., Vargas, R., MacMillan, R. A., Batjes, N. H.,
https://doi.org/10.5194/soil-7-217-2021
SOIL, 7, 217–240, 2021
238
L. Poggio et al.: Soil information for the globe
Leenaars, J. G. B., Ribeiro, E., Wheeler, I., Mantel, S., and
Kempen, B.: SoilGrids250m: Global gridded soil informa-
tion based on machine learning, PLoS ONE, 12, e0169748,
https://doi.org/10.1371/journal.pone.0169748, 2017b.
Heuvelink, G. B. M., Angelini, M. E., Poggio, L., Bai, Z., Batjes,
N. H., van den Bosch, R., Bossio, D., Estella, S., Lehmann, J.,
Olmedo, G. F., and Sanderman, J.: Machine learning in space
and time for modelling soil organic carbon change, Eur. J. Soil
Sci., 1–17, https://doi.org/10.1111/ejss.12998, 2020.
Hounkpatin, O. K., de Hipt, F. O., Bossa, A. Y., Welp, G., and
Amelung, W.: Soil organic carbon stocks and their determining
factors in the Dano catchment (Southwest Burkina Faso), Catena,
166, 298–309, 2018.
IPBES: Global assessment report on biodiversity and ecosystem
services of the Intergovernmental Science – Policy Platform
on Biodiversity and Ecosystem Services, edited by: Brondizio,
E. S., Settele, J., Díaz, S., and Ngo, H. T., Report, IPBES
Secretariat, Bonn, Germany, available at: https://www.ipbes.net/
global-assessment-report-biodiversity-ecosystem-services (last
access: 21 May 2021), 2019.
IUSS Working Group WRB: World Reference Base for Soil Re-
sources 2014, in: Update 2015 – International Soil Classification
System for Naming Soils and Creating Legends for Soil Maps,
no. 106 in World Soil Resources Reports, Food and Agriculture
Organization of the United Nations (FAO), Rome, Italy, avail-
able at: http://www.fao.org/3/i3794en/I3794en.pdf (last access:
21 May 2021), 2015.
Ivushkin, K., Bartholomeus, H., Bregt, A. K., Pulatov, A.,
Kempen, B., and de Sousa, L.: Global mapping of soil
salinity
change,
Remote
Sens.
Environ.,
231,
11260,
https://doi.org/10.1016/j.rse.2019.111260, 2019.
Janssen, P. H. M. and Heuberger, P. S. C.: Calibration of process-
oriented models, Ecol. Model., 83, 55–66, 1995.
Karger, D., Conrad, O., Böhner, J., Kawohl, T., Kreft, H., Soria-
Auza, R., Zimmermann, N., Linder, H., and Kessler, M.:
CHELSA climatologies at high resolution for the earth’s land
surface areas (Version 1.1), World Data Center for Climate, 2016.
Kempen, B., Brus, D. J., and de Vries, F.: Operationalizing dig-
ital soil mapping for nationwide updating of the 1
:
50 000
soil map of the Netherlands, Geoderma, 241–242, 313–329,
https://doi.org/10.1016/j.geoderma.2014.11.030, 2015.
Kempen, B., Dalsgaard, S., Kaaya, A. K., Chamuya, N.,
Ruipérez-González,
M.,
Pekkarinen,
A.,
and
Walsh,
M.
G.:
Mapping
topsoil
organic
carbon
concentrations
and
stocks
for
Tanzania,
Geoderma,
337,
164–180,
https://doi.org/10.1016/j.geoderma.2018.09.011, 2019.
Keskin, H. and Grunwald, S.: Regression kriging as a workhorse
in the digital soil mapper’s toolbox, Geoderma, 326, 22–41,
https://doi.org/10.1016/j.geoderma.2018.04.004, 2018.
Khaledian, Y. and Miller, B. A.: Selecting appropriate machine
learning methods for digital soil mapping, Appl. Math. Model.,
81, 401–418, 2020.
Kuhn, M.: A Short Introduction to the caret Package, R Found Stat
Comput, https://cran.microsoft.com/snapshot/2015-08-17/web/
packages/caret/vignettes/caret.pdf (last access: 21 May 2021),
10 pp., 2015.
Lark, R. and Bishop, T.: Cokriging particle size fractions of the soil,
Eur. J. Soil Sci., 58, 763–774, 2007.
Liu, F., Zhang, G.-L., Song, X., Li, D., Zhao, Y., Yang, J.,
Wu, H., and Yang, F.: High-Resolution and Three-Dimensional
Mapping of Soil Texture of China, Geoderma, 361, 114061,
https://doi.org/10.1016/j.geoderma.2019.114061, 2020.
Luo, Y., Ahlström, A., Allison, S., Batjes, N., Brovkin, V., Car-
valhais, N., Chappell, A., Ciais, P., Davidson, E., Finzi, A.,
Georgiou, K., Guenet, B., Hararuk, O., Harden, J., He, Y.,
Hopkins, F. M., Jiang, L., Koven, C., Jackson, R., Jones, C.,
Lara, M., Liang, J., McGuire, A., Parton, W. J., Peng, C., Ran-
derson, J., Salazr, A., Sierra, C., Smoth, M., Tian, H., Todd-
Brown, K., Torn, M., van Groeningen, K., Wang, Y. P., Westm,
O., Wei, Y., Wieder, W., Xia, J., Xia, X., Xu, X., and Zhu,
T.: Towards more realistic projections of soil carbon dynamics
by Earth System Models, Global Biogeochem. Cy., 30, 40–56,
https://doi.org/10.1002/2015GB005239, 2016.
Ma, Y., Minasny, B., McBratney, A., Poggio, L., and Fa-
jardo,
M.:
Predicting
soil
properties
in
3D:
Should
depth
be
a
covariate?,
Geoderma,
383,
114794,
https://doi.org/10.1016/j.geoderma.2020.114794, 2021.
Malhotra, A., Todd-Brown, K., Nave, L. E., Batjes, N. H.,
Holmquist, J. R., Hoyt, A. M., Iversen, C. M., Jack-
son, R. B., Lajtha, K., Lawrence, C., Vindušková, O.,
Wieder, W., Williams, M., Hugelius, G., and Harden, J.:
The landscape of soil carbon data: emerging questions,
synergies and databases, Prog. Phys. Geog., 43, 707–719,
https://doi.org/10.1177/0309133319873309, 2019.
Mallavan, B., Minasny, B., and Mcbratney, A.: Homosoil, a
Methodology for Quantitative Extrapolation of Soil Information
Across the Globe, in: Digital Soil Mapping, Progress in Soil Sci-
ence, edited by: Boettinger, J. L., Howell, D. W., Moore, A. C.,
Hartemink, A. E., and Kienast-Brown, S., Springer, Dordrecht,
The Netherlands, 137–150, https://doi.org/10.1007/978-90-481-
8863-5_12, 2010.
McBratney, A. B., Mendonça Santos, M. L., and Minasny,
B.:
On
digital
soil
mapping,
Geoderma,
117,
3–52,
https://doi.org/10.1016/S0016-7061(03)00223-4, 2003.
Meinshausen, N.: Quantile regression forests, J. Mach. Learn. Res.,
7, 983–999, 2006.
Meyer, H. and Pebesma, E.: Predicting into unknown space? Es-
timating the area of applicability of spatial prediction models,
ArXiv, abs/2005.07939, 2020.
Meyer, H., Reudenbach, C., Hengl, T., Katurji, M., and
Nauss, T.: Improving performance of spatio-temporal ma-
chine learning models using forward feature selection and
target-oriented validation, Environ. Modell. Softw., 101, 1–9,
https://doi.org/10.1016/j.envsoft.2017.12.001, 2018.
Minasny, B. and McBratney, A. B.: Digital soil mapping: A
brief history and some lessons, Geoderma, 264, 301–311,
https://doi.org/10.1016/j.geoderma.2015.07.017, 2016.
Mora-Vallejo,
A.,
Claessens,
L.,
Stoorvogel,
J.,
and
Heuvelink,
G.
B.
M.:
Small
scale
digital
soil
map-
ping
in
Southeastern
Kenya,
Catena,
76,
44–53,
https://doi.org/10.1016/j.catena.2008.09.008, 2008.
Moulatlet, G. M., Zuquim, G., Figueiredo, F. O. G., Lehtonen,
S., Emilio, T., Ruokolainen, K., and Tuomisto, H.: Using dig-
ital soil maps to infer edaphic affinities of plant species in
Amazonia: Problems and prospects, Ecol. Evol., 7, 8463–8477,
https://doi.org/10.1002/ece3.3242, 2017.
SOIL, 7, 217–240, 2021
https://doi.org/10.5194/soil-7-217-2021
L. Poggio et al.: Soil information for the globe
239
Nauman, T. W. and Duniway, M. C.: Relative prediction intervals
reveal larger uncertainty in 3D approaches to predictive digi-
tal soil mapping of soil properties with legacy data, Geoderma,
347, 170–184, https://doi.org/10.1016/j.geoderma.2019.03.037,
2019.
Nijbroek, R., Piikki, K., Söderström, M., Kempen, B., Turner, K.,
Hengari, S., and Mutua, J.: Soil Organic Carbon Baselines for
Land Degradation Neutrality: Map Accuracy and Cost Tradeoffs
with Respect to Complexity in Otjozondjupa, Namibia, Sustain-
ability, 10, 1610, https://doi.org/10.3390/su10051610, 2018.
NRCS National Soil Survey Center: Gridded Soil Survey
Geographic
(gSSURGO)
Database
User
Guide,
Version
2.2, Tech. Rep., National Soil Survey Center, available
at:
http://www.nrcs.usda.gov/wps/PA_NRCSConsumption/
download?cid=nrcs142p2_051847&ext=pdf (last access: 21
May 2021), 2016.
Nussbaum, M., Spiess, K., Baltensweiler, A., Grob, U., Keller, A.,
Greiner, L., Schaepman, M. E., and Papritz, A.: Evaluation of
digital soil mapping approaches with large sets of environmental
covariates, SOIL, 4, 1–22, https://doi.org/10.5194/soil-4-1-2018,
2018.
Omuto, C., Nachtergaele, F., and Vargas Rojas, R.: State of the Art
Report on Global and Regional Soil Information: Where are we?
Where to go?, Report, Food and Agriculture Organization of the
United Nations (FAO), available at: http://www.fao.org/docrep/
017/i3161e/i3161e.pdf (last access: 21 May 2021), 2012.
Oreskes,
N.:
Evaluation
(Not
Validation)
of
Quantita-
tive
Models,
Environ.
Health
Persp.,
106,
1453–1460,
https://doi.org/10.1289/ehp.98106s61453, 1998.
Pachepsky, Y. and Rawls, W. J. (Eds.): Development of Pedotransfer
Functions in Soil Hydrology, in: Developments in Soil Science,
Elsevier, 2004.
Pekel,
J.-F.,
Cottam,
A.,
Gorelick,
N.,
and
Belward,
A. S.: High-resolution mapping of global surface wa-
ter
and
its
long-term
changes,
Nature,
540,
418–422,
https://doi.org/10.1038/nature20584, 2016.
Pelletier, J., Broxton, P., Hazemberg, P., Zeng, X., Troch, P., Niu, G.,
Williams, Z., Brunke, M., and Gochis, D.: Global 1 km Gridded
Thickness of Soil, Regolith, and Sedimentary Deposit Layers,
https://doi.org/10.3334/ORNLDAAC/1304, 2016.
Piikki, K., Söderström, M., and Stadig, H.: Local adapta-
tion of a national digital soil map for use in precision
agriculture, Advances in Animal Biosciences, 8, 430–432,
https://doi.org/10.1017/s2040470017000966, 2017.
Poggio,
L.
and
Gimona,
A.:
3D
mapping
of
soil
tex-
ture
in
Scotland,
Geoderma
Regional,
9,
5–16,
https://doi.org/10.1016/j.geodrs.2016.11.003, 2017a.
Poggio, L. and Gimona, A.: Assimilation of optical and
radar remote sensing data in 3D mapping of soil proper-
ties over large areas, Sci. Total Environ., 579, 1094–1110,
https://doi.org/10.1016/j.scitotenv.2016.11.078, 2017b.
Poggio, L., Gimona, A., and Brewer, M. J.: Regional scale map-
ping of soil properties and their uncertainty with a large num-
ber of satellite-derived covariates, Geoderma, 209–210, 1–14,
https://doi.org/10.1016/j.geoderma.2013.05.029, 2013.
R Core Team: A Language and Environment for Statistical
Computing, R Foundation for Statistical Computing, Vienna,
Austria, available at: http://www.R-project.org/ (last access:
21 May 2021), 2020.
Reuter, H. I. and Hengl, T.: Worldgrids – a public repository of
global soil covariates, in: Digital Soil Assessments and Be-
yond – Proceedings of the 5th Global Workshop on Digital
Soil Mapping, 287–292, https://tinyurl.com/yyx7vd9s (last ac-
cess: 21 May 2021), 2012.
Ribeiro, E., Batjes, N., and Van Oostrum, A.: Towards the stan-
dardization and harmonization of world soil data, in: Pro-
cedures Manual Report 2018, World Soil Information Ser-
vice (WoSIS), Wageningen, The Netherlands, Report 2018/01,
https://doi.org/10.17027/isric-wdcsoils.20180001, 2018.
Robinson, N., Regetz, J., and Guralnick, R. P.: EarthEnv-DEM90: A
nearly-global, void-free, multi-scale smoothed, 90 m digital ele-
vation model from fused ASTER and SRTM data, ISPRS J. Pho-
togramm., 87, 57–67, 2014.
Rockstroem, J., Falkenmark, M., Lannerstad, M., and Karlberg,
L.: The planetary water drama: Dual task of feeding humanity
and curbing climate change, Geophys. Res. Lett., 39, L15401,
https://doi.org/10.1029/2012gl051688, 2012.
Rossiter, D. G.: Maps and Models Are Never Valid, but They Can
Be Evaluated (with Responses), Pedometron, 41, 19–28, avail-
able at: http://www.pedometrics.org/Pedometron/Pedometron41.
pdf (last access: 21 May 2021), 2017.
Roudier, P., Burge, O. R., Richardson, S. J., McCarthy, J. K., Greal-
ish, G., and Ausseil, A.-G.: National Scale 3D Mapping of Soil
pH Using a Data Augmentation Approach, Remote Sens.-Basel,
12, 2872, https://doi.org/10.3390/rs12182872, 2020.
Sanderman, J., Hengl, T., and Fiske, G. J.: Soil carbon debt of
12 000 years of human land use, P. Natl. Acad. Sci. USA, 114,
9575–9580, https://doi.org/10.1073/pnas.1706103114, 2017.
Savtchenko, A., Ouzounov, D., Ahmad, S., Acker, J., Leptoukh, G.,
Koziana, J., and Nickless, D.: Terra and Aqua MODIS products
available from NASA GES DAAC, Adv. Space Res., 34, 710–
714, 2004.
Sayre, R., Dangermond, J., Frye, C., Vaughan, R., Aniello, P.,
Breyer, S., Cribbs, D., Hopkins, D., Nauman, R., Derrenbacher,
W., Wright, D. J., Brown, C., Convis, C., Smith, J. H., Ben-
son, L., VanSistine, P., Warner, H., Cress, J. J., Danielson, J. J.,
Hamann, S. L., Cecere, T., Reddy, A. D., Burton, D., Grosse, A.,
True, D., Metzger, M., Hartmann, J., Moosdorf, N., Durr, H., Pa-
ganini, M., Defourny, P., Arino, O., Maynard, S., Anderson, M.,
and Comer, P.: A new map of global ecological land units – an
ecophysiographic stratification approach, Association of Ameri-
can Geographers, Washington D.C., USA, 2014.
Schoeneberger, P., Wysicki, D., Benham, E., and Broderson, W.:
Field book for describing and sampling soils (ver. 3.0), Natural
Resources Conservation Service, National Soil Survey Center,
Lincoln, Nebraska, USA, available at: http://www.nrcs.usda.gov/
Internet/FSE_DOCUMENTS/nrcs142p2_052523.pdf (last ac-
cess: 21 May 2021), 2012.
Shangguan, W., Dai, Y., Duan, Q., Liu, B., and Yuan, H.: A global
soil data set for earth system modeling, J. Adv. Model. Earth Sy.,
6, 249–263, https://doi.org/10.1002/2013MS000293, 2014.
Shrestha,
D.
L.
and
Solomatine,
D.
P.:
Machine
learn-
ing
approaches
for
estimation
of
prediction
interval
for the model output, Neural Networks, 19, 225–235,
https://doi.org/10.1016/j.neunet.2006.01.012, 2006.
Smith, P., Soussana, J.-F., Angers, D., Schipper, L., Chenu, C.,
Rasse, D., Batjes, N. H., van Egmond, F., McNeill, S., Kuhn-
ert, M., Arias-Navaro, C., Olesen, J. E., Chirinda, N., Fornara,
https://doi.org/10.5194/soil-7-217-2021
SOIL, 7, 217–240, 2021
240
L. Poggio et al.: Soil information for the globe
D., Joosse, P., Wollenberg, L., Alvaro-Fuentes, J., and Cobena,
A.: How to measure, report and verify soil carbon change to re-
alise the potential of soil carbon sequestration for atmospheric
greenhouse gas remova, Global Change Biol., 26, 219–241,
https://doi.org/10.1111/gcb.14815, 2019.
Soussana, J.-F., Lutfalla, S., Ehrhardt, F., Rosenstock, T., Lamanna,
C., Havlík, P., Richards, M., Wollenberg, E., Chotte, J.-L.,
Torquebiau, E., Ciais, P., Smith, P., and Lal, R.: Matching pol-
icy and science: Rationale for the “4 per 1000 – soils for
food security and climate” initiative, Soil Till. Res., 188, 3–15,
https://doi.org/10.1016/j.still.2017.12.002, 2017.
Springmann, M., Clark, M., Mason-D’Croz, D., Wiebe, K.,
Bodirsky, B. L., Lassaletta, L., de Vries, W., Vermeulen, S. J.,
Herrero, M., Carlson, K. M., Jonell, M., Troell, M., DeClerck,
F., Gordon, L. J., Zurayk, R., Scarborough, P., Rayner, M.,
Loken, B., Fanzo, J., Godfray, H. C. J., Tilman, D., Rock-
ström, J., and Willett, W.: Options for keeping the food
system within environmental limits, Nature, 562, 519–525,
https://doi.org/10.1038/s41586-018-0594-0, 2018.
Stockmann, U., Padarian, J., McBratney, A., Minasny, B., de Brog-
niez, D., Montanarella, L., Hong, S. Y., Rawlins, B. G., and Field,
D. J.: Global soil organic carbon assessment, Global Food Secu-
rity, 6, 9–16, https://doi.org/10.1016/j.gfs.2015.07.001, 2015.
Stoorvogel, J. J., Bakkenes, M., Temme, A. J. A. M., Bat-
jes, N. H., and ten Brink, B.: S-World: a Global Soil Map
for Environmental Modelling, Land Degrad. Dev., 28, 22–33,
https://doi.org/10.1002/ldr.2656, 2017.
Strobl, C., Boulesteix, A.-L., Kneib, T., Augustin, T., and Zeileis,
A.: Conditional variable importance for random forests, BMC
Bioinformatics, 9, 307, https://doi.org/10.1186/1471-2105-9-
307, 2008.
Todd-Brown, K. E. O., Randerson, J. T., Post, W. M., Hoffman, F.
M., Tarnocai, C., Schuur, E. A. G., and Allison, S. D.: Causes
of variation in soil carbon simulations from CMIP5 Earth system
models and comparison with observations, Biogeosciences, 10,
1717–1736, https://doi.org/10.5194/bg-10-1717-2013, 2013.
Tóth, G., Jones, A., and Montanarella, L.: LUCAS Topsoil sur-
vey: methodology, data and results, Report, Land Resource Man-
agement Unit – Soil Action, European Commission Joint Re-
search Centre Institute for Environment and Sustainability, avail-
able at: https://esdac.jrc.ec.europa.eu/ESDB_Archive/eusoils_
docs/other/EUR26102EN.pdf (last access: 21 May 2021), 2013.
UNEP: The benefits of soil carbon – managing soils for multi-
ple, economic, societal and environmental benefits, United Na-
tions Environmental Programme, Nairobi, Kenya, 19–33, https:
//tinyurl.com/2vu4nsf4 (last access: 21 May 2021), 2012.
van Bussel, L. G. J., Grassini, P., Van Wart, J., Wolf, J., Claessens,
L., Yang, H., Boogaard, H., de Groot, H., Saito, K., Cassman,
K. G., and van Ittersum, M. K.: From field to atlas: Upscaling
of location-specific yield gap estimates, Field Crop. Res., 177,
98–108, https://doi.org/10.1016/j.fcr.2015.03.005, 2015.
van der Esch, S., Brink, B. T., Stehfest, E., Bakkenes, M., Sewell,
A., Bouwman, A., Meijer, J., Westhoek, H., and van den Berg,
M.: Exploring future changes in land use and land condition
and the impacts on food, water, climate change and biodiver-
sity: Scenarios for the UNCCD Global Land Outlook, Report,
UNCCD, available at: https://tinyurl.com/yagvs9vu (last access:
21 May 2021), 2017.
van Ittersum, M. K., Cassman, K. G., Grassini, P., Wolf, J., Tit-
tonell, P., and Hochman, Z.: Yield gap analysis with local
to global relevance – A review, Field Crop. Res., 143, 4–17,
https://doi.org/10.1016/j.fcr.2012.09.009, 2013.
Vitharana, U. W. A., Mishra, U., and Mapa, R. B.: National soil or-
ganic carbon estimates can improve global estimates, Geoderma,
337,
55–64,
https://doi.org/10.1016/j.geoderma.2018.09.005,
2019.
Wan, Z.: MODIS land surface temperature products users’ guide,
Institute for Computational Earth System Science, Univer-
sity of California, Santa Barbara, California, USA, avail-
able at: https://lpdaac.usgs.gov/documents/118/MOD11_User_
Guide_V6.pdf, 2006.
Wilson, A. M. and Jetz, W.: Remotely sensed high-resolution
global
cloud
dynamics
for
predicting
ecosystem
and
biodiversity
distributions,
PLoS
Biol.,
14,
e1002415,
https://doi.org/10.1371/journal.pbio.1002415, 2016.
WOCAT: Where the land is greener: Case studies and analysis of
soil and water conservation initiatives worldwide, CTA, UNEP,
FAO and CDE, Berne, 2007.
Wright, M. N. and Ziegler, A.: ranger: A Fast Implementation
of Random Forests for High Dimensional Data in C++ and R,
J. Stat. Softw., 77, 1–17, https://doi.org/10.18637/jss.v077.i01,
2017.
Yigini, Y. and Panagos, P.: Assessment of soil organic car-
bon stocks under future climate and land cover changes
in
Europe,
Sci.
Total
Environ.,
557–558,
838–850,
https://doi.org/10.1016/j.scitotenv.2016.03.085, 2016.
Yoo, A. B., Jette, M. A., and Grondona, M.: Slurm: Simple linux
utility for resource management, in: Workshop on Job Schedul-
ing Strategies for Parallel Processing, 44–60, 2003.
Zanter, K.: Landsat 4–7 Surface Reflectance Code (LEDAPS) Prod-
uct Guide, available at: https://www.usgs.gov/land-resources/nli/
landsat/landsat-normalized-difference-vegetation-index (last ac-
cess: 21 May 2021), 2019.
SOIL, 7, 217–240, 2021
https://doi.org/10.5194/soil-7-217-2021