Back to report · Publisher record

Poggio et al. (2021). SoilGrids 2.0: producing soil information for the globe with quantified spatial uncertainty

Page 1

Figures and page graphics, page 1

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.

Page 2

Figures and page graphics, page 2

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

Page 3

Figures and page graphics, page 3

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

Page 4

Figures and page graphics, page 4

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

Page 5

Figures and page graphics, page 5

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

Page 6

Figures and page graphics, page 6

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

Page 7

Figures and page graphics, page 7

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

Page 8

Figures and page graphics, page 8

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

Page 9

Figures and page graphics, page 9

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

Page 10

Figures and page graphics, page 10

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

Page 11

Figures and page graphics, page 11

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

Page 12

Figures and page graphics, page 12

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

Page 13

Figures and page graphics, page 13

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

Page 14

Figures and page graphics, page 14

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

Page 15

Figures and page graphics, page 15

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

Page 16

Figures and page graphics, page 16

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

Page 17

Figures and page graphics, page 17

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

Page 18

Figures and page graphics, page 18

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

Page 19

Figures and page graphics, page 19

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

Page 20

Figures and page graphics, page 20

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

Page 21

Figures and page graphics, page 21

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

Page 22

Figures and page graphics, page 22

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

Page 23

Figures and page graphics, page 23

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

Page 24

Figures and page graphics, page 24

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


Document Outline