Hengl et al. (2017). SoilGrids250m: Global gridded soil information based on machine learning

5

]. Spatial information

about soil water parameters is likely to become increasingly critical in areas affected by climate
change [

6

]. Soils and soil information are also particularly relevant for the Sustainable Devel-

opment goal target 15.3 of achieving Land Degradation Neutrality (LDN), as specified by the
United Nations Convention to Combat Desertification (UNCCD;

http://www.unccd.int

), and

are one of the main areas of interest of the FAO’s Global Soil Partnership initiative [

7

]. Fol-

berth et al. [

8

] have recently discovered that accurate soil information might be the key to pre-

dicting either buffering or amplifying impacts of climate change on food production.

To reduce the gap between soil data demand and availability, ISRIC (International Soil

Reference Information Centre)—World Soil Information released a Global Soil Information
system called

“SoilGrids”

. The first version of SoilGrids (predictions at 1 km spatial resolu-

tion released in 2014), was, at the time, a ‘proof of concept’ demonstrating that global compi-
lations of soil profiles can be used in an automated framework to produce complete and
consistent spatial predictions of soil properties and classes [

9

]. Since the launch of the system

in 2014, several colleagues have recognized and reported some of the limitations of the first
version of the system. Mulder et al. [

10

] observed, using more detailed soil profile data and

maps, that SoilGrids likely overestimated all low values for organic carbon content in France.
Likewise, Griffiths et al. [

11

] reported underestimation of the pH in comparison to UK

national data. The overestimation of low values happened mainly as an effect of limited fit-
ting success (so that both high and low values are smoothed out). In addition, many of the
artifacts visible in the Harmonized World Soil Database (HWSD) [

12

], which was used as

one of the covariates to produce the first version of SoilGrids, e.g. country borders, were
propagated to SoilGrids1km. Some users have also expressed concerns that the first version
of SoilGrids did not provide predictions for arid and desert areas and hence can be consid-
ered an incomplete product [

13

ISRICWorldSoil/SoilGrids250m/

To address these criticisms and concerns, we have re-designed and re-implemented Soil-

Grids with a particular emphasis on addressing methodological limitations of SoilGrids1km.
Hence, our main objective was to build a more robust system with improved output data qual-
ity; especially considering spatial detail and attribute accuracy of spatial predictions. We imple-
mented the following six key improvements:

1. We replaced linear models with tree-based, non-linear machine learning models to account

for non-linear relationships—especially for modeling soil property–depth relationships—
but also to be able to better represent local soil–covariate relationships. Predictions are now
primarily data-driven. Much less time is spent on choosing models, which also reduces the
complexity of producing updates.

2. We replaced single prediction models with an ensemble framework i.e. we use at least two

methods for each soil variable to reduce overshooting effects.

SoilGrids250m: Global gridded soil information

PLOS ONE | DOI:10.1371/journal.pone.0169748

February 16, 2017

2 / 40

used to generate SoilGrids250m predictions is fully

documented via:

https://github.com/

.

Funding: ISRIC is a non-profit organization

primarily funded by the Dutch government. The

funders had no role in study design, data collection

and analysis, decision to publish, or preparation of

the manuscript. GILAB DOO provided support in

the form of salaries for author AB, but did not have

any additional role in the study design, data

collection and analysis, decision to publish, or

preparation of the manuscript. The specific roles of

this author are articulated in the ‘author

contributions’ section.

Competing interests: Aleksandar Blagoti

ć

is

employee and web-developer of GILAB DOO. There

are no patents, products in development or

marketed products to declare. This does not alter

our adherence to all the PLOS ONE policies on

sharing data and materials.

3. We extended the initial list of covariates to include a wider diversity of MODIS land prod-

ucts and to better represent factors of soil formation. The spatial resolution of covariates
was increased from 1 km to 250 m with the expectation that finer resolution will help
increase the prediction accuracy.

4. We re-implemented the global soil mask using state-of-the-art land cover products [

14

ftp://ftp.soilgrids.org/data/

The current soil mask now includes all previously excluded dryland and sand dune areas so
that most of the land mask (

>

95%) is represented.

5. The global compilation of soil profiles and samples used for model training was also

extended. We added extra points for the Russian Federation, Brazil, Mexico and the Arctic
circle; and re-visited data harmonization issues.

6. We created and inserted expert-based pseudo-points for a selection of parameters to mini-

mize extrapolation effects in undersampled geographic areas lacking field observations,
such as deserts, semi-deserts, glaciers and permafrost areas.

We present here the technical development and accuracy assessment of the updated Soil-

Grids system at 250 m resolution. In the following sections we describe the workflows used to
generate spatial predictions and report results of model fitting and accuracy assessment based
on 10–fold cross-validation. We conclude the article by suggesting some possible applications
of this new data set and identifying possible future improvements. SoilGrids250m map layers
are available for download via

www.SoilGrids.org

under the Open Database License (ODbL).

GeoTiffs can also be obtained from

.

Methods and materials

Target variables

SoilGrids provides predictions for the following list of standard soil properties and classes [

9

]:

• Soil organic carbon content in ‰ (g kg

−1

),

• Soil pH in H

2

O and KCl solution,

• Sand, silt and clay (weight %),

• Bulk density (kg m

−3

) of the fine earth fraction (

<

2 mm),

• Cation-exchange capacity (cmol + /kg) of the fine earth fraction,

• Coarse fragments (volumetric %),

• Depth to bedrock (cm) and occurrence of

R

horizon,

• World Reference Base (WRB) class—at present, we map 118 unique soil classes, e.g. Plinthic

Acrisols, Albic Arenosols, Haplic Cambisols (Chromic), Calcic Gleysols and similar [

15

http://sis.agr.gc.ca/cansis/taxa/

This is about four times as many classes as in the previous version of SoilGrids,

• United States Department of Agriculture (USDA) Soil Taxonomy suborders—i.e. 67 soil

classes [

16

].

We generated predictions at seven standard depths for all numeric soil properties (except

for depth to bedrock and soil organic carbon stock): 0 cm, 5 cm, 15 cm, 30 cm, 60 cm, 100 cm
and 200 cm, following the vertical discretisation as specified in the GlobalSoilMap specifica-
tions [

17

]. Averages over (standard) depth intervals, e.g. 0–5 cm or 0–30 cm, can be derived by

taking a weighted average of the predictions within the depth interval using numerical

SoilGrids250m: Global gridded soil information

PLOS ONE | DOI:10.1371/journal.pone.0169748

February 16, 2017

3 / 40

integration, such as the trapezoidal rule:

1

b

a

Z

b

a

f

ð

x

Þ

dx

1

ð

b

a

Þ

1

2

X

N

1

k

¼

1

x

k

þ

1

x

k

f

ð

x

k

Þ þ

f

ð

x

k

þ

1

Þ

ð

1

Þ

where

N

is the number of depths,

x

k

is the

k

-th depth and

f

(

x

k

) is the value of the target variable

(i.e., soil property) at depth

x

k

. For example, for the 0–30 cm depth interval, with soil pH values

at the first four standard depths equal to 4.5, 5.0, 5.3 and 5.0, the pH is estimated as

1

30

2

5

0

ð

Þ

4

:

5

þ

5

:

0

ð

Þ þ

15

5

ð

Þ

5

:

0

þ

5

:

3

ð

Þ þ

30

15

ð

Þ

5

:

3

þ

5

:

0

ð

Þ

½

=

30

0

:

5

¼

5

:

083

(

Fig 1

).

Based on predictions of soil organic carbon content, bulk density, and coarse fragments, we

also derived soil organic carbon stock (tha

−1

) for the six GlobalSoilMap standard depth inter-

vals following the standard approach [

].

shows an example of observed vs predicted

values and corresponding derived soil organic carbon stock for 0–1 m and 1–2 m depths.

Model fitting and spatial prediction of depth to bedrock is based also on water well drilling

data. Model fitting and spatial prediction of soil depth to bedrock variables is explained in
detail in Shangguan et al. [

19

].

We set the reference soil surface at the air/soil boundary, as per FAO [

20

], hence all soil

material is included. Some national soil survey teams (and also earlier versions of the FAO

Fig 1. Standard soil depths following the GlobalSoilMap.net specifications and example of numerical
integration following the trapezoidal rule.

doi:10.1371/journal.pone.0169748.g001

SoilGrids250m: Global gridded soil information

PLOS ONE | DOI:10.1371/journal.pone.0169748

February 16, 2017

4 / 40

standards) define 0 cm depth at the start of the mineral soil, i.e. just below the

O

or the

P

(peat)

horizon. Consider for example the following sample soil profile from Canada [

21

]:

hor

top

bottom

bd

orgcarb

LFH

-12

0

0.07

48.1

Ae

0

11

1.3

0.6

AB

11

25

1.53

0.4

Bt

25

44

1.62

0.4

which shows that the vertical coordinates of the organic layer of this soil site are negative (

LFH

indicates Litter—Fermentation—Humus);

orgcarb

indicates soil organic carbon,

bd

is the

bulk density and

top

and

bottom

are the upper and lower horizon depth in cm). Therefore,

to avoid vertical mismatches between different national systems, all systems that put the zero
level at the start of the mineral soil have been adjusted to a reference with the zero level at the
air/soil boundary. For the example soil profile from Canada this means that 12 cm was added
to all top and bottom values (in the example above, there is a significant discontinuity in values
in organic carbon that drops from 48.1% to 0.6% within 12 cm of depth).

Fig 2. Example of soil variable-depth curves: Original sampled soil profiles (black rectangles) vs predicted SoilGrids values at seven standard
depths (broken red line), and predicted soil organic carbon stock for depth intervals 0–100 and 100–200 cm. Locations of points from the USDA
National Cooperative Soil Survey Soil Characterization database: mineral soil

S1991CA055001

(-122.37˚W, 38.25˚N), and an organic soil profile

S2012CA067002

(-121.62˚W, 38.13˚N).

doi:10.1371/journal.pone.0169748.g002

SoilGrids250m: Global gridded soil information

PLOS ONE | DOI:10.1371/journal.pone.0169748

February 16, 2017

5 / 40

Input profile data

For model building, we used soil profile data from ca. 150,000 unique sites spread over all con-
tinents (

Fig 3

; see acknowledgments for a full list). These have been imported, cleaned and

merged into a single global compilation of soil points with unique column names and IDs.

Preparation of the global compilation of standardized soil training points took several

months of work. The translation and cleaning up of soil properties and soil classes took a large
amount of time. About 15–20% of the original soil profile data was only reported using a
national classification system, e.g. the Canadian and Brazilian classification systems. Since
some information is better than none, where possible we translated national classification sys-
tems to the two international (World Reference Base and USDA) classification systems. For
translation we used published correlation tables either reported in Krasilnikov et al. [

22

] or

reported on the agency websites; see e.g. correlation of Canadian Soil Taxonomy published
(

) and correlation of the Brazilian classification system (

http://

www.pedologiafacil.com.br/classificacao.php

). We also consulted numerous local soil classifi-

cation experts and requested their feedback and corrections in the (online) correlation tables
(distributed via Google spreadsheets). Some national classification systems, such as the Austra-
lian soil classification system, are simply too different from the USDA and WRB systems to
allow satisfactory correlation. These data were therefore not used. The full list of correlation
tables is available from ISRIC’s github account at

https://github.com/ISRICWorldSoil

.

Another time-consuming operation was merging laboratory measurements and field obser-

vations and their harmonization to a standard format. In some cases missing values in the
original tables had been coded as

"0"

values, which can have a serious influence on prediction

models; in other cases we implemented and applied functions to locate and correct typos and
other gross errors. Some variables, such as soil organic carbon, needed to be converted either
from soil organic matter (e.g. divide by 1.724) and/or by removing CaCO

3

(Calcium carbon-

ates) from total carbon. Nevertheless, the majority of soil variables from various national soil

Fig 3. Input profile data: World distribution of soil profiles used for model fitting (about 150,000 points shown on the map; see
acknowledgments for a complete list of data sets used). Yellow points indicate pseudo-observations. For the majority of points shown on this map,
laboratory data can be accessed from ISRIC’s World Soil Information Service (WoSIS) at

http://wfs.isric.org/geoserver/wosis/wfs

.

doi:10.1371/journal.pone.0169748.g003

SoilGrids250m: Global gridded soil information

PLOS ONE | DOI:10.1371/journal.pone.0169748

February 16, 2017

6 / 40

profile data bases appeared to be compatible and relatively easy to merge—soil scientists across
continents do measure similar things, but often express the results using different measure-
ment units, vocabularies and standards.

We imported all original tables

as-is

, next documented all conversion functions through

R

scripts (available via ISRIC’s github account), to accommodate reproducible research and facil-
itate that conversion functions may, in the future, be further modified and improved. The
majority of the points (excluding LUCAS points and other data sets with specific restricting
terms of use) and legends used for model building and for producing SoilGrids are also avail-
able for public use via ISRIC’s WoSIS Web Feature Service (

http://www.isric.org/data/wosis

)

and/or the ISRIC’s institutional github account.

Expert-based pseudo-observations

Even though the input training point data are extensive and cover most continents and cli-
matic zones, some large areas that have extreme climatic conditions and/or have very
restricted access, are significantly undersampled. This occurs largely in the following five types
of areas:

1. Semi-arid and arid lands, deserts and sand dunes,

2. Mountain tops, steep slopes of mountains and similar inaccessible areas,

3. Areas covered by ice and/or snow, i.e. glaciers,

4. Inaccessible tropical forest,

5. Areas governed by totalitarian and hostile regimes, with military conflicts or war.

It might seem obvious to soil surveyors that there is no soil organic carbon in the top 2 m of

the active sand dunes of the Sahara, but any model fitted without observations in the Sahara
could result in dubious extrapolation and questionable predictions. In addition, relationships
across transitional areas—from semi-arid zones to deserts—can be difficult to represent with-
out enough points at both edges of the feature space. Some sand dunes in the USA have been
actually sampled and analyzed in the laboratory. For example, Lei [

23

] has shown that sand

dunes in the Mojave desert have an average pH of 8.1, 98% sand and 0% organic carbon.
Again, although it might seem obvious that deserts consist mainly of sand, and that steep
slopes without vegetation are either very shallow or show bedrock at the surface, the model is
not aware of such expert knowledge and hence such features need to be ‘numerically repre-
sented’ in the calibration dataset. We therefore decided, instead of masking out all such areas
from soil mapping, to insert pseudo-observations and fill gaps in the feature space for the first
four of the five types of areas listed above, i.e. to add pseudo-observations to the training data-
set, which we then use for model building.

We used the following data sources to delineate sand dunes, bare rock and glaciers and pro-

duce their respective land masks:

• Sand dunes mask—To delineate the global distribution of sand dunes we used mean annual

long-term surface temperature generated from the MODIS LST data product (MOD11A2),
long-term MODIS Mid-Infrared (MIR) band (MCD43A4) and a slope map. After visual
inspection of the border of the Sahara desert, it was clear that sand dunes can be relatively
accurately delineated using MIR reflectance, mean daily annual temperature (

>

25˚C) and a

slope map (

<

25 rad).

• Bare rock mask—To delineate bare rock we also used the MODIS MIR band (MCD43A4)

and a slope map. Bare rock or dominantly rocky areas show high MIR surface reflectance

SoilGrids250m: Global gridded soil information

PLOS ONE | DOI:10.1371/journal.pone.0169748

February 16, 2017

7 / 40

and are associated with steep slopes (

>

32 rad). To the initial mask map estimated using

MODIS MIR band and slope map, we also added bare rock areas from more detailed maps
available for some countries, such as Iceland and northern Europe [

19

• Glaciers mask—To represent global distribution of glaciers we used the GLIMS Geospatial

Glacier Database [

24

For each of the three masks we then generated randomly 100–400 points based on the rela-

tive global extent and assigned soil properties and soil classes accordingly (e.g. in the case of
WRB’s Protic Arenosols for sand dunes, Lithic and Rendzic Leptosols for bare rock areas,
Cryosols for areas adjacent to glaciers; in the case of USDA’s Psamments for sand dunes, Orth-
ents for bare rock areas and Turbels for glaciers; for sand dunes we also inserted estimated val-
ues of 0 soil organic carbon, 98% sand and 0% coarse fragments). For model training for
predicting soil classes we also used pseudo-observations generated from the best available soil
polygon maps: for poorly accessible tropical forest areas, such as Indonesia, we used the Land
information system of Kalimantan [

25

], and for northern latitudes, i.e. to represent permafrost

soils, the Northern Circumpolar Soil Carbon Database was used [

26

].

When inserting pseudo-observations we tried to follow three simple rules of thumb to min-

imize any negative effects:

• keep the relative percentage of pseudo-points small i.e. try not to exceed 1–2% of the total

number of training points,

• only insert pseudo-points for which the actual ground value is known with high confidence,

e.g. sand content in sand dune areas,

• if polygon maps are used to insert pseudo-observations, we tried to use the most detailed soil

polygon maps and focus on polygons with very high thematic purity.

Soil covariates

As covariate layers for producing SoilGrids250m predictions we used an extensive stack of
covariates, which are primarily based on remote sensing data. These include (see e.g.

Fig 4

):

• DEM-derived surfaces—slope, profile curvature, Multiresolution Index of Valley Bottom

Flatness (VBF), deviation from Mean Value, valley depth, negative and positive Topographic
Openness and SAGA Wetness Index—all based on the global merge of SRTMGL3 DEM and
GMTED2010 [

27

]. All DEM derivatives were computed using

SAGA GIS

28

],

• Long-term averaged monthly mean and standard deviation of the MODIS Enhanced Vege-

tation Index (EVI). Derived using a stack of MOD13Q1 EVI images [

29

• Long-term averaged mean monthly surface reflectances for MODIS bands 4 (NIR) and 7

(MIR). Derived using a stack of MCD43A4 images [

30

• Long-term averaged monthly mean and standard deviation of the MODIS land surface tem-

perature (daytime and nighttime). Derived using a stack of MOD11A2 LST images [

31

],

• Long-term averaged mean monthly hours under snow cover based on a stack of MOD10A2

8-day snow occurrence images [

32

],

• Land cover classes (cultivated land, forests, grasslands, shrublands, wetlands, tundra, artifi-

cial surfaces and bareland cover) for the year 2010 based on the GlobCover30 product by the
National Geomatics Center of China [

14

]. Upscaled to 250 m resolution and expressed in

percent of pixel coverage,

SoilGrids250m: Global gridded soil information

PLOS ONE | DOI:10.1371/journal.pone.0169748

February 16, 2017

8 / 40

• Monthly precipitation images derived as the weighted average between the WorldClim

monthly precipitation [

33

] and GPCP Version 2.2 [

34

• Long-term averaged mean monthly hours under snow cover. Derived using a stack of

MOD10A2 8-day snow occurrence images,

Fig 4. Examples of covariates used to generate SoilGrids: TWI is the Topographic Wetness Index (values multiplied by 100), EVI is the MODIS
Enhanced Vegetation Index (values multiplied by 10,000), s.d. LST is the long-term standard deviation of MODIS Land Surface Temperatures
(values in Celsius degrees). Location: San Francisco bay area, California. Size of the bounding box is 300 by 300 km.

doi:10.1371/journal.pone.0169748.g004

SoilGrids250m: Global gridded soil information

PLOS ONE | DOI:10.1371/journal.pone.0169748

February 16, 2017

9 / 40

• Lithologic units (acid plutonics, acid volcanic, basic plutonics, basic volcanics, carbonate

sedimentary rocks, evaporite, ice and glaciers, intermediate plutonics, intermediate volca-
nics, metamorphics, mixed sedimentary rocks, pyroclastics, siliciclastic sedimentary rocks,
unconsolidated sediment) based on Global Lithological Map GLiM [

35

http://oas.gsfc.nasa.gov/floodmap/

• Landform classes (breaks/foothills, flat plains, high mountains/deep canyons, hills, low

hills, low mountains, smooth plains) based on the USGS’s Map of Global Ecological Land
Units [

36

].

• Global Water Table Depth in meters; after Fan et al. [

37

],

• Long-term averaged mean monthly MODIS Flood Water based on the NRT Global MODIS

Flood Mapping Flood Water product (

),

• Landsat-based estimated distribution of Mangroves; after Giri et al. [

38

• Average soil and sedimentary-deposit thickness in meters; after Pelletier et al. [

39

].

These covariates were selected to represent factors of soil formation according to Jenny

40

]: climate, relief, living organisms, water dynamics and parent material. Out of the five

main factors, water dynamics and living organisms (especially vegetation dynamics) are not
trivial to represent as these operate over long periods of time and often exhibit chaotic behav-
iour. Using reflectance bands such as the mid-infrared MODIS bands from a single day, would
have little use to soil mapping for areas with dynamic vegetation, i.e. with strong seasonal
changes in vegetation cover. To account for seasonal fluctuation and for inter-annual varia-
tions in surface reflectance, we instead used long-term temporal signatures of the soil surface
derived as monthly averages from long-term MODIS imagery (15 years of data). We assume
here that, for each location in the world, long-term average seasonal signatures of surface
reflectance or vegetation index provide a better indication of soil characteristics than only a
single snapshot of surface reflectance. Computing temporal signatures of the land surface
requires a considerable investment of time (comparable to the generation of climatic images vs
temporary weather maps), but it is possibly the only way to represent the cumulative influence
of living organisms on soil formation.

For processing the covariates we used a combination of Open Source GIS software, primar-

ily

SAGA GIS

28

],

R

packages

raster

[

41

],

sp

42

GSIF

and

GDAL

43

] for reprojecting,

mosaicking and merging tiles.

SAGA GIS

and

GDAL

were found to be highly suitable for pro-

cessing large data as parallelization of computing was relatively easy to implement.

We updated the 1 km global soil mask map using the most detailed 30 m resolution global

land cover map from 2010. This was combined with the global water mask [

44

] and the global

sea mask map based on the SRTM DEM [

45

] to produce one consistent global soil mask that

includes all land areas, expect for: (a) fresh water bodies such as lakes and rivers, and (b) per-
manent ice.

Spatial prediction framework

Spatial prediction, i.e. fitting of models and generation of maps, was fully implemented via the

R

environment for statistical computing. The process of generating SoilGrids predictions con-

sists of four main steps (see

Fig 5

):

• overlay points and covariates and prepare regression matrix,

• fit spatial prediction models,

SoilGrids250m: Global gridded soil information

PLOS ONE | DOI:10.1371/journal.pone.0169748

February 16, 2017

10 / 40

Fig 5. The (data-driven) statistical framework used for generating SoilGrids. SoilGrids are primarily based on publicly released soil profile
compilations, NASA’s MODIS and SRTM data products and Open Source software compiled with the ATLAS library:

R

(including contributed

packages), and Open Source Geospatial Foundation (OSGeo) supported software tools.

doi:10.1371/journal.pone.0169748.g005

SoilGrids250m: Global gridded soil information

PLOS ONE | DOI:10.1371/journal.pone.0169748

February 16, 2017

11 / 40

• apply spatial prediction models using tiled raster stacks (covariates),

• assess accuracy using cross-validation.

For practical purposes, we implemented these steps separately for each of the following

groups of soil variables:

• WRB soil groups and USDA soil suborders were modelled using ensemble models based

on

nnet::multinom

(which fits multinomial log-linear models via neural networks)

46

] and

ranger::ranger

(fits random forest) functions [

47

]. We mapped probabilities

of occurrence for each individual soil class (118 probability maps for WRB and 67 for
USDA),

• Soil properties (organic carbon, bulk density, CEC, pH, soil texture fractions and coarse

fragments) were modelled as 3D variables using an ensemble of

ranger::ranger

and

xgboost::xgboost

(fits Gradient Boosting Tree) [

48

]. Soil depth is used as a covariate,

so that the resulting models predict values of a target variable for any given depth, i.e. in
3D,

• Depth to bedrock was also modelled using

ranger::ranger

and

xgboost::xgboost

functions, but the output is a 2D map.

To optimize the model tuning parameters we consistently used the

caret::train

func-

tion [

49

], which is also suited for big data. The fine-tuning of the parameters is summarized in

the following three steps:

1. Randomly subset the regression matrix to e.g. 15,000 observations (usually 5–10% of the

total size),

2. Fit and validate a list of models for a combination of tuning parameters,

3. Select the optimal parameters (i.e. those that produce the lowest RMSE using repeated

cross-validation) and fit the final model using all observations.

Models for WRB and USDA classes are defined as:

R> TAXNWRB

*

DEMMRG5 + SLPMRG5 +

. . .

+ ASSDAC3

where

DEMMRG5 + SLPMRG5 +

. . .

+ ASSDAC3

are the covariate layers,

TAXNWRB

is the

observed taxonomic class in the WRB system (target variable). An example of a soil property
model is given by:

R> ORCDRC

*

Depth + DEMMRG5 +

. . .

+ ASSDAC3

where

DEMMRG5 +

. . .

+ ASSDAC3

are the covariate layers,

ORCDRC

is the value of organic

carbon observed (target variable), and

Depth

is the sampling / observation depth.

For each variable we fitted a separate model and merged predictions from at least two mod-

els to minimize overshooting effects [

50

]. The merging of predictions is done by using the

average model accuracy estimated during the fine-tuning of model parameters, i.e. as a
weighted average [

50

]:

f

ð

x

Þ ¼

P

M
k

¼

1

w

k

f

k

ð

x

Þ

P

M
k

¼

1

w

k

;

w

k

¼

1

s

2

k

;

CV

ð

2

Þ

where

f

ð

x

Þ

is the final ensemble prediction,

M

is the number of models,

w

k

is the model weight

and

s

2

k

;

CV

is the model squared prediction error obtained using cross-validation. In practice,

SoilGrids250m: Global gridded soil information

PLOS ONE | DOI:10.1371/journal.pone.0169748

February 16, 2017

12 / 40

both

ranger::ranger

and

xgboost::xgboost

report about the same error in most

cases, hence the final prediction is often close to the unweighted average.

We also applied post-processing, mainly to remove artifacts: in the case of soil classes, we

filter out all classes theoretically impossible to occur in a given area, such as Gypsisols in arctic
climatic zones, using a simple soil-climate matrix (documented on the project github). For tex-
ture fractions we also applied a standardization function to ensure that all predictions are
between 0 and 100, and that the fractions sum up to 100%, e.g.:

Sand

c

½

%

¼

Sand

ð

Sand

þ

Silt

þ

Clay

Þ

100

ð

3

Þ

where Sand

c

is the corrected sand content.

SoilGrids can be considered as a Big Data project, especially in terms of data volumes and

variety. The total size of all input and output data used to generate SoilGrids exceeds 30 TiB, so
that a first step in preparing SoilGrids250m was to obtain a Synology 12-Bay NAS storage
server with 60 TiB space. Handling such a large data set presented major challenges consider-
ing computational complexity and network bandwidth limitations. To optimize computing
performance, especially spatial overlay, model fitting, predictions and export of predictions,
we used exclusively parallelized versions of functions. For prediction, parallelization is already
implemented internally via the

ranger

or

xgboost

software; for other processes we primar-

ily used the

snowfall

package [

51

All processing was implemented on a single dedicated high performance server with 256

GiB RAM, 8 TiB hard disk space, 48 cores (Intel Xeon 2xE5-2690v3 24c/48t 2.6–3.5 GHz) and
running on Ubuntu 15.10 (Willy Werewolf) OS and R-cran 3.2.3 using ATLAS (Automatically
Tuned Linear Algebra Software) 3.11.38 library. Even after parallelization, producing predic-
tions for all soil variables and all depths took 10+ days of continuous computing, i.e. about 12
thousand CPU hours (about 90% of the computing time is invested in generating predictions).
Because the current system is fully scalable, the next update of SoilGrids could be completed in
shorter time frames, e.g. by boosting the number of computer cores, although this might also
greatly increase the production costs.

The tiling system

For tiling, we used the Equi7 Grid system [

52

] which splits the global land mass into seven sep-

arate planar grids (Europe and Asia are split into two land masses with some small overlap).
The Equi7 Grid system was selected for several practical reasons [

52

]:

1. The projections of the Equi7 Grid are equidistant and hence suitable for various geographic

analyses, especially for derivation of buffer distances and for hydrological DEM modeling,
i.e. to derive all DEM-based soil covariates,

2. Areal and shape distortions stemming from the Equi7 Grid projection are relatively small,

yielding a small grid oversampling factor,

3. The Equi7 Grid systems ensures an efficient raster data storage while suppressing inaccura-

cies during spatial transformation. Especially for high-resolution global data, these are
important features.

The global soil mask at 250 m resolution contains about 1.6 billion pixels (Africa: 330 mil-

lion, Europe: 110 million, North America: 230 million, South America: 210 million, Antartica:
0.05 million, Oceania: 140 million, Asia: 360 million). We provide the final outputs in both the
Equi7 Grid system and in geographical WGS84 coordinates. Final global mosaics in the

SoilGrids250m: Global gridded soil information

PLOS ONE | DOI:10.1371/journal.pone.0169748

February 16, 2017

13 / 40

WGS84 system were produced by reprojecting all pixels using

GDAL

warp and translate func-

tions [

43

]. The ground resolution of 250 m corresponds to a geographical resolution of 1/480

decimal degrees. An image representing the whole world at this resolution comprises 172

k

col-

umns and 72k rows.

Final predictions are available both as mosaics and as 1˚ tiles (16,360 tiles to represent the

world land mask); tiles are considered more suitable for users interested in regional and
national data, and mosaics (at resolutions of 20 km, 1 km and 250 m) are deemed suitable for
global modellers.

Accuracy assessment

For accuracy assessment of both numeric and categorical variables we used 10–fold repeated
cross-validation. Each model is re-fitted 10 times using 90% of the data and predictions
derived from the fitted models are compared with observations of the remaining 10%. For
each of the 14 numeric soil properties we derived the coefficient of determination (

R

2

—the

amount of variation explained by the model), mean error (ME) and root mean squared error
(RMSE). The amount of variation explained by the model is derived as:

R

2

¼

1

SSE

SST

100%

ð

4

Þ

where

SSE

is the sum of squared errors at cross-validation points and

SST

is the total sum of

squares. A coefficient of determination close to 1 indicates a perfect model, i.e. 100% of vari-
ation has been explained by the model. Numeric variables with skew distributions were log-
transformed prior to modeling and hence for these variables we report the amount of varia-
tion explained by the model after log-transformation. Also for the cross-validation correla-
tion plots we used either log or linear scale depending on whether log-transformation was
applied.

For predictions of soil WRB and USDA classes we calculated the map purity (0–100%) for

the dominant soil class at cross-validation points and weighted kappa metrics [

53

] as imple-

mented in the

psych

package. For the predicted probabilities of soil class occurrences (0–1

probability values) we also derived the area under the receiver operating characteristic curve
(AUC) and the True Positive Rate (TPR) statistic as implemented in the

ROCR

package [

54

,

55

]. Values of TPR range from 0 to 1. Values of AUC close to 1 show high prediction perfor-

mance, while values around 0.5 and below are considered poor.

For soil WRB and USDA classes we also generated global maps of the scaled Shannon

Entropy Index using the per-class probability maps [

56

57

]:

H

s

ð

x

Þ ¼

X

K

k

¼

1

p

k

ð

x

Þ

log

K

ð

p

k

ð

x

ÞÞ

ð

5

Þ

where

K

is the number of possible classes, log

K

is the logarithm to base

K

and

p

k

is probability

of class

k

. The scaled Shannon Entropy Index (

H

s

) is in the range from 0–1, where 0 indicates

no ambiguity (one of the

p

k

equals one and all others are zero) and 1 indicates maximum

confusion (all

p

k

equal

1

K

) [

58

]. Note that the scaled Shannon Entropy Index should not be

confused with classification accuracy assessment:

H

s

is an internal accuracy measure derived

from the model and not based on comparison of predictions with (cross-)validation data,
such as the purity and kappa metrics. For Shannon index of 0 at some location accuracy
could still be completely wrong because the soil class at that location could actually be a dif-
ferent one.

SoilGrids250m: Global gridded soil information

PLOS ONE | DOI:10.1371/journal.pone.0169748

February 16, 2017

14 / 40

Results

Model fitting

Summary results of model fitting are given in Figs

and

and Tables

and

. The

ranger

package reports model fitting success via the R-square based on Out-of-bag (OOB) samples,
i.e. the amount of variation explained by the model, which ranged from a low of 0.59 for coarse
fragments to a high of 0.85 for soil pH. R-square estimated using

xgboost

(derived using

repeated cross-validation) was lower, ranging from 0.37 for coarse fragments to 0.60 for soil

Fig 6. Fitted variable importance plots for target variables. Generated as an average of predictions using the

ranger

and

xgboost

packages (for soil

types results are based on the

ranger

model only).

DEPTH.f

is depth from soil surface,

T

MOD3

and

N

MOD3

are mean monthly temperatures daytime

and nighttime (red color),

TWI

,

DEM

,

VBF

and

VDP

are DEM-parameters (bisque color),

M

MOD4

are mean monthly MODIS NIR band reflectances (cyan

color),

P

MRG3

are mean monthly precipitation (blue color),

E

MOD5

are mean monthly EVI derivatives (dark green color),

VW

MOD1

are monthly MODIS

Precipitable Water Vapor images (orange color),

C

GLC5

are land cover classes (light green color), and

ASSDAC3

is the average soil and sedimentary-

deposit thickness (brown color).

doi:10.1371/journal.pone.0169748.g006

SoilGrids250m: Global gridded soil information

PLOS ONE | DOI:10.1371/journal.pone.0169748

February 16, 2017

15 / 40

pH. On average, the two packages report R-square values between 0.4–0.8 with an overall aver-
age of 0.60. This number corresponds closely to our results produced using 10–fold cross vali-
dation with repeated fitting. Comparing these new results to average R-square values of 0.38
for the original SoilGrids1km predictions reveals a significant improvement of close to + 50%.

The

train

function of the package

caret

usually picked a relatively high

Mtry

parame-

ter (number of variables randomly sampled as candidates at each split) as optimal for soil
properties: the optimized values ranged from 18 for coarse fragments to 22 for all other soil
properties. Higher

Mtry

is recommend for cases where the number of covariates is large and

Fig 7. Examples of relationships for target variables and the most important covariates: (top row) bulk density in kg m

−

3

, (middle row) soil pH,

and (bottom row) soil organic carbon in permilles (on log scale). Plots show target variables and the top three most important covariates as reported
by the random forest model.

DEPTH.f

is the observed depth from soil surface,

T09MOD3

is mean monthly temperature for September,

TMDMOD3

is mean

annual temperature,

PRSMRG3

is total annual precipitation,

M04MOD4

is mean monthly MODIS NIR band reflectance for April,

P07MRG3

is mean monthly

precipitation for July,

T01MOD3

is mean monthly temperature for January, and

T02MOD3

is mean monthly temperature for February.

doi:10.1371/journal.pone.0169748.g007

SoilGrids250m: Global gridded soil information

PLOS ONE | DOI:10.1371/journal.pone.0169748

February 16, 2017

16 / 40

multiple variables influence the target variables with equal importance [

59

]. For the Gradient

Boosting Tree method,

train

always selected the same combination of tuning parameters for

all soil properties:

nrounds = 100

,

max_depth = 3

,

eta = 0.4

,

gamma = 0

,

colsam-

ple_bytree = 0.8

and

min_child_weight = 1

. This may be because we limited the

combinations of tuning parameters to 10 to speed up processing speed. Higher values for

xgboost

tuning parameters are indicative of higher-level complexity of the model: many

Table 1. SoilGrids average prediction error for key soil properties based on 10–fold cross-validation. N = “Number of samples used for training”, ME =
“Mean Error”, MAE = “Mean Absolute Error”, RMSE = “Root Mean Squared Error” and R-square = “Coefficient of determination” (amount of variation explained
by the model). For variables with a skew distribution, such as organic carbon, coarse fragments and CEC, the accuracy statistics are also provided on log-
scale

.

Variable name

N

Min

Max

ME

MAE

RMSE

R-square

RMSE

R-square

Soil organic carbon
(gravimetric)

605,054

0

520

-0.292

10.2

32.8

63.5%

0.715

68.8%

pH index
(H

2

O solution)

604,019

2.1

11.0

-0.002

0.4

0.5

83.4%

Sand content
(gravimetric)

616,762

1%

94%

-0.037

9.0

13.1

78.6%

Silt content
(gravimetric)

613,750

2%

74%

0.023

6.7

9.8

79.4%

Clay content
(gravimetric)

625,159

2%

68%

-0.102

6.6

9.5

72.6%

Coarse fragments
(volumetric)

303,139

0%

89%

-0.104

5.5

10.9

55.9%

1.185

64.3%

Bulk density
(fine earth fraction)

140,596

250

2870

-1.574

108.3

164.7

75.8%

Cation-exchange capacity
(fine earth fraction)

393,585

0

234

-0.071

5.5

10.3

64.5%

0.483

67.0%

Depth to bedrock
(in cm)

1,580,798

0

125,000

-29

678

835

54.0%

1.12

42.8%

doi:10.1371/journal.pone.0169748.t001

Table 2. Mapping performance of SoilGrids250m compared to summary results for SoilGrids1km [

9

]. Amount of variation explained by models (

Eq 4

),

i.e. prediction accuracy for soil types was determined using 10–fold cross-validation. GSIF = “Global Soil Information Facilities”.

Variable name

Type

Units

GSIF
code

Amount of var. explained

(SoilGrids1km)

Amount of var. explained

(SoilGrids250m)

Relative

improvement

Soil organic carbon

3D

g kg

−

1

ORCDRC

22.9%

68.8%

200%

pH index
(H

2

O solution)

3D

10

−

1

PHIHOX

50.5%

83.4%

65%

Sand content
(gravimetric)

3D

kg kg

−

1

SNDPPT

23.5%

78.6%

234%

Silt content
(gravimetric)

3D

kg kg

−

1

SLTPPT

34.9%

79.4%

127%

Clay content
(gravimetric)

3D

kg kg

−

1

CLYPPT

24.4%

72.6%

198%

Coarse fragments
(volumetric)

3D

cm

3

cm

−

3

CRFVOL

-

64.3%

-

Bulk density
(fine earth fraction)

3D

kg m

−

3

BLD

31.8%

75.8%

138%

Cation-exchange
capacity
(fine earth fraction)

3D

cmol +

/kg

CEC

29.4%

67.0%

128%

Depth to bedrock

2D

cm

BDT

-

42.8%

-

doi:10.1371/journal.pone.0169748.t002

SoilGrids250m: Global gridded soil information

PLOS ONE | DOI:10.1371/journal.pone.0169748

February 16, 2017

17 / 40

relationships between soil properties and covariates are non-linear and a greater number of
splits is possibly required to represent this complexity.

Fig 6

shows the top 15 soil covariates for each target variable. This indicates that, for exam-

ple, spatial pattern of soil pH is primarily influenced by precipitation and surface reflectance
(MODIS Medium-Infrared band 6 for months April and May especially). Also, for most vari-
ables depth emerges as the most important covariate, especially for soil organic carbon, bulk
density and coarse fragments. For soil types and soil textures, DEM-parameters, i.e. soil form-
ing factors of relief, especially flow-based DEM-indices, emerge as second-most dominant
covariates. These results largely correspond with conventional soil survey knowledge (survey-
ors have been using relief as a key guideline to delineate soil bodies for decades), but it is
encouraging to have these findings supported by statistical modeling of real data on a global
scale.

Although lithology is not in the list of top 15 most important predictors, spatial patterns of

lithologic classes can often be distinctly recognized in the output predictions. This is especially
true for soil texture fractions and coarse fragments. In general, for predicting soil chemical
properties, climatic variables (especially precipitation) and surface reflectance seem to be the
most important, while for soil classes and soil physical properties it is a combination of relief,
vegetation dynamics and parent material.

Fig 7

shows some individual relationships between target variables and several of the most

important covariates. For soil pH we observe that the relationship with total annual rainfall is
close to linear; for soil organic carbon and depth the relationship is linear on a log-log scale.
Many such individual correlations can also be interpreted and understood in terms of pedo-
logic knowledge. For example, higher MIR reflectance may be associated with high concentra-
tion of salts in soil and hence higher pH; higher rainfall and cooler climates often result in
higher organic carbon content because the speed of organic matter accumulation is higher
than the speed of decomposition. For the majority of soil variables, however, relationships are
not clearly linear and often many soil covariates are equally important.

We have also investigated possibilities for using kriging of residuals to improve predictions

of soil properties. Because the majority of spatial variation has been explained by covariates
and machine learning models, it appears that no significant spatial autocorrelation structure
can be observed for residuals (i.e. almost all variograms show pure nugget effect structure) at
distances

<

300 km for almost all continents and all variables. Although locally, where the

point density is high, kriging of residuals could still be beneficial for mapping of CEC and
depth to bedrock, overall kriging of residuals for global land mass does not seem to be neces-
sary nor is it practical to implement for billions of pixels: it would only marginally improve the
accuracy of predictions at high computing costs.

Accuracy assessment

Table 1

shows summary results of cross-validation for soil properties (global assessment). In

all cases there is no large overestimation of values, although for organic carbon and CEC the
models seem to somewhat under-estimate the overall mean. For log-transformed variables we
applied the accuracy assessment in the log-transformed space which yields asymmetric predic-
tion intervals after back-tranformation. For example, predictions for organic carbon are
±0.715 in log-space, which means that the 90% probability prediction interval for a case where
the soil organic carbon prediction equals 20‰ (2%) is 6–65‰; for a case where the soil organic
carbon prediction equals 150‰ it is 46–485‰. Prediction intervals are hence still fairly wide,
which might make SoilGrids of limited usability for detailed spatial modeling e.g. at farm level.
Note also that because there is significant spatial clustering of the training points, it is possible

SoilGrids250m: Global gridded soil information

PLOS ONE | DOI:10.1371/journal.pone.0169748

February 16, 2017

18 / 40

that the validation results might be somewhat more optimistic than if we had validated predic-
tions by using points collected following some (objective) probability sampling, as described in
Brus et al. [

60

]. On the other hand, the cross-validation results do not show any serious sys-

tematic over- or underestimation (ME close to zero), which is also visible from the correlation
plots (

Fig 8

).

Table 2

shows results for SoilGrids250m in comparison with the previous system at 1 km res-

olution. Improvements in average RMSE are between 30–80% and can largely be attributed to
the use of machine learning algorithms in place of multiple linear regression, but also to invest-
ments in preparing finer resolution covariates and additional and improved soil profile data.

Fig 8. Correlation (density) plots produced as a result of 10–fold cross-validation. See also

Table 1

for more details.

doi:10.1371/journal.pone.0169748.g008

SoilGrids250m: Global gridded soil information

PLOS ONE | DOI:10.1371/journal.pone.0169748

February 16, 2017

19 / 40

The most challenging variables to model with this set of covariates are coarse fragments

and depth to bedrock, although in no case is the R-square

<

50%. Nevertheless, the RMSE is

still relatively high in comparison to many local soil mapping projects. Users should thus be
aware that the uncertainty levels are still relatively high. There are also still problems with over-
estimation of low values, clearly visible for example in the case of soil organic carbon content.
Overall, predictions for most properties are unbiased, i.e. most predictions are fairly symmet-
ric around the 1:1 line (

Fig 8

).

For soil classes, out-of-bag average prediction accuracy, reported by the packages, was

between 20–28% for the WRB classification system and between 34–48% for the USDA sys-
tem. The 10–fold cross-validation results showed that the weighted kappa for WRB classes is
42%, with map purity 28%; for USDA classes the kappa is 57%, while the map purity is 48%.
Although WRB classes seem to be somewhat more challenging to model than USDA subor-
ders, this comparison should be considered within the context of: (a) the number of classes,
and (b) similarity between classes. The WRB classification contains about two times more
classes than USDA suborders, and many WRB classes with highest confusion fall in taxo-
nomically similar groups. Further evaluation of classification accuracy has shown that, at the
level of WRB soil groups, map purity jumps to 60%, i.e. it becomes comparable to the USDA
system. Remaining WRB soil groups with map purity

<

50% are Planosols, Phaeozems and

Ferrasols.

A more detailed assessment of prediction accuracy derived using the

ROCR

package, i.e. per

each individual class, shows that the average TPR is about 0.93 for USDA soil suborders
(

Table 3

), and about 0.90 for WRB classes (

Table 4

). Also maps of the scaled Shannon Entropy

index (

Fig 9

) indicate that produced soil class maps for USDA soil classification system are less

uncertain than for the WRB system: WRB classification is critically uncertain for Australia and
India, parts of Africa and highlands of Latin America. Maps of uncertainty closely reflect
extrapolation areas and could be potentially very useful for planning new soil surveys aimed at
mapping soil types. For example,

Fig 10

shows that the highest confusion (lowest prediction

accuracy) is systematically connected with distribution of river valleys, urban areas and hill-
slopes.

In summary, the cross-validation results for predicting class probabilities indicate relatively

high correspondence between prediction probabilities and observed soil types, which is also
confirmed visually by overlaying observed classes and prediction probabilities. Nevertheless, it
appears from Tables

3

and

4

that for some classes, such as Cambisols, Luvisols, Fluvisols and

Planosols in the WRB system, and Aquepts, Fluvents and Aquents in the USDA Soil Taxo-
mony system, the confusion of predictions with other classes is still relatively high.

Discussion

In the following sections we address some remaining discussion points and suggest ways to
improve SoilGrids and embark on new research directions. Although we have reached current
effective limits imposed by software capabilities and availability of remote sensing data sources,
the accuracy of SoilGrids could still be improved. Globally, by adding more covariates based
on the most recent remote sensing data (see

Fig 11

), and locally, by combining global predic-

tions with local prediction models. Global models could be further improved especially by
revising (even re-designing) each of the three main components of the system:

• Soil training data,

• Statistical / Machine Learning models, and

• Covariate layers.

SoilGrids250m: Global gridded soil information

PLOS ONE | DOI:10.1371/journal.pone.0169748

February 16, 2017

20 / 40

Table 3. Classification accuracy for predicted USDA class probabilities based on 10–fold cross-validation, ordered according to number of occur-
rences. ME = “Mean Error”, TPR = “True Positive Rate”, AUC = “Area Under Curve”, N = “Number of occurrences”, USDA = “United States Department of
Agriculture” soil classification system. The 1st and 2nd most probable classes are taken from the confusion matrix.

Name

ME (%)

TPR

AUC

N

1st class

2nd class

Udalfs

0.0

0.88

0.93

6326

Udalfs

Udults

0.0

0.91

0.95

4997

Udults

Udalfs

Udolls

0.1

0.91

0.93

3901

Udolls

Udalfs

Ochrepts

0.1

0.89

0.91

2720

Ochrepts

Udalfs

Aqualfs

0.0

0.89

0.91

2594

Aqualfs

Udalfs

Aquolls

0.1

0.89

0.90

2450

Udolls

Aquolls

Udox

0.0

0.93

0.95

2229

Ustox

Udox

Ustolls

-0.2

0.95

0.97

2042

Ustolls

Borolls

0.1

0.97

0.98

2029

Borolls

Albolls

Ustox

0.1

0.93

0.95

2024

Ustox

Udox

Orthents

0.1

0.88

0.89

1911

Orthents

Udults

Aquepts

0.0

0.87

0.88

1734

Aquolls

Aquepts

Psamments

0.1

0.90

0.92

1725

Psamments

Udults

Fluvents

0.1

0.84

0.85

1579

Udults

Udalfs

Orthods

0.1

0.97

0.98

1538

Orthods

Ochrepts

Udepts

0.1

0.90

0.91

1429

Udepts

Udults

Aquents

-0.1

0.84

0.85

1342

Aquepts

Udalfs

Ustalfs

-0.1

0.95

0.96

1332

Ustalfs

Ustolls

Xerolls

0.0

0.97

0.98

1319

Xerolls

Xeralfs

Argids

-0.1

0.98

0.99

907

Argids

Xerolls

Turbels

0.1

0.99

1.00

787

Turbels

Orthels

0.0

0.97

0.98

648

Ochrepts

Orthels

Xeralfs

-0.3

0.97

0.98

615

Xeralfs

Xerolls

Usterts

-0.2

0.97

0.98

590

Usterts

Ustolls

Albolls

-0.2

0.92

0.93

589

Borolls

Aquolls

Xerepts

-0.3

0.99

588

Xerepts

Xeralfs

Arents

-0.2

0.99

554

Arents

Ustox

Aquults

-0.2

0.94

380

Udults

Aquults

Cambids

-0.1

0.99

362

Cambids

Argids

Humults

-0.1

0.98

348

Humults

Udults

Hemists

-0.2

0.93

347

Ochrepts

Hemists

Torrox

-0.3

0.98

0.99

334

Ustox

Torrox

Saprists

-0.4

0.93

319

Saprists

Udalfs

Histels

-0.3

0.99

302

Histels

Turbels

Aquods

0.0

0.94

301

Orthods

Udults

Calcids

-0.1

0.98

0.99

301

Argids

Calcids

Ustults

0.0

0.99

286

Ustults

Ustalfs

Fibrists

0.0

0.96

0.97

250

Fibrists

Udults

Udands

0.0

0.99

234

Udands

Udox

Xerands

0.0

0.99

231

Xerands

Xerolls

Aquerts

-0.2

0.95

226

Aqualfs

Udalfs

Xererts

0.0

0.96

184

Xererts

Xerolls

Uderts

0.0

0.95

0.96

177

Udults

Uderts

Ustepts

0.0

0.97

175

Ustolls

Ustepts

Cryands

0.0

0.99

161

Cryands

Ochrepts

(Continued )

SoilGrids250m: Global gridded soil information

PLOS ONE | DOI:10.1371/journal.pone.0169748

February 16, 2017

21 / 40

Increasing and improving the quality and quantity of the training data

The most fruitful avenue for improving the current predictions is likely in improving the qual-
ity and quantity of soil profile data. ISRIC has invested decades in obtaining, digitizing, clean-
ing up and standardizing soil profile data. A large portion of these data (about 80,000 unique
points) is publicly available via ISRIC’s Web Feature Service WoSIS (

http://wfs.isric.org/

geoserver/wosis/wfs

) [

61

]; remaining soil profile data sets not publicly available via ISRIC’s

WoSIS WFS can be obtained by contacting the corresponding original data providers as listed
in the Acknowledgment section. This collection of soil profile data is of similar scope and util-
ity when compared to other international data initiatives in meteorology (e.g. Global Historical
Climatology Network) and biodiversity (

http://gbif.org

).

Although the training data shown in

Fig 3

appear to be quite dense, there are still large gaps

in terms of representation of the feature space. Tropics, wetlands, semi-arid to hyper-arid
areas and mountains are still largely under-represented. There are undoubtedly millions of soil
field observations in the world unused for global soil mapping activities that could be collated
and used to improve predictions. FAO’s Global Soil Partnership (

http://www.fao.org/

globalsoilpartnership/

) has set as one of its main objectives the preparation of an international

compilation of reference soil profiles to help catalyze using soil data for decision making.
Hence, there are already some initiatives in this direction.

Harmonization of soil laboratory data and soil descriptive variables is another area that will

need to be improved. For example, we had to standardize soil depths for several databases by
re-aligning 0 depth to soil surface. Some soil databases only contain information about the
mineral soil and put the zero level at the start of the mineral soil. But such soils might have an

Table 3. (Continued)

Name

ME (%)

TPR

AUC

N

1st class

2nd class

Cryepts

0.0

0.98

150

Ochrepts

Cryepts

Humods

0.0

0.92

149

Orthents

Orthods

Cryods

-0.1

0.99

1.00

133

Orthods

Cryods

Torrerts

-0.1

0.98

106

Ustolls

Torrerts

Cryolls

-0.2

0.99

79

Borolls

Cryolls

Gelods

-0.7

1.00

78

Turbels

Gelods

Gypsids

-0.1

0.99

70

Argids

Gypsids

Vitrands

-0.3

0.98

62

Vitrands

Ochrepts

Torrands

-0.3

0.99

60

Xerolls

Torrands

Durids

-0.3

0.99

59

Argids

Xerolls

Xerults

-0.3

0.97

53

Xeralfs

Humults

Rendolls

-0.5

0.93

41

Udalfs

Ochrepts

Salids

-0.6

0.94

37

Argids

Fluvents

Cryalfs

-0.7

1.00

32

Ochrepts

Borolls

Folists

-0.5

0.99

30

Orthods

Cryods

Gelands

-1.0

0.97

26

Gelods

Turbels

Perox

-0.7

0.99

21

Udults

Perox

Aquands

-0.7

0.98

19

Xerolls

Udands

Ustands

-0.9

1.00

17

Ustalfs

Orthents

Aquox

-1.0

0.98

16

Udults

Udox

Cryids

0.99

8

Argids

Borolls

Gelepts

0.83

6

Ochrepts

Turbels

doi:10.1371/journal.pone.0169748.t003

SoilGrids250m: Global gridded soil information

PLOS ONE | DOI:10.1371/journal.pone.0169748

February 16, 2017

22 / 40

Table 4. Classification accuracy for predicted WRB class probabilities based on 10–fold cross-validation, ordered according to number of occur-
rences. ME = “Mean Error”, TPR = “True Positive Rate”, AUC = “Area Under Curve”, N = “Number of occurrences”, WRB = “World Reference Base” soil clas-
sification system. The 1st and 2nd most probable classes are taken from the confusion matrix.

Name

ME (%)

TPR

AUC

N

1st class

2nd class

Haplic Cambisols

0.1

0.78

0.81

5619

Haplic Cambisols

Haplic Cambisols (Dystric)

Haplic Luvisols

0.1

0.86

0.88

2975

Haplic Luvisols

Haplic Cambisols

Haplic Acrisols

0.0

0.93

0.94

2607

Haplic Acrisols

Haplic Ferralsols

-0.1

0.94

0.95

1887

Haplic Ferralsols

Haplic Acrisols

Haplic Fluvisols

-0.1

0.88

0.89

1776

Haplic Fluvisols (Calcaric)

Haplic Fluvisols (Eutric)

Haplic Calcisols

0.0

0.93

0.94

1745

Haplic Calcisols

Calcaric Regosols

Haplic Kastanozems

0.0

0.96

0.97

1718

Haplic Kastanozems

Haplic Luvisols

Gleyic Luvisols

0.0

0.92

0.93

1686

Albic Luvisols

Gleyic Luvisols

Aric Regosols

-0.2

0.91

0.92

1488

Calcaric Regosols

Haplic Leptosols

Haplic Chernozems

0.1

0.96

0.97

1394

Haplic Chernozems

Haplic Kastanozems

Albic Luvisols

-0.1

0.94

0.95

1389

Gleyic Luvisols

Haplic Luvisols

Calcaric Regosols

-0.2

0.92

0.93

1379

Aric Regosols

Haplic Leptosols

Haplic Podzols

0.1

0.96

0.97

1359

Haplic Podzols

Haplic Cambisols

Haplic Cambisols (Dystric)

0.0

0.90

0.91

1334

Haplic Cambisols

Haplic Podzols

Haplic Cambisols (Calcaric)

-0.1

0.92

1173

Haplic Cambisols

Haplic Calcisols

Haplic Phaeozems

-0.1

0.90

0.91

1114

Haplic Phaeozems

Haplic Chernozems

Haplic Lixisols

0.0

0.92

0.93

1094

Haplic Lixisols (Chromic)

Haplic Lixisols

Haplic Leptosols

0.0

0.91

0.92

1092

Haplic Leptosols

Haplic Leptosols (Eutric)

Haplic Gleysols

0.0

0.88

0.89

1054

Haplic Gleysols (Eutric)

Haplic Cambisols

Haplic Vertisols

0.1

0.93

1040

Haplic Vertisols (Eutric)

Haplic Vertisols

Haplic Arenosols

0.1

0.91

0.92

935

Haplic Arenosols

Haplic Cambisols

Ferralic Arenosols

0.1

0.96

0.97

920

Ferralic Arenosols

Haplic Ferralsols

Haplic Cryosols

0.0

0.99

1.00

884

Haplic Cryosols

Haplic Cambisols

Haplic Cambisols (Eutric)

0.1

0.84

0.85

857

Haplic Cambisols

Haplic Luvisols

Haplic Alisols

0.1

0.94

0.95

827

Haplic Acrisols

Haplic Cambisols

Luvic Phaeozems

-0.4

0.90

0.91

741

Luvic Phaeozems

Haplic Luvisols

Rendzic Leptosols

0.0

0.94

0.95

695

Haplic Cambisols

Rendzic Leptosols

Haplic Fluvisols (Calcaric)

0.0

0.94

692

Haplic Fluvisols

Haplic Calcisols

Petric Calcisols

0.1

0.97

0.98

679

Petric Calcisols

Haplic Calcisols

Haplic Regosols (Eutric)

0.1

0.86

0.87

677

Haplic Cambisols

Haplic Luvisols

Lithic Leptosols

0.0

0.93

655

Haplic Ferralsols

Haplic Acrisols

Umbric Gleysols

0.0

0.91

0.92

621

Mollic Gleysols

Calcic Gleysols

Mollic Gleysols

0.0

0.91

575

Umbric Gleysols

Calcic Gleysols

Haplic Vertisols (Eutric)

0.1

0.95

568

Haplic Vertisols

Haplic Kastanozems

Haplic Gypsisols

0.1

0.98

565

Haplic Gypsisols

Aric Regosols

Haplic Solonetz

0.1

0.92

539

Gleyic Solonetz

Solodic Planosols

Calcic Gleysols

0.0

0.90

0.91

514

Umbric Gleysols

Mollic Gleysols

Haplic Nitisols (Rhodic)

0.1

0.94

0.95

492

Haplic Ferralsols

Haplic Acrisols

Haplic Fluvisols (Eutric)

0.1

0.91

465

Haplic Fluvisols

Haplic Ferralsols

Haplic Lixisols (Chromic)

0.0

0.97

441

Haplic Lixisols

Haplic Ferralsols

Calcic Vertisols

-0.1

0.93

437

Calcic Vertisols

Haplic Vertisols

Calcic Kastanozems

-0.1

0.95

415

Haplic Kastanozems

Haplic Luvisols

Leptic Regosols

0.1

0.96

404

Petric Calcisols

Haplic Luvisols

Haplic Luvisols (Chromic)

-0.1

0.93

396

Haplic Luvisols

Haplic Ferralsols

Haplic Solonchaks

0.1

0.95

383

Haplic Solonchaks (Sodic)

Haplic Solonchaks

(Continued )

SoilGrids250m: Global gridded soil information

PLOS ONE | DOI:10.1371/journal.pone.0169748

February 16, 2017

23 / 40

Table 4. (Continued)

Name

ME (%)

TPR

AUC

N

1st class

2nd class

Luvic Chernozems

-0.1

0.93

377

Haplic Kastanozems

Luvic Phaeozems

Acric Ferralsols

0.0

0.97

0.98

371

Haplic Ferralsols

Acric Ferralsols

Fibric Histosols

0.0

0.96

371

Fibric Histosols

Haplic Acrisols

Calcic Luvisols

0.0

0.92

369

Haplic Cambisols

Haplic Luvisols

Calcic Chernozems

0.0

0.94

358

Calcic Chernozems

Haplic Cambisols

Aluandic Andosols

0.0

0.97

341

Aluandic Andosols

Haplic Cambisols

Luvic Calcisols

0.0

0.95

322

Haplic Calcisols

Haplic Kastanozems

Protic Arenosols

0.0

0.99

322

Protic Arenosols

Haplic Leptosols

Haplic Albeluvisols

-0.1

0.97

321

Haplic Albeluvisols

Haplic Cambisols

Mollic Solonetz

-0.1

0.94

312

Gleyic Solonetz

Haplic Kastanozems

Haplic Acrisols (Ferric)

-0.1

0.95

310

Haplic Acrisols

Haplic Cambisols

Haplic Planosols (Eutric)

0.0

0.89

310

Haplic Podzols

Haplic Acrisols

Haplic Gleysols (Eutric)

-0.1

0.93

306

Haplic Gleysols

Haplic Acrisols

Ferralic Cambisols

0.0

0.93

301

Haplic Ferralsols

Haplic Acrisols

Cryic Histosols

-0.1

0.99

299

Cryic Histosols

Haplic Cryosols

Gleyic Solonetz

-0.2

0.94

282

Mollic Solonetz

Haplic Solonetz

Haplic Cambisols (Humic)

0.0

0.92

0.93

273

Haplic Cambisols

Haplic Acrisols

Leptic Phaeozems

0.1

0.97

267

Leptic Phaeozems

Haplic Luvisols

Haplic Regosols (Dystric)

0.0

0.87

262

Haplic Cambisols

Haplic Acrisols

Haplic Leptosols (Eutric)

0.0

0.93

261

Haplic Leptosols

Haplic Calcisols

Acric Plinthosols

0.1

0.97

251

Haplic Acrisols

Haplic Ferralsols

Hemic Histosols

0.0

0.97

250

Hemic Histosols

Albic Luvisols

Endogleyic Cambisols

0.0

0.87

249

Haplic Cambisols

Haplic Acrisols

Haplic Cambisols (Chromic)

0.0

0.93

242

Haplic Cambisols

Haplic Ferralsols

Vertic Cambisols

0.1

0.88

241

Haplic Ferralsols

Haplic Cambisols

Leptic Luvisols

0.2

0.96

0.97

229

Haplic Luvisols

Leptic Phaeozems

Solodic Planosols

0.0

0.96

222

Haplic Solonetz

Haplic Kastanozems

Hypoluvic Arenosols

0.0

0.96

205

Hypoluvic Arenosols

Haplic Arenosols

Leptic Cambisols

0.1

0.94

199

Haplic Luvisols

Petric Calcisols

Umbric Ferralsols

0.1

0.96

192

Haplic Ferralsols

Haplic Acrisols

Gleyic Podzols

0.1

0.95

176

Gleyic Podzols

Haplic Acrisols

Turbic Cryosols

0.0

0.99

1.00

168

Haplic Cryosols

Turbic Cryosols

Vitric Andosols

0.1

0.97

166

Haplic Cambisols

Aluandic Andosols

Haplic Acrisols (Humic)

0.0

0.96

164

Haplic Acrisols

Haplic Cambisols

Haplic Fluvisols (Arenic)

0.0

0.98

163

Haplic Fluvisols

Ferralic Arenosols

Stagnic Luvisols

0.0

0.93

163

Haplic Cambisols

Haplic Luvisols

Mollic Leptosols

0.0

0.90

162

Petric Calcisols

Haplic Leptosols

Haplic Acrisols (Alumic)

0.0

0.98

156

Haplic Acrisols

Haplic Ferralsols

Plinthic Acrisols

0.0

0.94

152

Haplic Acrisols

Plinthic Acrisols

Calcic Solonetz

0.0

0.93

149

Haplic Calcisols

Haplic Kastanozems

Haplic Ferralsols (Xanthic)

0.0

0.96

146

Haplic Ferralsols

Haplic Acrisols

Vertic Luvisols

-0.1

0.93

0.94

140

Haplic Cambisols

Haplic Luvisols

Haplic Lixisols (Ferric)

-0.1

0.96

134

Haplic Lixisols

Haplic Acrisols

Mollic Vertisols

0.0

0.96

133

Mollic Vertisols

Haplic Cambisols

Haplic Solonchaks (Sodic)

-0.1

0.97

130

Haplic Solonchaks

Haplic Arenosols

Sapric Histosols

-0.1

0.90

128

Haplic Cambisols

Fibric Histosols

Haplic Ferralsols (Rhodic)

-0.1

0.96

125

Haplic Ferralsols

Haplic Acrisols

(Continued )

SoilGrids250m: Global gridded soil information

PLOS ONE | DOI:10.1371/journal.pone.0169748

February 16, 2017

24 / 40

organic layer as well. Since the thickness of the organic horizon of these soil profiles is not
reported, their vertical coordinates could not be corrected. There are many situations like this
that require careful analysis of harmonization steps, so that any serious over or under-estima-
tion can be avoided.

It is also fundamentally important that we do not limit ourselves to legacy soil profile data

only. The soil science community needs to actively begin investing in collecting new soil pro-
file field observations, especially in the previously mentioned ecological and climatic zones
that have been under-sampled. For example, the AfSIS project (

http://africasoils.net

) has spent

already half a decade on collecting new samples for Africa. We believe that there is great poten-
tial in undertaking various types of feature space distribution analysis (see e.g. Minasny et al.
[

62

] and Fitzpatrick et al. [

63

]) and optimizing new sampling of additional soil profiles using,

for example, Latin Hypercube sampling principles. By adding only a few hundred new points
that are carefully allocated in extrapolation areas, the accuracy of predictions is likely to
improve more rapidly than if we double the number of points in areas already well repre-
sented. Collection of the new samples could even be implemented via crowd labour or crowd-
sourcing systems so that also local soil surveyors / enthusiasts could get involved (we are
currently testing using MySoil observations contributed by non-specialists, kindly donated to
SoilGrids by the British Geological Survey).

Table 4. (Continued)

Name

ME (%)

TPR

AUC

N

1st class

2nd class

Calcic Gypsisols

-0.4

0.96

124

Calcaric Regosols

Haplic Calcisols

Haplic Cambisols (Sodic)

-0.1

0.98

120

Haplic Cambisols

Ferralic Arenosols

Haplic Calcisols (Sodic)

-0.4

0.98

115

Haplic Calcisols

Haplic Cambisols

Haplic Fluvisols (Dystric)

-0.2

0.93

107

Haplic Fluvisols

Haplic Ferralsols

Haplic Gleysols (Dystric)

-0.3

0.91

100

Haplic Gleysols

Haplic Ferralsols

Gypsic Solonchaks

-0.4

0.98

98

Haplic Gypsisols

Calcaric Regosols

Haplic Luvisols (Ferric)

-0.2

0.95

98

Haplic Luvisols

Haplic Lixisols

Haplic Arenosols (Calcaric)

-0.3

0.94

97

Haplic Arenosols

Haplic Calcisols

Umbric Albeluvisols

-0.2

0.99

1.00

97

Umbric Albeluvisols

Haplic Albeluvisols

Alic Nitisols

-0.1

0.98

70

Haplic Acrisols

Alic Nitisols

Haplic Andosols

-0.2

0.93

67

Aluandic Andosols

Haplic Luvisols

Haplic Planosols (Dystric)

0.0

0.92

62

Ferralic Arenosols

Haplic Ferralsols

Luvic Stagnosols

0.0

0.99

61

Haplic Cambisols

Gleyic Luvisols

Haplic Umbrisols

0.0

0.93

57

Haplic Cambisols

Haplic Acrisols

Albic Arenosols

-0.1

0.93

54

Haplic Acrisols

Haplic Arenosols

Lixic Plinthosols

-0.2

0.94

49

Haplic Ferralsols

Haplic Acrisols

Leptic Umbrisols

-0.4

0.97

40

Haplic Luvisols

Haplic Leptosols

Petric Durisols

-0.3

1.00

39

Petric Durisols

Haplic Phaeozems

Cutanic Alisols

-0.4

0.98

34

Haplic Cambisols

Cutanic Alisols

Endogleyic Planosols

-0.2

0.93

34

Haplic Acrisols

Haplic Luvisols

Haplic Regosols (Sodic)

-0.5

0.97

34

Haplic Vertisols

Leptic Regosols

Luvic Planosols

-0.4

0.88

29

Haplic Luvisols

Haplic Cambisols

Calcic Histosols

-0.9

0.94

18

Haplic Acrisols

Haplic Gleysols

Vetic Acrisols

-1.6

0.90

15

Haplic Acrisols

Haplic Ferralsols

Histic Albeluvisols

-2.3

1.00

13

Umbric Albeluvisols

Fibric Histosols

Vitric Cryosols

-2.3

1.00

13

Vitric Cryosols

Haplic Cambisols

doi:10.1371/journal.pone.0169748.t004

SoilGrids250m: Global gridded soil information

PLOS ONE | DOI:10.1371/journal.pone.0169748

February 16, 2017

25 / 40

Improving the modeling framework

A major improvement from SoilGrids1km to SoilGrids250m is that we now consistently use
machine learning techniques to generate predictions. In the previous version of SoilGrids we
used various types of (Generalized) Linear Models in combination with natural splines to
model soil property-depth relationships, but this resulted in soil property-depth relationships
that were the same across the globe, which is unrealistic and suboptimal. To tackle such prob-
lems we now use dominantly tree-based models—random forest and gradient tree boosting—
to account for local relationships between soil variables and covariates.

Fig 2

(left) shows that,

Fig 9. Maps of scaled Shannon Entropy index (

Eq 5

) for USDA and WRB soil classification maps.

doi:10.1371/journal.pone.0169748.g009

SoilGrids250m: Global gridded soil information

PLOS ONE | DOI:10.1371/journal.pone.0169748

February 16, 2017

26 / 40

indeed, predictions produced using tree-based models adjust locally to observed values. The
current version of SoilGrids is thus, we contend, able to better represent both global and local
patterns.

It has already been demonstrated that random forest can outperform linear models, espe-

cially in being able to better represent complex non-linear relationships in large data sets [

64

–

66

]. Likewise, gradient tree boosting has already won several Kaggle.com competitions (Kaggle

is a platform for predictive modeling and analytics competitions on which companies and
researchers post their data and statisticians and data miners from all over the world compete
to produce the best models). However, tackling the complexities of data size has been a major
challenge. In the case of SoilGrids, the regression matrices had up to one million point pairs
with over 150 covariates, hence their size and complexity well exceeds what can be handled
with desktop computers. Ultimately, we decided to primarily rely on three

R

packages—

caret

[

49

ranger

[

47

] and

xgboost

48

]—that have proven to be capable of processing

huge raster stacks. By using these three open source packages and a single dedicated server

Fig 10. Example of scaled Shannon Entropy index for USDA and WRB soil classification maps with a zoom in on USA state Illinois near the
city of Chicago. This figure uses the same legend as used in

Fig 9

.

doi:10.1371/journal.pone.0169748.g010

SoilGrids250m: Global gridded soil information

PLOS ONE | DOI:10.1371/journal.pone.0169748

February 16, 2017

27 / 40

(current costs of about $800 / month) we were able to optimize and fit all models needed to
generate SoilGrids within a few hours, and to generate all predictions for the entire world
within 12 days.

Machine learning (ML) greatly simplifies model fitting: basically, a soil surveyor does not

need to suggest or impose any relationships—the analyst only needs to list a target variable
and covariates, and machine learning does the ‘magic’ of optimizing model parameters. On
the one hand this is an attractive property because using the ML framework for global soil
mapping allows mapping hundreds of soil variables in parallel with little human interaction.
On the other hand it has also risks and limitations:

• ML is sensitive to noise and errors in the data. Even a few typos in the input values can result

in significant blunders in output maps,

• The computational intensity of ML, when compared to fitting linear models or similar, is an

order of magnitude greater. As the number of training points grows, the computational load

Fig 11. List of some remote sensing data of relevance for global soil mapping projects (i.e. with a near to global
coverage and with remote sensing technology of interest to soil mapping). Landsat 8 is part of the Landsat Data
Continuity Mission (LDCM) maintained by NASA and the United States Geological Survey (USGS). ALOS Global Digital
Surface Model is a product of the Japanese Aerospace Exploration Agency. Sentinel–1,2 is the Earth observation
mission developed by the European Space Agency as part of the Copernicus Programme. WorldDEM

™

is a commercial

product distributed by Airbus Defence and Space.

doi:10.1371/journal.pone.0169748.g011

SoilGrids250m: Global gridded soil information

PLOS ONE | DOI:10.1371/journal.pone.0169748

February 16, 2017

28 / 40

grows exponentially. At some stage, it becomes currently infeasible and overly expensive to
compute predictions using machine learning,

• Extrapolation of models fitted using ML remains risky. Without using pseudo-points to fill-

in data gaps in feature space for some parameters, machine learning can potentially produce
worse maps (on average for most of the soil mask) than linear models,

• Because the sampling locations are clearly biased towards agricultural areas, and because

most of the training points come from the developed world (especially USA), it is very well
possible that SoilGrids predictions are significantly biased in undersampled parts of the
world. In principle, the best solution to this problem is to continuously add more points
from undersampled areas, especially in Africa (tropical soils and wetlands) and the Russian
Federation,

• With ML approaches it is difficult to derive spatially explicit measures of the prediction accu-

racy. We calculated accuracy measures using 10–fold cross-validation, but these are only
global measures.

• ML approaches have a high degree of “black box” modeling and it is difficult to incorporate

knowledge of soil forming processes in the prediction algorithm. But perhaps we can also
learn from ML models by closer inspection and interpretation of how dominant covariates
influence soil property and soil class predictions.

Could machine learning put soil mappers out of work? Probably not. Solid knowledge of

soil science, spatial statistics and/or geostatistics in projects such as SoilGrids is needed more
than ever. For example, it is clear that in order to improve SoilGrids, more focus will need to
be put on improving the feature space representation (adding extra samples) and on improv-
ing visualization and interpretation of complex relationships. Such improvements are not pos-
sible without understanding principles of spatial sampling and soil-environment relationships.
Expert knowledge on soil-landscape relations and soil distribution remains important to eval-
uate the results and assess if predicted spatial patterns make sense from a pedological view-
point. Even though the existing machine learning methods have proven to show improved
predictive performance, much work remains to make them more robust, less sensitive to blun-
ders, incorporate soil-landscape process knowledge and make them more suited for input data
of variable accuracy.

With the current version of SoilGrids, we have also not yet adequately addressed the prob-

lem of vertical soil stratigraphy. At this stage, we remain unable to properly model how some
soil horizons show smooth transition of soil properties, and some show clear and abrupt dis-
continuities (as in geological layers on a meso-scale). In the next update of SoilGrids we hope
to improve modeling and prediction of occurrence of diagnostic soil horizons (e.g. Histic,
Nitic, Albic etc) in 3D, so that transitions between horizons can be represented more
accurately.

Preparation and conversion of soil class input data could also be much improved. Several

research groups [

67

68

] are now looking into automating soil classification (i.e. by using auto-

mated or semi-automated soil classification software). Eberhardt [

69

], for example, demon-

strated using German soil profile data that soil classification can be completely automatized.
Future versions of SoilGrids could also try to derive soil classes by applying exact rules per
pixel, instead of trying to predict them from point data. This might be an ambitious project—
often the classification systems (keys and rules of classification) can be very detailed and
require a comprehensive combination of diagnostic properties, laboratory data, soil-moisture
and temperature regimes, etc. in order to deduce the correct classification. This is without

SoilGrids250m: Global gridded soil information

PLOS ONE | DOI:10.1371/journal.pone.0169748

February 16, 2017

29 / 40

considering the sensitivity of such classifiers to data gaps and uncertainties. Incorporating
uncertainty into such complex soil classification algorithms is yet another challenge. So far, we
have managed to produce global maps of the scaled Shannon Entropy index (

Fig 9

) that clearly

indicate under-represented areas. A sensible approach to improving predictions of soil types
would be to set the sampling intensity proportional to the Shannon Entropy index or
completely focus on areas where the Shannon Entropy index is

>

80%. In that sense, there

seems to be slightly more work needed for the WRB classification system than for the USDA
system.

We have also so far explicitly avoided trying to model posterior distributions of target var-

iables, i.e. map uncertainty for each soil variable. Although tools for modeling uncertainty in
ML methods already exist (see e.g. Meinshausen [

70

]), these are hundreds of times more

computationally intensive and will probably need to be re-implemented in some high-per-
formance computing infrastructure. One future objective is to implement a framework to
model uncertainties of all predictions using a robust statistical framework, such as quantile
regression forests, but this might be highly challenging, especially when the data volumes
grow larger.

Another opportunity for improvement lies in using spatiotemporal modeling [

71

http://ncsslabdatamart.sc.egov.usda.gov/

72

] vs

purely spatial modeling. Stockmann et al. [

2

] recently made progress in modeling global soil

organic carbon dynamics, mainly using time-series of MODIS land cover images, but numer-
ous challenges remain:

• There might not be enough well-distributed soil profile data in the time-domain that support

fitting of spatiotemporal (and/or dynamic) models. As we move back further in the past,
there are fewer and fewer observations, so potential time-domain gaps are possibly an order
of magnitude more serious than spatial data gaps,

• Some soil properties such as soil water content, soil temperature, and even soil nutrients,

change not simply within seasons, but also within weeks or days. At this stage, global fitting
of spatiotemporal models for such variables that vary at short time scales might remain unat-
tainable (until new global soil monitoring networks are established),

• Legacy soil profile data exhibit a significant noise (diversity of methods, laboratories) so that,

for example for soil organic carbon, where temporal dynamics are slow, it will be difficult to
detect real changes in time in a situation where the signal-to-noise ratio is low,

• It is almost impossible to properly validate spatiotemporal predictions produced for past

periods of time. There are very few and sparse validation soil data collected using objective
probability sampling designs (as described in Brus et al. [

60

]). Eventually, we might never

know how accurate our models are in predicting the past status of soil from 50 or 100 years
ago. One possible solution to this problem is linking soil science more directly with paleon-
tology and archeology, but this will probably not work for all soil variables.

Predicting at resolutions finer than 250 meter

Because the algorithms and software we have used in this work are already optimized for
processing large data, this opens a possibility to further speed up model fitting and predic-
tion and to generate predictions at ever finer resolutions.

Fig 11

identifies some new remote

sensing data land products of relevance to global soil mapping. Note that some remote sens-
ing products, such as Landsat 8 and ASTER (distributed as scenes), require significant pro-
cessing capacities before they can be assembled and prepared for use in global soil mapping.
Nevertheless, considering the amount of remote sensing data available publicly today, we

SoilGrids250m: Global gridded soil information

PLOS ONE | DOI:10.1371/journal.pone.0169748

February 16, 2017

30 / 40

anticipate that the Open Source software used in this work will soon (12–24 months) be able
to support generation of 30 m resolution SoilGrids, provided that enough resources exist to
cover the costs of preparing soil covariates and producing global predictions at these fine
resolutions.

Presently, the biggest challenges for upgrading SoilGrids to finer resolution are the

resources required to prepare all required remote sensing input data and computational capac-
ity needed to make fine resolution predictions globally. The software seems to be much less of
a problem. Although

R

has been often criticized for not being suited for large GIS layers, our

experience with SoilGrids has convinced us that, with proper combination of parallelization
and tiling of objects, and by using packages implemented in

C++

or similar, equally efficient

computing can be achieved with

ranger

and

xgboost

(hence within

R

) or by using soft-

ware such as

h2o

(based on

Java

). The remaining bottleneck of

R

we experienced was the

size of models produced using random forest—the objects often exceeded 5–10 GiB and as
such require significant RAM during predictions. Such memory problems in

R

could possibly

be solved via the following two strategies:

1. Disk caching: by using the

ff

or a similar packages to save the forests on disk,

2. Efficient tree representation: transform trees to a simpler structure with the same output.

In the case of random forest, the number of trees required for a given accuracy depends on

the number of rows and columns, i.e. the number of observations (

n

) and covariates (

p

). Usu-

ally, for many rows only few trees are required, while for

p

n

problems (for example in

genetics) many more trees are needed. It should be generally fine to reduce the number of
trees to fewer than 300 but this could be at the expense of loss in accuracy. Lopes [

73

] shows a

framework, based on bootstrapping, to detect an optimal number of trees given some error
threshold. For example, in many cases, even 150 trees is sufficient to achieve stable results after
which a trade-off between computation time and accuracy offers no additional advantages.
We have not tried fine-tuning the number of trees per property (we consistently use 300 trees
as a practical compromise between precision and computing time) because this would have
been an additional load to the project.

Another serious challenge to producing finer resolution SoilGrids is the current lack of

adequately detailed geological data, i.e. data to represent the underlying lithology and min-
eralogy. We have thus far used the Global Lithological Map (GLiM) [

35

] as the key layer to

represent parent material, but this layer is probably even coarser than 1 km resolution
remote sensing covariates, and still contains numerous artifacts such as country/state bor-
ders. Although the OneGeology initiative is of obvious interest to global soil mapping proj-
ects, it has not, so far, delivered any globally consistent and complete information on parent
material. Likewise, the latest most accurate DEM of the world (WorldDEM

™) is an order of

magnitude more accurate and more detailed than the SRTM DEM [

74

] and as such would

be an ideal covariate for many regional and global soil mapping projects. However, it will
likely remain a commercial product available to larger business only (civil engineering and
mineral exploration), and hence of limited use to global soil mappers. In that sense, USA’s
NASA and USGS, with its MODIS, Landsat and similar civil-applications missions will
likely remain the main source of spatial covariate data to support global soil mapping
initiatives.

Other potentially useful covariates for predicting soil properties and classes could be maps

of paleolithic i.e. pre-historic climatic conditions of soil formation, e.g. glacial landscapes and
processes, past climate conditions and similar. These could likely become significant predic-
tors of many current soil characteristics. Information on pre-historic climatic conditions and

SoilGrids250m: Global gridded soil information

PLOS ONE | DOI:10.1371/journal.pone.0169748

February 16, 2017

31 / 40

land use is unfortunately often not available, especially not at detailed cartographic scales,
although there are now several global products that represent, for example, dynamics of land
use / changes of land cover (see e.g. HYDE data set by Klein et al. [

75

]) through the past 1500+

years. As the spatial detail and completeness of such pre-historic maps increases, they will
become potentially interesting covariates for global soil modeling.

Merging global and local: A system for automated soil image fusion

SoilGrids is not expected to be as accurate or relevant as locally produced maps and models
that make use of considerably greater amounts of local point data and finer local covariates.
This is especially the case for OECD countries that can draw upon orders of magnitude
more soil profile data than were used in this work (for illustration, it is estimated that
German Federal agencies alone have in possession 2–3 million complete soil profiles).
Comparison of SoilGrids with similar national or continent-wide products shows that there
is a general match in spatial patterns for many physical and chemical soil properties,
although there are still substantial differences (

Fig 12

). This indicates that promising possi-

bilities exist for further combination of local and global predictions (see further
discussion).

For both Tasmania and California, SoilGrids seems to show somewhat smoother predic-

tions, with some smoothing of higher and lower values, which is especially visible in the cross-
histogram scatter plots (

Fig 12

). SoilGrids tends to overestimate soil pH for parts of Tasmania

covered with rainforests mainly. There were not many ground observations to support the pre-
diction models for those areas, hence some systematic deviation could be expected and will
likely occur in other similar areas as well. We did not run a systematic comparison of values
for all soil properties, but

Fig 12

indicates that merging SoilGrids250m with 100m resolution

predictions using higher density of local soil profiles could help to gradually improve accura-
cies locally and to fill gaps in locally generated predictions.

Mulder et al. [

10

] correctly recognized that, in many areas in the world, locally produced

predictions of soil properties could likely be significantly more accurate than SoilGrids. Our
hope is, nevertheless, that SoilGrids250m will be used by national and regional soil data pro-
duction teams with, or as a supplement to, local data, and that ultimately most users will use
merged (ensemble) global-local predictions for final decision making. We especially recom-
mend the following two frameworks for combining global and local data:

1. SoilGrids predictions as covariate layers for producing finer resolution local predictions of

soil properties (i.e. as an input for downscaling),

2. Ensemble predictions = SoilGrids + local soil spatial prediction models combined.

Option 2, i.e. produce ensemble predictions for smaller areas for which finer resolution

and/or higher quality soil covariates are available, is possibly the most attractive option consid-
ering that local and global predictions can then be generated independently. In that sense, Soil-
Grids could also be considered to be just one (the coarsest) component of a global soil
variation curve (

Fig 13

). But how many components to use to represent soil variation? Are two

components enough? How to optimally merge components where the accuracy is unknown
(not enough ground data for validation)? These will be areas of further research. In that con-
text, Malone et al. [

77

] recently made progress in testing and developing methods for merging

predictions from polygon-based maps and maps derived using spatial predictions. However,
running such models in an automated way for large areas (i.e. a system for an automated soil
image fusion) might take years before an operational system for global soil data fusion is fully
functional.

SoilGrids250m: Global gridded soil information

PLOS ONE | DOI:10.1371/journal.pone.0169748

February 16, 2017

32 / 40

Conclusions

Soil has long been considered one of the least developed global environmental layers with data
available only at coarse resolutions and with limited accuracy [

78

,

79

]. ISRIC—World Soil

Information has a vision and a mission to produce soil information and map products that are
globally complete and consistent, scientifically robust, open, transparent and reproducible,
continuously improved and updated, easy to discover and access, easy to use and meaningful
to users. With this next generation SoilGrids250m we hope to continue to demonstrate prog-
ress in the production and distribution of improved global soil map products and to motivate,
especially non-soil scientists, to use these new soil data in their models and spatial planning,
i.e. directly as input for generation of soil functional properties and agro-ecological variables
and indicators to support decision making. With its Open Data license and web-services, we
aim to serve quality soil information freely and universally for science, society and a sustain-
able future.

Fig 12. Comparison between predicted soil pH: (above) SoilGrids (our predictions) for part of California and predictions based on the SSURGO
data set (for 0–200 cm depth interval) developed by the National Cooperative Soil Survey, (below) SoilGrids (our predictions) for Tasmania and
predictions based on the Soil and Landscape Grid of Australia [

76

] (for 0–5 cm depth interval). The correlation coefficients between the two data

sources are 0.79 and 0.71, respectively. Crosses on the map indicate soil profiles used for generating SoilGrids.

doi:10.1371/journal.pone.0169748.g012

SoilGrids250m: Global gridded soil information

PLOS ONE | DOI:10.1371/journal.pone.0169748

February 16, 2017

33 / 40

We have demonstrated, using a series of cross-validation tests, that the new version of

SoilGrids represents a significant improvement upon the previous products at 1 km resolu-
tion, especially in terms of spatial detail and attribute accuracy. Future work is required to
determine if these improvements in accuracy could also help produce more accurate Global
Gridded Crop Models (GGCMs) that allow for more reliable estimates of impact of climate
change and land degradation on food production [

8

]. Data accessibility problems with Soil-

Grids have also been addressed: SoilGrids are now available for viewing in fusion with satel-
lite imagery via the data portal SoilGrids.org (

Fig 14

). SoilGrids rasters can also be

downloaded via FTP for smaller areas; at point locations through the SoilInfo App and the
REST SoilGrids. There should be fewer and fewer obstacles for ecologists, agronomists,
hydrologists, climatologists, foresters and spatial planners to discover, obtain and use soil
data in their daily work.

Fig 13. SoilGrids can be considered the ‘coarsest’ component of the global soil variation ‘signal’ curve. Other components, e.g. finer
products based on local / more detailed 250–100 m resolution imagery, could be added to produce a merged product.

doi:10.1371/journal.pone.0169748.g013

SoilGrids250m: Global gridded soil information

PLOS ONE | DOI:10.1371/journal.pone.0169748

February 16, 2017

34 / 40

Acknowledgments

SoilGrids are based on numerous soil profile data sets that have been kindly contributed by
various national and international agencies: the USA National Cooperative Soil Survey Soil
Characterization database (

) and profiles from the

USA National Soil Information System (

http://soils.usda.gov/technical/nasis/

), Land Use/Land

Cover Area Frame Survey (LUCAS) Topsoil Survey database [

80

] (

http://eusoils.jrc.ec.europa.

eu/content/lucas-2009-topsoil-data

), Africa Soil Profiles database [

81

], Australian National

Soil Information by CSIRO Land and Water [

82

83

], Mexican National soil profile database

http://dx.doi.org/10.1016/j.gfs.2015.

84

] provided by the Mexican Instituto Nacional de Estadı´stica y Geografı´a / CONABIO, Bra-

zilian national soil profile database [

85

] provided by the University of S

ão Paulo, Chinese

National Soil Profile database [

86

] provided by the Institute of Soil Science, Chinese Academy

of Sciences, soil profile archive from the Canadian Soil Information System [

87

] and Forest

Ecosystem Carbon Database (FECD), ISRIC-WISE [

88

], The Northern Circumpolar Soil Car-

bon Database [

26

], eSOTER profiles [

89

], SPADE [

90

], Unified State Register of soil resources

RUSSIA (Version 1.0. Moscow—2014), National Database of Iran provided by the Tehran
University, points from the Dutch Soil Information System (BIS) prepared by Wageningen
Environmental Research, and others. We are also grateful to USA’s NASA, USGS and USDA
agencies for distributing vast amounts of remote sensing data (especially MODIS land prod-
ucts and SRTM DEM), and to the Open Source software developers of the packages

ranger

,

Fig 14. Basic design and functionality of SoilGrids.org: Soil web-mapping browser that provides interactive viewing of 3D soil layers.
Reference administrative data, basic functionality and output data license of SoilGrids.org are primarily based on OpenStreetMap.

doi:10.1371/journal.pone.0169748.g014

SoilGrids250m: Global gridded soil information

PLOS ONE | DOI:10.1371/journal.pone.0169748

February 16, 2017

35 / 40

xgboost

,

caret

,

raster

,

GDAL

,

SAGA GIS

and similar, and without which SoilGrids

would most likely not be possible. Every effort has been made to trace copyright holders of the
materials used to produce SoilGrids spatial predictions. Should we, despite all our efforts have
overlooked contributors please contact ISRIC and we shall correct this unintentional omission
without any delay and will acknowledge any overlooked contributions and contributors in
future updates.

Author Contributions

Conceptualization: TH JMJ GBMH MK WS MAG RAM JGBL IW BK.

Data curation: MRG WS XG RAM NHB JGBL SM.

Formal analysis: TH GBMH MK WS MNW BK.

Funding acquisition: TH GBMH.

Investigation: TH JMJ GBMH MRG WS BK.

Methodology: TH JMJ GBMH MRG MK WS MNW BBM RAM SM BK.

Project administration: TH JMJ MK.

Resources: JMJ MK AB BBM ER.

Software: TH JMJ MRG MK AB WS MNW ER BK.

Supervision: TH GBMH.

Validation: MRG XG MAG RV RAM NHB JGBL ER IW SM.

Visualization: TH JMJ AB MK ER.

Writing – original draft: TH JMJ GBMH WS MNW BBM RAM NHB IW BK.

Writing – review & editing: TH JMJ GBMH MK WS MNW BBM RAM NHB IW BK.

References

1.

Scharlemann JPW, Tanner EVJ, Hiederer R, Kapos V. Global soil carbon: understanding and managing
the largest terrestrial carbon pool. Carbon Management. 2014; 5(1):81–91. doi:

10.4155/cmt.13.77

2.

Stockmann U, Padarian J, McBratney A, Minasny B, de Brogniez D, Montanarella L, et al. Global soil
organic carbon assessment. Global Food Security. 2015; 6:9–16.

07.001

.

3.

Aksoy E, Yigini Y, Montanarella L. Combining Soil Databases for Topsoil Organic Carbon Mapping in
Europe. PLoS ONE. 2016; 11(3):1–17. doi:

10.1371/journal.pone.0152098

4.

Shani U, Ben-Gal A, Tripler E, Dudley LM. Plant response to the soil environment: An analytical model
integrating yield, water, soil type, and salinity. Water resources research. 2007; 43(8). doi:

10.1029/

2006WR005313

5.

Shepherd KD, Shepherd G, Walsh MG. Land health surveillance and response: A framework for evi-
dence-informed land management. Agricultural Systems. 2015; 132:93–106. doi:

10.1016/j.agsy.2014.

09.002

6.

Pogson M, Hastings A, Smith P. Sensitivity of crop model predictions to entire meteorological and soil
input datasets highlights vulnerability to drought. Environmental Modelling & Software. 2012; 29(1):37–
43. doi:

10.1016/j.envsoft.2011.10.008

7.

Montanarella L, Vargas R. Global governance of soil resources as a necessary condition for sustainable
development. Current Opinion in Environmental Sustainability. 2012; 4(5):559–564. doi:

cosust.2012.06.007

SoilGrids250m: Global gridded soil information

PLOS ONE | DOI:10.1371/journal.pone.0169748

February 16, 2017

36 / 40

8.

Folberth C, Skalsky´ R, Moltchanova E, Balkovi

č

J, Azevedo LB, Obersteiner M, et al. Uncertainty in soil

data can outweigh climate impact signals in global crop yield simulations. Nature Communications.
2016; 7. doi:

10.1038/ncomms11872

9.

Hengl T, Mendes de Jesus J, MacMillan RA, Batjes NH, Heuvelink GBM, Ribeiro E, et al. SoilGrids1km
—Global Soil Information Based on Automated Mapping. PLoS ONE. 2014; 9(e105992). doi:

PMID:

Mulder V, Lacoste M, Richer-de Forges A, Martin M, Arrouays D. National versus global modelling the
3D distribution of soil organic carbon in mainland France. Geoderma. 2016; 263:16–34. doi:

10.1016/j.apsoil.2015.06.018

geoderma.2015.08.035

11.

Griffiths RI, Thomson BC, Plassart P, Gweon HS, Stone D, Creamer RE, et al. Mapping and validating
predictions of soil bacterial biodiversity using European and national scale datasets. Applied Soil Ecol-
ogy. 2016; 97:61–68. doi:

12.

FAO/IIASA/ISRIC/ISS-CAS/JRC. Harmonized World Soil Database (version 1.2). Rome: FAO; 2012.

13.

Domisch S, Amatulli G, Jetz W. Near-global freshwater-specific environmental variables for biodiversity
analyses in 1 km resolution. Scientific data. 2015; 2. doi:

10.1038/sdata.2015.73

PMID:

26647296

14.

Chen J, Chen J, Liao A, Cao X, Chen L, Chen X, et al. Global land cover mapping at 30 m resolution: A
POK-based operational approach. ISPRS Journal of Photogrammetry and Remote Sensing. 2015;
103:7–27. doi:

10.1016/j.isprsjprs.2014.09.002

15.

IUSS Working Group WRB. World reference base for soil resources 2006: a framework for international
classification, correlation and communication. World soil resources reports No. 103. Rome: Food and
Agriculture Organization of the United Nations; 2006.

16.

U S Department of Agriculture. Keys to Soil Taxonomy. 12th ed. U.S. Government Printing Office;
2014.

17.

Arrouays D, Grundy MG, Hartemink AE, Hempel JW, Heuvelink GBM, Young Hong S, et al. Chapter
Three—GlobalSoilMap: Toward a Fine-Resolution Global Grid of Soil Properties. In: Sparks DL, editor.
Soil carbon. vol. 125 of Advances in Agronomy. Academic Press; 2014. p. 93–134.

18.

Nelson DW, Sommers LE. Total carbon, organic carbon, and organic matter. In: Page AL, Miller RH,
Keeney DR, editors. Methods of soil analysis, Part 2. 2nd ed. Agron. Monogr. 9. Madison, WI: ASA
and SSSA; 1982. p. 539–579.

19.

Shangguan W, Hengl T, de Jesus JM, Yuan H, Dai Y. Mapping the global depth to bedrock for land sur-
face modeling. Journal of Advances in Modeling Earth Systems. 2016; p. n/a–n/a.

20.

FAO. Guidelines for soil profile description. 4th ed. Rome: Food and Agriculture Organization of the
United Nations; 2006.

21.

Shaw C, Bhatti JS, Sabourin KJ. An Ecosystem Carbon Database for Canadian Forests. Information
report. Edmonton, Alberta: Canadian Forest Service, Northern Forestry Centre; 2005.

22.

Krasilnikov P, Marti JJI, Arnold R, Shoba S, editors. A Handbook of Soil Terminology, Correlation and
Classification. Earthscan LLC; 2009.

23.

Lei SA. Soil Properties of the Kelso Sand Dunes in the Mojave Desert. The Southwestern Naturalist.
1998; 43(1):47–52.

24.

Raup B, Racoviteanu A, Khalsa SJS, Helm C, Armstrong R, Arnaud Y. The GLIMS geospatial glacier
database: a new tool for studying glacier change. Global and Planetary Change. 2007; 56(1):101–110.
doi:

10.1016/j.gloplacha.2006.07.018

25.

Mantel S, Wo¨sten H, Verhagen J, et al. Biophysical land suitability for oil palm in Kalimantan, Indonesia.
ISRIC Report 2007/01. ISRIC—World Soil Information; 2007.

26.

Hugelius G, Tarnocai C, Broll G, Canadell JG, Kuhry P, Swanson DK. The Northern Circumpolar Soil
Carbon Database: spatially distributed datasets of soil coverage and soil carbon storage in the northern
permafrost regions. Earth System Science Data. 2013; 5(1):3–13.

27.

Danielson JJ, Gesch DB. Global multi-resolution terrain elevation data 2010 (GMTED2010). Open-File
Report 2011–1073. U.S. Geological Survey; 2011.

28.

Conrad O, Bechtel B, Bock M, Dietrich H, Fischer E, Gerlitz L, et al. System for Automated Geoscientific
Analyses (SAGA) v. 2.1.4. Geoscientific Model Development. 2015; 8(7):1991–2007. doi:

10.5194/

gmd-8-1991-2015

29.

Savtchenko A, Ouzounov D, Ahmad S, Acker J, Leptoukh G, Koziana J, et al. Terra and Aqua MODIS
products available from NASA GES DAAC. Advances in Space Research. 2004; 34(4):710–714. doi:

10.1016/j.asr.2004.03.012

30.

Mira M, Weiss M, Baret F, Courault D, Hagolle O, Gallego-Elvira B, et al. The MODIS (collection V006)
BRDF/albedo product MCD43D: temporal course evaluated over agricultural landscape. Remote Sens-
ing of Environment. 2015; 170:216–228. doi:

10.1016/j.rse.2015.09.021

SoilGrids250m: Global gridded soil information

PLOS ONE | DOI:10.1371/journal.pone.0169748

February 16, 2017

37 / 40

31.

Wan Z. MODIS land surface temperature products users’ guide. ICESS, University of California; 2006.

32.

Hall DK, Riggs GA. Accuracy assessment of the MODIS snow products. Hydrological Processes. 2007;
21(12):1534–1547.

33.

Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A. Very high resolution interpolated climate sur-
faces for global land areas. International Journal of Climatology. 2005; 25:1965–1978. doi:

10.1002/joc.

1276

34.

Huffman GJ, Bolvin DT. GPCP version 2 combined precipitation data set documentation. vol. 1. Labora-
tory for Atmospheres, NASA Goddard Space Flight Center and Science Systems and Applications;
2009.

35.

Hartmann J, Moosdorf N. The new global lithological map database GLiM: A representation of rock
properties at the Earth surface. Geochemistry, Geophysics, Geosystems. 2012; 13(12):n/a–n/a. doi:

10.1029/2012GC004370

36.

Sayre R, Dangermond J, Frye C, Vaughan R, Aniello P, Breyer S, et al. A new map of global ecological
land units—an ecophysiographic stratification approach. Washington, DC: USGS / Association of
American Geographers; 2014.

37.

Fan Y, Li H, Miguez-Macho G. Global patterns of groundwater table depth. Science. 2013; 339
(6122):940–943. doi:

10.1126/science.1229881

PMID:

23430651

38.

Giri C, Ochieng E, Tieszen LL, Zhu Z, Singh A, Loveland T, et al. Status and distribution of mangrove
forests of the world using earth observation satellite data. Global Ecology and Biogeography. 2011; 20
(1):154–159. doi:

10.1111/j.1466-8238.2010.00584.x

39.

Pelletier JD, Broxton PD, Hazenberg P, Zeng X, Troch PA, Niu GY, et al. A gridded global data set of
soil, immobile regolith, and sedimentary deposit thicknesses for regional and global land surface model-
ing. Journal of Advances in Modeling Earth Systems. 2016; 8. doi:

10.1002/2015MS000526

40.

Jenny H. Factors of soil formation: a system of quantitative pedology. Dover books on Earth sciences.
Dover Publications; 1994.

41.

Hijmans RJ, van Etten J. raster: Geographic data analysis and modeling; 2012. Available from:

http://

CRAN.R-project.org/package=raster

.

42.

Pebesma EJ, Bivand RS. Classes and methods for spatial data in R. R news. 2005; 5(2):9–13.

43.

Mitchell T, Developers G. Geospatial Power Tools: GDAL Raster & Vector Commands. Locate Press;
2014.

44.

Carroll ML, Townshend JR, DiMiceli CM, Noojipady P, Sohlberg RA. A new global raster water mask at
250 m resolution. International Journal of Digital Earth. 2009; 2( 4):291–308. doi:

10.1080/

17538940902951401

45.

Rabus B, Eineder M, Roth A, Bamler R. The shuttle radar topography mission—a new class of digital
elevation models acquired by spaceborne radar. ISPRS Journal of Photogrammetry and Remote Sens-
ing. 2003; 57(4):241–262. doi:

10.1016/S0924-2716(02)00124-7

46.

Venables WN, Ripley BD. Modern applied statistics with S. 4th ed. New York: Springer-Verlag; 2002.

47.

Wright MN, Ziegler A. ranger: A Fast Implementation of Random Forests for High Dimensional Data in
C++ and R. Journal of Statistical Software. 2016; p. 18.

48.

Chen T, Guestrin C. XGBoost: A Scalable Tree Boosting System. ArXiv e-prints. 2016;.

49.

Kuhn M. Building Predictive Models in R Using the caret Package. Journal of Statistical Software. 2008;
28(1):1–26.

50.

Sollich P, Krogh A. Learning with ensembles: How over-fitting can be useful. In: Proceedings of the
1995 Conference. vol. 8; 1996. p. 190.

51.

Knaus J, Porzelius C, Binder H, Schwarzer G. Easier parallel computing in R with snowfall and sfClus-
ter. The R Journal. 2009; 1(1):54–59.

52.

Bauer-Marschallinger B, Sabel D, Wagner W. Optimisation of global grids for high-resolution remote
sensing data. Computers & Geosciences. 2014; 72:84–93. doi:

10.1016/j.cageo.2014.07.005

53.

Cohen J. Weighted kappa: Nominal scale agreement provision for scaled disagreement or partial credit.
Psychological bulletin. 1968; 70(4):213. doi:

10.1037/h0026256

PMID:

19673146

54.

Sing T, Sander O, Beerenwinkel N, Lengauer T. ROCR: visualizing classifier performance in R. Bioin-
formatics. 2005; 21(20):3940–3941. doi:

10.1093/bioinformatics/bti623

PMID:

16096348

55.

Sing T, Sander O, Beerenwinkel N, Lengauer T.

ROCR

: Visualizing the Performance of Scoring Classifi-

ers; 2016. Available from:

http://CRAN.R-project.org/package=ROCR

.

56.

Shannon CE. Communication in the presence of noise. Proceedings of the IRE. 1949; 37(1):10–21. doi:

10.1109/JRPROC.1949.232969

SoilGrids250m: Global gridded soil information

PLOS ONE | DOI:10.1371/journal.pone.0169748

February 16, 2017

38 / 40

57.

Borda M. Fundamentals in Information Theory and Coding. Springer Berlin Heidelberg; 2011. Available
from:

https://books.google.nl/books?id=Lyte2yl1SPAC

.

58.

Kempen B. Updating soil information with digital soil mapping [PhD thesis]. Wageningen University.
Wageningen; 2011. Available from:

http://edepot.wur.nl/187198

.

59.

Wright MN, Ziegler A. ranger: A Fast Implementation of Random Forests for High Dimensional Data in
C++ and R. ArXiv e-prints. 2015;.

60.

Brus DJ, Kempen B, Heuvelink GBM. Sampling for validation of digital soil maps. European Journal of
Soil Science. 2011; 62(3):394–407. doi:

10.1111/j.1365-2389.2011.01364.x

61.

Ribeiro E, Batjes N, Leenaars JGB, van Oostrum A, editors. Towards the standardization and harmoni-
zation of world soil data. ISRIC Report 2015/03. Wageningen, the Netherlands: ISRIC—World Soil
Information; 2015.

62.

Minasny B, McBratney AB. A conditioned Latin hypercube method for sampling in the presence of ancil-
lary information. Computers & Geosciences. 2006; 32(9):1378–1388. doi:

10.1016/j.cageo.2005.12.

009

63.

Fitzpatrick MC, Gotelli NJ, Ellison AM. MaxEnt versus MaxLike: empirical comparisons with ant species
distributions. Ecosphere. 2013; 4(5):art55–. doi:

10.1890/ES13-00066.1

64.

Statnikov A, Wang L, Aliferis C. A comprehensive comparison of random forests and support vector
machines for microarray-based cancer classification. BMC Bioinformatics. 2008; 9(1):319. doi:

PMID:

Kuhn M, Johnson K. Applied predictive modeling. Springer; 2013.

66.

Hengl T, Heuvelink GBM, Kempen B, Leenaars JGB, Walsh MG, Shepherd KD, et al. Mapping Soil
Properties of Africa at 250 m Resolution: Random Forests Significantly Improve Current Predictions.
PLoS ONE. 2015; 10(e0125814). doi:

10.1371/journal.pone.0125814

PMID:

26110833

67.

La´ng V, Fuchs M, Waltner I, Miche´li E. Soil taxonomic distance, a tool for correlation: As exemplified by
the Hungarian Brown Forest Soils and related WRB Reference Soil Groups. Geoderma. 2013;
192:269–276. doi:

10.1016/j.geoderma.2012.07.023

68.

Brevik EC, Calzolari C, Miller BA, Pereira P, Kabala C, Baumgarten A, et al. Soil mapping, classification,
and pedologic modeling: History and future directions. Geoderma. 2016; 264:256–274. doi:

10.1016/j.spasta.2015.04.001

geoderma.2015.09.022

69.

Eberhardt E. Deriving WRB soil types using National German Soil Survey Guidelines. In: NCSS
National Conference. NCSS Standards. Asheville, NC: United States Department of Agriculture USDA-
NRCS; 2011. p. 1–2.

70.

Meinshausen N. Quantile regression forests. The Journal of Machine Learning Research. 2006; 7:983–
999.

71.

Kilibarda M, Hengl T, Heuvelink G, Gra¨ler B, Pebesma E, Per

č

ec Tadi

ć

M, et al. Spatio-temporal inter-

polation of daily temperatures for global land areas at 1 km resolution. Journal of Geophysical
Research: Atmospheres. 2014; 119(5):2294–2313.

72.

Gasch C, Hengl T, Gra¨ler B, Meyer H, Magney T, Brown D. Spatio-temporal interpolation of soil water,
temperature, and electrical conductivity in 3D+T: The Cook Agronomy Farm data set. Spatial Statistics.
2015; 14(Part A):70–90. doi:

73.

Lopes ME. Measuring the Algorithmic Convergence of Random Forests via Bootstrap Extrapolation.
Davis CA: Department of Statistics, University of California; 2015.

74.

Baade J, Schmullius C. Uncertainties of a TanDEM-X derived Digital Surface Model 2014—A case
study from the Roda catchment, Germany. In: Geoscience and Remote Sensing Symposium
(IGARSS), 2014 IEEE International; 2014. p. 4327–4330.

75.

Klein Goldewijk K, Beusen A, Van Drecht G, De Vos M. The HYDE 3.1 spatially explicit database of
human-induced global land-use change over the past 12,000 years. Global Ecology and Biogeography.
2011; 20(1):73–86. doi:

10.1111/j.1466-8238.2010.00587.x

76.

Kidd D, Webb M, Malone B, Minasny B, McBratney A. Eighty-metre resolution 3D soil-attribute maps for
Tasmania, Australia. Soil Research. 2015; 53(8):932–955.

77.

Malone BP, Minasny B, Odgers NP, McBratney AB. Using model averaging to combine soil property
rasters from legacy soil maps and from point data. Geoderma. 2014; 232:34–44. doi: