Reference soil groups map of Ethiopia based on legacy data and machine learning-technique: EthioSoilGrids 1.0

Up-to-date digital soil resource information and its comprehensive understanding are crucial to supporting crop production and sustainable agricultural development. Generating such information through conventional approaches consumes time and resources, and is difficult for developing countries. In Ethiopia, the soil resource map that was in use is qualitative, dated (since 1984), and small scaled (1 : 2 M), which limit its practical applicability. Yet, a large legacy soil profile dataset accumulated over time and the emerging machine-learning modeling approaches can help in generating a high-quality quantitative digital soil map that can provide better soil information. Thus, a group of researchers formed a Coalition of the Willing for soil and agronomy data-sharing and collated about 20 000 soil profile data and stored them in a central database. The data were cleaned and harmonized using the latest soil profile data template and 14 681 profile data were prepared for modeling. Random forest was used to develop a continuous quantitative digital map of 18 World Reference Base (WRB) soil groups at 250 m resolution by integrating environmental covariates representing major soil-forming factors. The map was validated by experts through a rigorous process involving senior soil specialists or pedologists checking the map based on purposely selected district-level geographic windows across Ethiopia. The map is expected to be of tremendous value for soil management and other land-based development planning, given its improved spatial resolution and quantitative digital representation.

1 Introduction

Soils are important resources that support the development and production of various economic, social, and ecosystem services, and are useful in climate change mitigation and adaptation (Baveye et al., 2016). Data on the physical and chemical characteristics of soils and their spatial distribution are needed to define and plan their functions over time and space, which are important steps toward sustainable use and management of soils (Elias, 2016; Hengl et al., 2017).

In Ethiopia, soil surveys and mapping have been conducted at various scales with varying scopes, approaches, methodologies, qualities, and levels of detail (Abayneh, 2001; Abayneh and Berhanu, 2007; Berhanu, 1994; Elias, 2016; Zewdie, 2013). The most recent countrywide digital soil mapping efforts focused primarily on soil characteristics (Ali et al., 2020; Iticha and Chalsissa, 2019; Tamene et al., 2017), although soil class maps are equally important for allocating a particular soil unit for specific use (Leenaars et al., 2020a; Wadoux et al., 2020). Many attempts have been made to improve digital soil information systems (Hengl et al., 2021, 2017, 2015; Poggio et al., 2021). However, the initiatives were based on limited and unevenly distributed soil profile data (e.g., 1.15 soil profiles per 1000 km² for Ethiopia), which restricts the accuracy and applicability of the products.

In Ethiopia, thousands of soil profile data have been collected since the 1960s (Erkossa et al., 2022), but these data were scattered across different institutions and individuals (Ali et al., 2020). Furthermore, countrywide quantitative and gridded spatial soil-type information does not exist (Elias, 2016). The Ethiopian Soil Information System (EthioSIS) project attempted to develop a countrywide digital soil map focusing on topsoil characteristics, including plant nutrient content, but overlooked soil resource mapping (Ali et al., 2020; Elias, 2016), despite a strong need for a high-resolution soil resource map (Mulualem et al., 2018).

Ethiopia has an area of about 1.14×10⁶ km² consisting of varied environments, making its soils extremely heterogeneous. Capturing the heterogeneity using conventional soil survey and mapping approaches is an expensive and time-consuming endeavor (Hounkpatin et al., 2018). This can be circumvented by using available legacy soil profile data accumulated over decades and by tapping into the potential of advanced analytical techniques to develop high-resolution digital soil maps (Hounkpatin et al., 2018; Kempen, 2012, 2009). Therefore, the objectives of this study were to (i) develop a national legacy soil profile dataset that can be used as an input for various digital soil mapping exercises, and (ii) generate an improved 250 m digital Reference Soil Groups (RSGs) map of Ethiopia.

2 Methods

2.1 The study area

The study area covered the entire area of Ethiopia (1.14×10⁶ km²) located between 3 and 15° N, and between 33 and 48° E (Fig. 1). The topography of the country is marked by a large altitudinal variation, ranging from 126 m below sea level at Dalol in the northeast to 4620 m at Ras Dashen Mountain in the northwest (Billi, 2015; Enyew and Steeneveld, 2014). Ethiopia's wide range of topography, climate, parent material, and land use types created conditions for the formation of different soil types (Abayneh, 2005; Berhanu and Ochtman, 1974; Donahue, 1972; Mesfin, 1998; Nyssen et al., 2019; Virgo and Munro, 1978; Zewdie, 2013, 1999). More than 33 % of the country is covered by the central, upper, and highland complex (Abegaz et al., 2022), which embraces Africa's most prominent mountain system (Hurni, 1998).

Figure 1 Location map of Ethiopia (inset) and overview map of Esri World Topographic Map.

The country's complex topography strongly determines both rainfall and temperature patterns, by modifying the influence of the large-scale ocean–land–atmosphere pattern, thus creating diverse localized climates. Spatially, rainfall is characterized by a general decreasing trend in the direction from west to east, north, northeast, south and southeast. The lowlands in the southeast and northeast, covering approximately 55 % of the country's land area, are characterized by arid and semi-arid climates. Annual rainfall ranges from less than 300 mm in the southeastern and northwestern lowlands to over 2000 mm in the southwestern highlands (southern portion of the western highlands). The eastern lowlands get rain twice a year, in April–May and October–November, with two dry periods in between. The total annual precipitation in this region varies from less than 500 to 1000 mm. The driest of all regions is the Denakil Plain, which receives less than 500 mm of rain and sometimes none (Fazzini et al., 2015). Temperatures are also greatly influenced by the rapidly changing altitude, and the mean monthly values vary from ∼ 35 °C in the northeast lowlands to less than 7.5 °C over the north and central highlands.

The country is characterized by a wide variety of geological formations (Abyneh, 2005; Alemayehu et al., 2014; Elias, 2016; Zewdie, 2013). These include (i) recent and old volcanic activities; (ii) the highlands consisting of igneous rocks (mainly basalts); (iii) steep-sided valleys characterized by strong colluvial and alluvial deposits; (iv) metamorphic rocks exposed by denudation process; and (v) various sedimentary rocks such as limestone and sandstone in the relatively lower areas.

Diverse biophysical factors affecting the spatial distribution of vegetated land cover which in turn, both as single and combined factors, result in diverse soil types and properties across Ethiopia's landscapes (Hurni, 1998; Nyssen et al., 2019; WLRC-AAU, 2018). The spatiotemporal vegetation cover of the country has been characterized by a long history of land use and land cover changes (WLRC-AAU, 2018). In terms of the type and spatial coverage of major land use and land cover classes, woody vegetation (forest, woodland, and shrub and bush lands) covers about 57 % of the country in accordance with the national 2016 map (WLRC-AAU, 2018). This is followed by cultivated land (20 %) and grasslands (12 %). Barren lands are estimated to cover about 1/10 of the area of the country while other minor lands with ecological significance (i.e., wetlands, water bodies, and sub-afro-alpine and afro-alpine) cover about 1.2 % of the country's land mass.

2.2 Legacy soil profile data collation and preparation

The soil profile data generated over decades through various soil survey missions were kept in a variety of formats with limited accessibility. There has been no institution with a mandate to coordinate the generation, collation, harmonization, and sharing of soil profile data. This led to the formation of a group of individuals and institutions who were willing to exchange soil and agronomy data. Established in 2018, the group known as the Coalition of the Willing (CoW) was committed to addressing the challenges posed by the lack of soil and agronomy data access and sharing in the country (Tamene et al., 2021).

The CoW conducted a national soil and agronomy data ecosystem mapping which revealed that a plethora of legacy soil resource datasets exist across different institutions and individuals (Ali et al., 2020). The assessment also revealed that a sizable proportion of the data holders were willing to share the data in their custody, provided that some regulations were put in place to administer the data. The CoW developed and approved internal data-sharing guidelines (CoW, 2020), and facilitated data collation campaigns, which involved both formal and informal approaches to data holders.

Through a data collation campaign, soil profile data collected between the 1970s and 2021 were acquired from over 88 diverse sources (Ali et al., 2020; Tamene et al., 2021). Initially, 8000 profile data points were collated and subjected to improved modeling techniques to create a provisional WRB reference soil group map of Ethiopia. This was presented to various partners and data-holding institutions to demonstrate the power of data sharing. This created awareness and enabled us to mobilize and collate over 20 000 legacy soil profile data. These data were then added to the national data repository.

The data had varying levels of completeness in terms of soil field and environmental descriptions and laboratory analysis. These required a rigorous expert-based quality assessment and standardization before being compiled into a harmonized format. The expanded version of the Africa Soil Profile (AfSP) database (Leenaars et al., 2014) template was used for standardizing and harmonizing the data. Out of the collated soil profile data, 14 681 georeferenced data points were extracted based on completeness and cleanness for the purposes of modeling. The cleaned soil profile data set contained, at least, the reference soil group (RSG) nomenclature as outlined in the WRB legend. While the original soil profile records were set in different coordinate systems, all were projected into the adopted standard georeferencing system, namely, WGS84, decimal degrees in the QGIS (3.20.2) environment (QGIS Development Team, 2021). To verify their position, soil profile locations were plotted using a standard WGS84 coordinate system to verify that points matched the site description, geomorphological settings, and at the very least the source project boundary outline.

The accuracy of the data depends on the quality and reliability of the survey data themselves, which in turn requires expert knowledge and experience in soil description and classification (Leenaars et al., 2020a). In this study, data cleaning, validation, reclassification, and verification were carried out by a team of prominent national pedologists and soil surveyors, including those involved in the generation of some of the soil profile data themselves (Fig. 2).

Figure 2 Schematic presentation of data acquisition and workflow.

In addition, the Ministry of Agriculture (MoA) soil survey and mapping experts and other volunteers validated the legacy soil profile observations. This led to the reclassification of the soil types as deemed necessary. Such validation and reclassification involved re-examining the geomorphological setup of the soil profile locations using Google Earth as well as reviewing the site and soil descriptions and the corresponding laboratory data, and reviewing the proposed soil type. The harmonized datasets in the database were used as input soil profile data for modeling and mapping IUSS WRB reference soil groups.

2.3 Preparation and selection of environmental covariates

2.3.1 Covariate acquisition and preparation

In order to develop spatially continuous soil class and/or type maps, data on environmental covariates that represent directly or indirectly the soil-forming factors have to be integrated with soil profile data (Hengl and MacMillan, 2019). Environmental covariates are spatially explicit proxies of soil-forming factors based on the soil–environment relationship (McBratney et al., 2003; Shi et al., 2018). Acquisition and preparation of covariates represent a crucial step in digital soil mapping using machine-learning algorithms (McBratney et al., 2003). In this study, 68 potential candidate environmental variables representing soil-forming factors (climate, organisms, relief, parent material, and time) were derived from diverse remote sensing products and thematic maps (Hengl and MacMillan, 2019; McBratney et al., 2003).

Relief and topography-related covariates were derived from a 90 m Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) (Vågen, 2010). Climate-related variables including long-term mean, minimum, maximum, and standard deviation temperature as well as precipitation data for the period between 1983 and 2016 (Dinku et al., 2014) were acquired from Enhancing National Climate Services (ENACTS-NMA) initiatives with 4 km resolutions (Dinku et al., 2014). Moderate-resolution imaging spectroradiometer (MODIS) imagery raw bands and derived indices (Vågen, 2010) were downloaded from USGS EarthExplorer (https://earthexplorer.usgs.gov/, last access: 12 November 2021) to represent vegetation-related factors. National geological (Tefera et al., 1996) and land use and land cover (WLRC-AAU, 2018) thematic maps of Ethiopia were gathered to represent parent material and organisms, respectively.

Downscaling (disaggregating) or upscaling (aggregating) of rasters was also performed to match the target resolution. A 250 m spatial resolution was chosen to accommodate both the spatial resolution of the major covariate inputs and make it applicable for large-scale analysis. All layers were masked for buildings and water bodies by the national boundary of Ethiopia and a stacked layer was created using the raster package (R Core Team, 2020) to extract covariate values at the locations of soil profiles. One-hot encoding using the dummyVars function available in Caret package (Kuhn, 2008) was used to pre-process and convert categorical covariates into a binary vector. Each element of the binary vector represents the presence or absence of that category. One-hot encoding is beneficial because it enables machine-learning algorithms to interpret categorical variables as numerical features. The covariate pre-processing, visual inspection for inconsistencies, and resampling to a target grid of 250 m were conducted in QGIS [3.20.2] (QGIS Development Team, 2021), SAGA GIS [7.8.2] (Conrad et al., 2015) and R [version 4.05] (R Core Team, 2020) software packages. All input data were projected to a common Lambert azimuthal equal-area projection with the latitude of origin at 8.65 and center of meridian at 39.64, which is the center point for Ethiopia. This projection was selected since it is effective in minimizing area distortions over land. Each covariate was adjusted to have an identical spatial resolution, extent, and projection using two resampling methods. Continuous covariates were resampled using the bilinear spline method, whereas categorical covariates were resampled using the nearest neighbor method.

2.3.2 Covariate selection

Selecting an optimal set of covariates to effectively represent the soil–environment relationship is a key step in digital soil mapping (DSM) since improper selection of covariates will affect the quality of model outputs (Shi et al., 2018). In this study, near-zero variance assessment was conducted using the nearZeroVar function available in the R caret package (Kuhn, 2008) to identify and remove environmental variables that have little or no variance. In addition, preliminary random forest model training was performed to assess and identify covariates having high variable importance. After expert judgment, a total of 27 environmental variables (24 continuous and 3 categorical) were selected for modeling and predicting RSGs.

2.4 Modeling and mapping soil types or reference soil groups

2.4.1 Model tuning and quantitative evaluation

In digital soil mapping, machine-learning techniques have been extensively used to determine the relationship between soil types and environmental variables (McBratney et al., 2003). Many machine-learning models were developed in the past decades for digital soil mapping to spatially predict soil classes based on existing soil data and soil-forming environmental covariates (Heung et al., 2016). Random forest (RF), a tree-based ensemble method, is one of the most promising machine-learning techniques available for digital soil mapping (Breiman, 2001; Heung et al., 2016). RF has gained popularity due to its high overall accuracy and has been widely used in predictive soil mapping (Brungard et al., 2015; Hengl et al., 2018). Examples of the main strengths of the RF model are its ability to handle numerical and categorical data without any assumption of the probability distribution, and its robustness against nonlinearity and overfitting (Breiman, 2001; Svetnik et al., 2003). While building the RF model, data were split into training (80 %) and testing (20 %) components using random sampling for training the model and evaluating its performance, respectively (Kuhn, 2008). Hyper-parameter optimization and repeated cross-validation on the training dataset were performed for optimal model application using the ranger method of the Caret package. The three tuning parameters for ranger method are mtry, splitrule, and .min.node.size. Generally this function is used to tune the parameters in modeling in an automated fashion, as this will automatically check all the possible tuning parameters and return the optimized parameters on which the model gives the best accuracy. Model tuning was performed with a repeated 10-fold cross-validation procedure applying multiple combinations of hyper-parameters for the ranger method. This is a fast implementation of RF particularly suited for high-dimensional data (Wright and Ziegler, 2017). Then the number of covariates used for the splits (mtry), splitting rules (splitrule), and minimum node size (min.node.size) were optimized. The parameter ntree was adjusted to 1000 in the model, and mtry values (10, 15, 20), min.node.size values (5, 10, 15), and splitrule values (“variance”, “extratrees”, and “maxstat”) were fed for the optimization procedure. The accuracy of the testing dataset was related to the model performance for the new dataset, indicating the capacity of the model to predict at the unsampled location. A confusion matrix was also used to calculate a cross-tabulation of observed and predicted classes with associated statistics, i.e., producer's accuracy and user's accuracy.

2.4.2 Software and computational framework

In this study, various open-source software packages that provide a comprehensive set of tools and diverse capabilities were used for data preparation, analysis, and visualization. Data pre-processing and preparation were performed using QGIS (QGIS Development Team, 2021) and SAGA GIS (Conrad et al., 2015). For statistical analysis and machine-learning modeling, R (R Core Team, 2020) and relevant libraries were installed on a Windows server, 2016 standard with 250 GB of working memory, to handle the challenges associated with large-scale data processing and analysis.

2.4.3 Expert evaluation of spatial patterns of the beta-version soil map

Visual inspection of the DSM output over the terrain was used to identify abnormalities and assess how effectively it depicts landscape components (Rossiter et al., 2022). For this, we employed an expert-based qualitative assessment of the model output. This technique was used to complement model-based accuracy assessment and confirm agreement soil specialists or pedologists checking the map based on purposely selected district-level geographic windows across Ethiopia, representing different agro-ecological zones known to have diverse soil occurrences, and that were familiar to the panel of experts. Accordingly, an expert validation workshop was conducted using the first version of the reference soil groups (RSGs) map. About 45 multi-disciplinary scientists including soil surveyors, pedologists, geologists, and geomorphologists were drawn from national and international research, development, and higher-learning institutions to review the draft RSG map in plenary discussions. This was followed by breakout sessions where groups of experts evaluated the map based on their experience and knowledge of soil–landscape relations of the country and examined geographic windows.

Most importantly, disagreements regarding RSG occurrence and patterns of the modeling outputs across topo-sequences and contrasting soil-forming factor sequences were identified and discussed. Further, inferences on parts of the DSM framework that require improvement were recommended. After finalizing the evaluation at the group-level assessment, each group presented the results in the plenary followed by a discussion to get feedback from other participants. Following the plenary discussions, the participants created a group of six senior pedologists to work on the recommendations including changing the quality mask layer, validating the additional data obtained during the event, and assessing the re-modeling outputs.

After the second model was re-run, the group of senior pedologists together with geospatial experts re-evaluated the output using the selected districts based on the feedback from the first review, which was mainly on areas where there were “minor” and “major” concerns. Consequently, some improvements were made, e.g., in the areas where Vertisols, Fluvisols, and Leptosols were overestimated. Further, underestimated RSGs (Alisols, Solonetz, Planosols, Acrisols, Lixisols, Phaeozems, and Gleysols) showed a slight increase in area coverage and pattern improvements. However, the total area of Leptosols and Cambisols increased from the first run due to the partial exclusion of the mask layer used in the first round of modeling. The mask layer used in the first run was criticized for quality issues as it excluded significant soil areas and due to its weakness in capturing non-soil areas such as rock outcrops, salt flats, swamps, and sand dunes. Nevertheless, the spatial patterns of these soils occurring across previously considered “non-soil areas” were examined by the panel of experts. In parallel, geospatial and soil experts checked the raster map of the RSGs in the GIS environment to ensure areas with “no concern” before re-running the model are kept the same or changes are accepted by the panel of experts. The map from the second run is presented in this paper as EthioSoilGrids version 1.0 product.

3 Results and discussion

3.1 Soil profile datasets

Using the IUSS WRB (2015), the preliminary identified 14 742 georeferenced legacy soil profiles were classified and/or reclassified into 23 RSGs. Nearly 90 % of the soil profile points represented Vertisols, followed by Luvisols, Cambisols, Leptosols, Fluvisols, and Nitisols, which were found to be the dominant soil types in Ethiopia (Fig. 3). The remaining 10 % represented the Regosols, Alisols, Andosols, Arenosols, Calcisols, Solonetz, Lixisols, Phaeozems, Solonchaks, Acrisols, Planosols, Gleysols, Umbrisols, Ferralsols, Gypsisols, Plinthosols, and Stagnosols.

According to this study, about 72 % of the IUSS WRB (2015) RSGs were confirmed to occur in Ethiopia. This reconfirms the characterization of Ethiopia as a land of soil diversity being endowed with a diverse range of soil types (Elias, 2016; Mishra et al., 2004). One of the limitations with legacy soil data in categorical mapping is the imbalanced soil samples, in that all classes are not equally represented (Wadoux et al., 2020). For this study, soil profiles with fewer than 30 observations were objectively excluded from the model after examining the accuracy and spatial distribution of each RSG. Five RSGs (Umbrisols, Ferralsols, Gypsisols, Plinthosols, and Stagnosols) were excluded from the model and the EthioSoilGrids version 1.0 map.

Figure 3 Number of soil profile points per WRB reference soil groups.

After excluding the built-up and water surface areas, the average soil profile density was 13.1 per 1000 km² (Fig. 4), but the actual density varied across the different parts of the country. The variation tends to follow river basins, sub-basins, and agricultural land-use type-based studies from which most of the legacy data were pulled. For instance, in 30 intervention districts of the Capacity Building for Scaling up of Evidence-Based Best Practices in Agricultural Production in Ethiopia (CASCAPE) project, the average profile density was about 87 profiles per 1000 km² for a total area of about 26 830 km² (Leenaars et al., 2020a). Similarly, semi-detailed soil mapping missions in 15 districts conducted through the Bilateral Ethiopia–Netherlands Effort for Food, Income and Trade (BENEFIT)-REALISE project generated about 217 observations per 1000 km² (Leenaars et al., 2020b).

A soil type and depth map compilation and updating mission at a 1 : 250 000 scale by the Water Land Resource Center (WLRC) of Addis Ababa University collated and used about 3949 legacy soil profiles for the entire country (Ali et al., 2020), which is approximately 3.5 profiles per 1000 km². Although the distribution is not even and the eastern lowlands are sparsely represented, the number of data used in this study is 8.5 times higher than the 1712 legacy soil profiles data currently existing in the Africa soil profile database (Batjes et al., 2020; Leenaars et al., 2014).

Figure 4 Spatial distribution of collated legacy soil profile data.

The distribution of the soil profiles across the 32 agro-ecological zones (AEZ) of Ethiopia revealed that all, except two – tepid per-humid mid-highland (0.13 % landmass) and very cold sub-humid sub-afro-alpine to afro-alpine (0.03 % landmass) – were represented by soil profile observations. Furthermore, about 95 % of the profile observations represented 91 % of the AEZ aerial coverage (Appendix A). The distribution of legacy soil profiles varied across AEZs. In general, the top-ranked lowland AEZs with roughly 56 % area coverage were represented by 23 % of the total profile observations, whereas top-ranked highland AEZs with 20 % area coverage received 47 % of profile observations. For instance, warm desert, warm moist, hot arid, and warm sub-moist lowlands with area coverage of around 20 %, 15 %, 11 %, and 10 %, were represented roughly by 3 %, 11 %, 2 %, and 7 % of the total profiles, respectively. Tepid moist mid-highlands (8 % area coverage), tepid sub-humid mid-highlands (7 % area coverage), and tepid sub-moist mid-highlands (5 % area coverage) each were represented by 20 %, 15 %, and 12 % of the profiles, respectively.

3.2 Modeling and mapping

3.2.1 Variable importance

The RSG spatial pattern is primarily influenced by long-term average surface reflectance, flow-based DEM indices, and precipitation. Figure 5 shows variables of importance for determining RSG spatial prediction. The top-ranked variables were (i) long-term MODIS near-infrared (NIR) reflectance, (ii) multiresolution index of valley bottom flatness, (iii) long-term mean day–land surface temperature, (iv) long-term mean soil moisture, (v) standard deviation of long-term precipitation, (vi) long-term mean precipitation, and (vii) topographic wetness index.

MODIS long-term mean spectral signatures showed high relative importance. According to Hengl et al. (2017), accounting for seasonal vegetation fluctuation and inter-annual variations in surface reflectance, long-term temporal signatures of the soil surface, derived as monthly averages from long-term MODIS imagery, were more effective. Furthermore, Hengl and MacMillan (2019) explained that long-term average seasonal signatures of surface reflectance provide a better indication of soil characteristics compared with only a single snapshot of surface reflectance.

The multi-resolution valley bottom flatness index, a DEM-derived topography index, is the second top-ranked covariate driving soil variability across Ethiopia. This hydrological/soil removal and accumulation or deposition index is used to distinguish valley floor and ridgetop landscape positions (Soil Science Division Staff, 2017) greatly responsible for multiple soil-forming processes to operate over a particular landscape, resulting in a wide range of soil development. The influence of topography on spatial soil variation is manifested in every landscape of Ethiopia (Belay, 1997; Mesfin, 1998; Nyssen et al., 2019; Zewdie, 2013).

Long-term daily mean land surface temperature, mean soil moisture, rainfall standard deviation, and mean annual rainfall were among the top-ranked covariates for predicting the spatial variation of RSGs across the country. In Ethiopia, different soil genesis studies revealed that climate has a significant influence on soil development and properties and is, therefore, responsible for the existence of widely varying soils in the country (Abayneh, 2005; Abayneh et al., 2006; Fikru, 1988, 1980; Zewdie, 2013).

Among the most important covariates for predicting RSGs in the Ethiopian highlands are monthly average soil moisture for January (ranked third), long-term average soil moisture (ranked fourth), and monthly average soil moisture for August (ranked fifth) (Leenaars et al., 2020a). In the current study, soil moisture was among the 10 top-ranked covariates in modeling and explaining long-distance soil type variability across the country.

Figure 5 Random forest covariate relative importance for modeling RSGs. Note: prep = precipitation; prep_sd = standard deviation of precipitation; tmax = maximum temperature; tmin = minimum temperature; trange = temperature range; tav_sd = standard deviation of average temperature; pet = potential evapotranspiration; lstd = land surface temperature–day; lstn = land surface temperature–night; soil_moist = soil moisture; soil_temp = soil temperature; DEM  = digital elevation model (elevation); twi = topographic wetness index; aspect = topographic aspect; curv = topographic curvature; conv = topographic convergence index; ls = slope length and steepness factor (ls_factor); morph = terrain morphometry; mrvbf = multiresolution index of valley bottom flatness; slope = slope class (%); ndvi = normalized difference vegetation index (NDVI); evi = enhanced vegetation index (EVI); lulc = land use/land cover; lithology = geology; ref1 = red band; ref2 = near-infrared; ref7 = mid-infrared.

In this study, lithology showed a relatively low influence on soil variability that may be due to the use of a coarse-scale and less detailed lithology map, which may not sufficiently capture the spatial variability of the parent materials.

3.2.2 Model performance

The parameter optimization process resulted in mtry = 20, split rule = extra trees and minimum node size = 5. The overall accuracy of the model was 56.24 % which ranged between 54.43 % and 58.1 % with a 95 % confidence interval. The kappa values based on the internal cross-validation and testing dataset showed that the overall model performance produced using 10-fold cross-validation with the repeated fitting was 48 %. Considering similar area-based digital soil class mapping efforts, the overall accuracy was in line with the accuracies that were typically reported for soil class maps developed with RF models (Leenaars et al., 2020a) and statistical methods (Heung et al., 2016). Table 1 shows the confusion matrix at validation/testing points, i.e., 20 % of the observation. Further, the matrix indicates the producer's accuracy (class representation of observed versus predicted) and user's accuracy were not similar for all RSGs. The map purity is in the order of Lixisols, Calcisols, Alisols, Phaeozems, Vertisols, Andosols, Solonchaks, Fluvisols, Arenosols, Leptosols, Luvisols, Nitisols, and Cambisols. However, Vertisols, Calcisols, and Andosols are the observed classes that are best represented by the map followed by Fluvisols, Alisols, Nitisols, Leptosols, Luvisols, and Cambisols.

Global soil grids at 250 m resolution used machine-learning algorithms to map the global WRB RSGs with map purity and weighted kappa of 28 % and 42 %, respectively (Hengl et al., 2017). The SoilGrids 250 m WRB soil groups/classes prediction output–spatial soil patterns were not evaluated based on expert knowledge while in this study we did an extensive back-and-forth qualitative assessment by a panel of pedologists. The quantitative accuracy in the present study (about 56 %) coupled with an expert-based qualitative evaluation of the predicted maps indicated the development and achievement of a substantially enhanced national product for users of spatial soil resource information. This finding is a step forward and acceptable considering that SoilGrids maps are not expected to be as accurate as locally produced maps and models that use many more local-point data and finer local variables (Mulder et al., 2016). Further, the data and findings in this study can help improve the soil maps of Africa as they partially address the concern by Hengl et al. (2017), who recognized that WRB RSGs modeling in the global SoilGrids 250 m is critically uncertain for parts of Africa. This is mainly attributed to limited access to more local point data by regional and global modeling initiatives, unlike the present study which accessed a large number of legacy soil profile datasets.

Table 1 Confusion matrix of random forest RSG prediction (at validation/testing observations).

3.2.3 Modeling and mapping: EthioSoilGrids version 1.0

The study identified 18 RSGs in Ethiopia, mapped at 250 m resolution (Fig. 6). The model prediction showed that seven soil reference groups including Cambisols, Leptosols, Vertisols, Fluvisols, Nitisols, Luvisols, and Calcisols covered nearly 98 % of the total land area of the country (Fig. 7). Five soil reference groups (Solonchaks, Arenosols, Regosols, Andosols, and Alisols) were estimated to cover about 2 % of the land area, while trace coverages of Solonetz, Planosols, Acrisols, Lixisols, Phaeozems, and Gleysols were also found in some pocket areas.

In terms of spatial distribution, Nitisols and Luvisols dominated the northwestern and southwestern highlands while the southeastern lowlands were dominantly covered by Cambisols, Calcisols, and Fluvisols with some Solonchaks. The Vertisols extensively cover the north and southwestern lowlands along with the Ethiopia–Sudan border areas and central highland plateaus. The probability of occurrence of each RSG was mapped (Appendix C) in each modeling spatial window (i.e., the cell size of 250 m × 250 m). The dominant RSGs were aggregated based on the most probable RSGs in each spatial modeling window. There was high correspondence between the seven top-ranked prediction probabilities and observed soil types as confirmed visually by overlaying observed classes and prediction probabilities.

The overall occurrence and the relative position of each of the RSGs along the topo-sequence and its association with other RSGs agree with previous works (Abayneh et al., 2006; Ali et al., 2010; Abdenna et al., 2018; Asmamaw and Mohammed, 2012; Belay, 2000, 1998, 1997, 1996; Driessen et al., 2001; Elias, 2016; FAO, 1984a; Fikre, 2003; Mitiku, 1987; Mohammed and Belay, 2008; Mohammed and Solomon, 2012; Mulugeta et al., 2021; Nyssen et al., 2019; Sheleme, 2017; Shimeles et al., 2007; Tolossa, 2015; Zewdie, 2013). However, in some cases, the position of the RSGs along the topo-sequence and the association with other RSGs require further investigation. The disparities observed might be attributed to the positional accuracy of legacy point observations, the modeling approach, and most importantly the level of detail and scale/resolution of the environmental variables used in this study. We used the currently available coarse-resolution national geological map and hence soil parent material might be inadequately represented in the model, which probably resulted in irregular RSG sequences. For instance, the main driving factors to establish and explain the soil-landscape variability in the May-Leiba catchment of northern Ethiopia were geology (soil parent material) and different mass movements (Van de Wauw et al., 2008). These factors led to Cambisols–Vertisols catenas on basalt and Regosols–Cambisols–Vertisols catenas on limestone formations. Similar studies identified parent material strongly determines the soil type (e.g., Vertisol, Luvisol, Cambisol) (Nyssen et al., 2019). In general, in areas where there is complex soil diversity and distribution of soils, one of the most important parameters is to identify parent material including effective techniques to capture and delineate mass movement bodies, and human-induced soil erosion and deposition areas (Leenaars et al., 2020a; Nyssen et al., 2019; Van de Wauw et al., 2008).

Figure 6 Major reference soil groups of Ethiopia (EthioSoilGrid V1.0).

Considering the third position of Cambisols in the order of frequency of occurrence of RSGs per point observations (following Vertisols and Luvisols), these soils seem to be over-represented on the map (ranked first) apparently at the expense of Vertisols and Luvisols, and to some extent in places of Leptosols and other RSGs. This might be attributed to the fact that Cambisols create a geographical continuation with Vertisols and/or Luvisols at the lower slopes and Leptosols/Regosols at the higher slopes, suggesting the presence of some bordering soil qualities in respective transitional zones (Ali et al., 2010; Asmamaw and Mohammed, 2012; Sheleme, 2017; Zewdie, 2013).

The proportion of area mapped as Cambisols (34 %) revealed new insights compared with the information from the most cited spatial soil maps: Cambisols ranked second (21 %), second (16 %), fourth (9 %), and fourth (8 %) as reported by Berhanu (1980), FAO (1984b, 1998), and SoilGrids – Hengl et al. (2017), respectively. This might be due to (i) the number and distribution of profile observations, which is more extensive than the previous ones; (ii) the type and level of details of covariates considered; (iii) variations and rearrangements in the keys for classification of the RSGs among soil classification versions used in previous studies and misclassification/confusion of Vertisols with Vertic Cambisols, as legacy soil profile data come from diverse sources.

Figure 7 The area coverage (in %) for the major WRB RSGs. Note: the remaining 10 RSGs-Arenosols (0.44 %), Regosols (0.35 %), Andosols (0.31 %), Alisols (0.16 %), Solonetzs (0.04 %), Planosols (0.04 %), Acrisols (0.02 %), Lixisols (0.02 %), Phaeozems (0.02 %), and Gleysols (0.01 %) were not plotted because of their relatively small area coverage.

3.3 Expert validation of the soil map

Expert knowledge of soil–landscape relations and soil distribution remains important for evaluating the predictive soil mapping results and assessing whether the predicted spatial patterns make sense from a pedological viewpoint (Hengl et al., 2017; Poggio et al., 2021; Rossiter et al., 2022). An important step in qualitative model evaluation is, therefore, expert assessment, whereby professionals with broad experience in soil survey and mapping can evaluate and improve the quality of the soil resource map. This can highlight areas of agreement or concern across the landscape (Rossiter et al., 2022). The expert validation workshop provided useful insights and tangible improvements to the development of the map. While the plenary discussion provided an overview of the approaches followed in developing the map, the group discussions helped to have an in-depth review of the selected polygons of the map assigned to them. Participants were split into five groups (with 8–10 members each) and chose up to 60 polygons representing areas with which at least one of the group members has sufficient information, including data sources. Overall, the groups checked a total of 126 polygons (Fig. 8), which were fairly evenly distributed across the country.

Figure 8 The spatial distribution of districts validated by stakeholders and feedback categories according to the level of concerns raised.

The group members displayed the polygons one by one in a GIS environment and discussed the predicted dominant and associated soil RSGs and labeled them in one of three confirmation categories: (1) confirmed with “no concern”, (2) confirmed with “minor concern”, and (3) confirmed with “major concern”. Confirmation with “no concern” was made when all members of a group agreed on the types, the relative coverage, and the patterns of the predicted soils within the polygon. Confirmation with “minor concern” was made when all or some of the team members agreed on the predicted soil types within the polygons but did not agree on the order of abundance or the probability occurrence of one or two soils including observed spatial patterns. Confirmation with “major concern” was made when all members of the team did not agree on the predicted soil type, or when the presence of another soil type, other than the predicted types, was noted.

All three groups rated the accuracy of the map at 60+ %; of the 126 polygons, they expressed no concern for 63 %, minor concern for 23 %, and major concern for 14 % of the polygons. Furthermore, differences in the prevalence of RSGs and patterns of the modeling outputs across different soil-forming factor sequences, as well as inferences about which areas of the DSM framework still need work, were identified and elaborated on by the expert input and are presented in the subsequent sections.

3.4 Evaluation of results, limitations, and future direction

Up-to-date soil resource spatial information is critically missing at the required scale and extent in Ethiopia. As a result, resource management strategies miss their targets. Furthermore, the absence of such data at a required resolution and extent forced developers of decision support tools to pick and use the data they can access and afford. As a result, model outputs appear more site-specific or representation becomes homogeneous over the very heterogeneous landscapes that exist in reality. On the other hand, in large areas and complex landscapes such as Ethiopia, it is very difficult to address the demand for reasonably accurate and detailed soil-type maps using a conventional approach due to the costs involved and to the resources and time this requires. For instance, given the vastness of the country and the heterogeneous landscapes, a new conventional soil survey mission requires at least 170 000 profile point observations to map the entire terrestrial land mass of Ethiopia at a scale of 1 : 250 000 with at least one observation per square centimeter. Moreover, the soil profile data requirement definitely could have been much higher as we increase the scale of mapping and density of observations. In the present study, machine-learning techniques combined with expert input were implemented to produce a countrywide soil resource map of Ethiopia at reasonably higher accuracy and with less time and cost compared with conventional methods. In addition, rescue, compilations, and standardization of about 14 681 geo-referenced legacy soil profiles that can be included in the National Soil Information System (NSIS) of Ethiopia and the World Soil Information Center will support future national, regional, and global DSM efforts. The approach used here demonstrates the power of data and analytics to map the soil resources of Ethiopia, and the output is an exemplary use case for similar digital content development efforts in Ethiopia and beyond.

Moreover, in this study the quality-monitoring processes and methods were followed to filter dubious soil profiles as well as soil classification and harmonization protocols. Thereafter, the study followed a robust modeling framework and generated new insights into the relative area coverage of WRB RSGs of Ethiopia. In addition, the study provided coherent and up-to-date digital quantitative gridded spatial soil resource information to support the successful implementation of various digital agricultural solutions and decision support tools (DSTs).

The spatially explicit limitation of the present study is revealed by expert-based qualitative evaluation of spatial patterns across objectively selected geographic windows and prominent contrasting landscapes of Ethiopia. This qualitative assessment indicated areas of concern in terms of how well EthioSoilGrids version 1.0 represents soil geography across a mosaic of the country's landscapes. For instance, in the northeastern lowlands of Ethiopia, mainly along the “Denakil” depression, Fluvisols, Cambisols, and Vertisols were found on the map in areas where normally other soil types were expected to occur. In this area, the expected prediction and area coverage of Leptosols has probably been overshadowed by Fluvisols and Cambisols. Similarly, in some parts of western Ethiopia landscapes, the prediction of Vertisols overshadows other RSGs, which resulted in an underestimation of the area coverage of Fluvisols (along the “Akobo”, “Gilo”, and “Baro” rivers and their tributaries) and Alisols. Likewise, in the central parts of northwestern Ethiopia, the prediction of Nitisols was overshadowed by Vertisols and Luvisols, resulting in a likely underestimation of the Nitisols area coverage.

The relatively low model performance and some classification errors in some of the examined geographic windows (e.g., the Denakil depression, along Akobo, Baro, and Gilo rivers and the Somali region) are probably due to the paucity of samples from those areas (Fig. 4), the inadequacy of the dataset by RSGs, and over-representation of the dataset by some RSGs, such as Vertisols, Luvisols, and Cambisols. Balanced datasets are ideal to allow decision tree algorithms to produce better classification but for datasets with uneven class size, the generated classification model might be biased toward the majority class (Hounkpatin et al., 2018; Wadoux et al., 2020). In addition, uncertainty around the quality of the covariates included, not the covariates considered in the modeling process including management, use of validation methods that do not sufficiently control the effect of clustered samples, and small sample size for some RSGs could have possibly biased the modeling results in some geographic areas.

To improve the modeling performance, future studies could explore (i) adding data for under-represented geographic areas, land uses, and covariate spaces; (ii) opportunities to include other covariates (parent material and management) that could capture the variability of the country's heterogeneous landscapes; (iii) dimension reduction of covariates; (iv) use of remedial measures for imbalances in sample sizes; (v) comparing different cross-validation methods; (vi) use of an ensemble modeling approach and/or robust modeling technique that accommodates neighborhood size and connectivity analyses; (vii) use of a better-resolution/quality mask layer to segregate non-soil areas (rock outcrops, salt flats, sand dunes, and water bodies) from mapping areas; and (viii) implementation of quantitative and qualitative comparisons of national, regional, and global legacy soil maps/soil grids with new DSM products in terms of how well DSM products represent soil geography. In addition, future digital soil mapping strategies in Ethiopia may require consideration of new soil sampling missions in under-represented areas; adoption of standard soil sampling, description guidelines, and soil classification systems including soil physicochemical and mineralogical analysis; and a combination of local soil nomenclature/classification systems with RSGs and development of a map of RSGs with qualifiers. At the moment the under-sampled and under-represented areas are the Somali region, the Denakil, and the western and northwestern border areas of Ethiopia (Fig. 4). Despite these limitations, and to the best of our knowledge, the EthioSoilGrids v1.0 product provides the most complete soil information available for Ethiopia.

4 Conclusions

Coherent and up-to-date countrywide digital soil information is essential to support digital agricultural transformation efforts. This study involved collation, cleaning, harmonization, and validation of the legacy soil profile datasets, involving soil scientists with different backgrounds individually and in groups. To develop the 250 m digital soil resource map, a machine-learning modeling approach and expert validation were applied to the harmonized soil database and environmental covariates affecting soil-forming processes. Accordingly, about 20 000 soil profile data were collated, out of which about 14 681 were used for the modeling and mapping of 18 RSGs out of the 23 RSGs identified. Although unevenly distributed, the legacy soil profile data used in the modeling covered most of the agro-ecologies of the country.

Among the 18 RSGs mapped, the highest number of observed (3935) profiles represent Vertisols, followed by Luvisols, Cambisols, and Leptosols, while Gleysols were represented with the lowest number (63) of profiles. The modeling revealed that the most important covariates for predicting RSGs in Ethiopia are MODIS long-term reflectance, multiresolution index of valley bottom flatness, land surface temperature, soil moisture, long-term mean annual rainfall, and wetness index of the landscape.

Our 10-fold spatial cross-validation result showed an overall accuracy of about 56 % with varying accuracy levels among RSGs. The modeling result revealed that seven major soil reference groups including Cambisols (34 %), Leptosols (20 %), Vertisols (18 %), Fluvisols (10 %) Nitisols (7 %), Luvisols (6 %), and Calcisols (3 %) covered nearly 98 % of the total land area of the country, while minor coverage of other RSGs (Solonchaks, Arenosols, Regosols, Andosols, Alisols, Solonetzs, Planosols, Acrisols, Lixisols, Phaeozems, and Gleysols) was also detected in some areas. Compared with the existing soil resource map, the coverage of the first three major soil groups has substantially increased, which is related to the increased availability of soil profile data covering larger areas of the country, implying that these soils were previously underestimated. Cambisols and Vertisols which together represent nearly half of the total land area are relatively young with inherent fertility, suggesting a high agricultural potential for the country. However, given their limitations, these and the other soil types require the implementation of suitable land, water, and crop management techniques to sustainably exploit their potential.

The EthioSoilGrids version 1.0 product from this first countrywide RSGs modeling effort requires complementary activities. These include modeling and mapping that should go beyond RSGs and need to include second-level classifications including principal and supplementary qualifiers. Furthermore, a soil atlas of Ethiopia with details of the soil physicochemical properties needs to be prepared together with the map, which the authors and/or others responsible need to prioritize in their future research endeavors.

Appendix A: Legacy soil profile data distribution

Table A1 Distribution of legacy soil profile data by agroecology zones.

Note: ^a total area of Ethiopia 1.14×10⁶ km²; ^b total number of profiles 14 681.

Appendix B: Environmental covariates

Table B1 List, description, spatial and temporal extent, and source of covariates used in modeling the reference soil groups.

^a Africa Soil Information Service (AfSIS). ^b Fifth generation European Centre for Medium-Range Weather Forecasts (ECMWF) atmospheric reanalysis of the global climate.

Appendix C: Probability of occurrence of reference soil groups

Figure C1 Occurrence probability maps of Cambisols, Leptosols, Vertisols, and Fluvisols.

Figure C2 Occurrence probability maps of Nitisols, Luvisols, and Calcisols.