Visible, near-infrared, and shortwave-infrared spectra as an input variable for digital mapping of soil organic carbon

Original Research Article Visible, near-infrared, and shortwave-infrared spectra as an input variable for digital mapping of soil organic carbon Vahid Khosravi a , * , Asa Gholizadeh a , Radka Kode sov a a , Prince Chapman Agyeman a , Mohammadmehdi Saberioon b , Lubo s Bor uvka a a Department of Soil Science and Soil Protection, Faculty of Agrobiology, Food and Natural Resources, Czech University of Life Sciences Prague, Kamycka 129, Suchdol, Prague, 16500, Czech Republic b Helmholtz Centre Potsdam GFZ German Research Centre for Geosciences, Section 1.4 Remote Sensing and Geoinformatics, Telegrafenberg, Potsdam, 14473, Germany article info Article history: Received 27 April 2023 Received in revised form 2 October 2024 Accepted 4 October 2024 Available online xxx Keywords: SOC modeling and mapping Interpolated spectra Machine learning Regression kriging Uncertainty abstract This study proposes a novel methodology to employ discrete point spectra as input variable for digital mapping of soil organic carbon (SOC). Accordingly, two SOC modeling approaches were used in three agricultural sites in Czech Republic: i) machine learning (ML) including partial least squares regression (PLSR), cubist, random forest (RF), and support vector regression (SVR), and ii) regression kriging (RK) by the combination of ordinary kriging (OK) and PLSR (PLSR-K), cubist (cubist-K), RF (RF-K), and SVR (SVRK). Models were developed on environmental predictor covariates (EPCs) and thirty genetic algorithms (GA)-selected visible, near-infrared, and shortwave-infrared (VNIReSWIR) wavelengths spectra, individually and combined. Thirty rasters were then created using interpolation of the selected spectra and served as the input variables ewith and without EPCs eto test and compare the developed models and SOC predictive maps with each other and with those retrieved from the third approach: iii) kriging using OK of the measured and ML-predicted SOC. The impact of employing selected wavelengths’spectra and EPCs on models' performance was investigated using independent test samples and the uncertainty associated with the produced maps. Using interpolated spectra as the only input variable yielded a relatively acceptable accuracy (Nov a Ves: RMSE ¼0.19%, Údrnice: RMSE ¼0.12%, Klu cov: RMSE ¼0.13%). In comparison, the interpolated spectra coupled with EPCs enhanced the results. Regarding the uncertainty, however, the ML-based SOC maps were more reliable, than RK-based ones. Furthermore, maps produced using both spectra and EPCs showed less uncertainty than those constructed on the individual datasets. ©2024 International Research and Training Center on Erosion and Sedimentation, China Water and Power Press, and China Institute of Water Resources and Hydropower Research. Publishing services by Elsevier B.V. on behalf of KeAi Communications Co. Ltd. This is an open access article under the CC BYNC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). 1. Introduction Soil organic carbon (SOC) is a dynamic property that plays a crucial role in fertility of the agriculture and forest ecosystems (Bhunia et al., 2019). Various soil-related processes and services, including food production and climate change mitigation, are being affected by this soil property (Szatm ari et al., 2021). It is, hence, imperative to continuously monitor and map SOC across the soil landscape. Conventional laboratory-based chemical measurement methods are expensive, laborious, and time-consuming, needing extra chemicals that might be of environmental concern. The traditional SOC polygon maps are also expensive and timedemanding to produce, difficult to update, and without sufficient spatial resolution (Mahmoudzadeh et al., 2020). This has sparked an increasing interest in indirect prediction and digital mapping of the SOC content eas a target variableeusing separated and hybrid implementation of non-geospatial techniques such as simple and multiple linear regression and machine learning (ML) approaches, and geospatial techniques such as geostatistical kriging. As a subgroup of artificial intelligence, ML has recently been under significant attention in soil sciences (Padarian et al., 2020)to quantify the relationship between soil properties, point, and *Corresponding author. E-mail address: [email protected] (V. Khosravi). Contents lists available at ScienceDirect International Soil and Water Conservation Research journal homepage: www.elsevier.com/locate/iswcr https://doi.org/10.1016/j.iswcr.2024.10.002 2095-6339/©2024 International Research and Training Center on Erosion and Sedimentation, China Water and Power Press, and China Institute of Water Resources and Hydropower Research. Publishing services by Elsevier B.V. on behalf of KeAi Communications Co. Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). International Soil and Water Conservation Research xxx (xxxx) xxx Please cite this article as: V. Khosravi, A. Gholizadeh, R. Kode sov aet al., Visible, near-infrared, and shortwave-infrared spectra as an input variable for digital mapping of soil organic carbon, International Soil and Water Conservation Research, https://doi.org/10.1016/ j.iswcr.2024.10.002 imaging spectroscopy, and weather and climate data. ML algorithms can deal with hidden patterns within high dimensional complex datasets as well as non-linear dependencies between variables. Moreover, geostatistical methods, such as kriging, are principally used to quantify changes in soil properties over various distances mostly by means of a semi-variogram as their powerful central tool (Oliver, 1987). During the modeling process, ML techniques account for the deterministic part of the total variations and spatial correlation is disregarded, while kriging unravels the spatially dependent stochastic part and independent use of each model will cause losing one other part (Keskin &Grunwald, 2018). To avoid this problem, a number of hybrid techniques have been developed, including universal kriging (Burgess &Webster, 1980), kriging with external drift (Goovaerts,1997), and regression kriging (RK), which is also called “kriging combined with regression” (Knotters et al., 1995). Kriging with external drift and universal kriging share the same formulation, that estimation of trend and residuals is done in a single system, where the prediction variance is also determined. RK, on the other hand, involves kriging of the residual values after regression is applied to derive a target variable and can be independent of the kriging (Minasny &McBratney, 2007). One main important drawback of kriging is that it produces one uncertain value for each location while cannot address the real uncertainty. This is mainly due to the smoothing effect, making the model's final outputs unreliable to make proper decisions. Simulation techniques can deal with these problems by generating many interpolated surfaces, all able to reproduce the samples spatial features (Giustini et al., 2019), creating better reality representations while removing the smoothing effects by focusing on reproduction of the semi-variogram models or global statistics (Goovaerts, 1997). Accordingly, simulation techniques are increasingly preferred over kriging. In the case of normally distributed continuous data types, the gaussian geostatistical simulation is usually more popular. The simulation is considered as “conditional”, if the simulated values honor the observed values of known sample points at their locations, i.e., are conditional on the observed variables. As an extensively used type of conditional simulation, sequential gaussian simulation (SGS) employs the previously simulated values for simulation of the successive grid points (Goovaerts, 2001). According to the SCORPAN function defined by McBratney et al. (2003), soil maps are dependent on factors of soil (S), climate (C), organisms (O), relief (R), parent material (P), age (A), and space (N). Most studies have used all or some of these factors to predict and map SOC (Minasny et al., 2013). The considerable progress in proximal sensing apparatus has made the measurement of soil factors easier and faster, leading to more efficient prediction and digital mapping of soil properties. Portable X-ray fluorescence (Kebonye et al., 2021;Mukhopadhyay et al., 2020;Yan et al., 2023), gamma-ray radiation (Wang et al., 2024;Zhang et al., 2020), and visible, near-infrared, and shortwave-infrared (VNIReSWIR) spectroscopy (Ben-Dor et al., 2022;Brodský et al., 2013) are among the most popular proximal sensing data employed to measure soil factors. Although, because these types of data are discontinuous point measurements, different scenarios have been implemented to incorporate them into digital soil mapping. Some studies have implemented point measurements as covariates in co-kriging of the desired soil properties (Kim et al., 2019). Other studies have used them indirectly to calibrate ML or RKbased prediction models and then employed (geo)statistical techniques to interpolate the models' output and obtain the target parameter at unsampled locations (Chakraborty et al., 2017; Kebonye et al., 2021). Ben-Dor &Banin (1995) convolved laboratory spectra of soil samples into Landsat thematic mapper bands (excluding the thermal bands) and found a rough but positive correlation between the proximal and remote sensing data, mainly due to high signal-to-noise ratio of the multispectral sensor. This methodology was employed by de Sousa Mendes et al. (2021) to create a new environmental variable, the Best Synthetic Soil Image, to enhance the predictive power of some soil attributes at three different depths. Point measurements, however, have not directly been used as input variables into digital soil mapping models to obtain the spatial variability of the target variable. This is mainly due to their discrete nature and a large number of containing variables eespecially in the case of soil spectra. This study proposes a novel solution to fill the above-mentioned lacuna by producing rasters via ordinary kriging (OK) of the genetic algorithm (GA)-selected wavelengths of the measured spectra and employing them as input predictor variables alone and coupled with environmental predictor covariates (EPCs), to predict and map SOC in unsampled areas (first objective). To achieve this, three different approaches were used: i) ML approach; to develop models based on four ML algorithms (i.e., partial least squares regression (PLSR), cubist, random forest (RF), and support vector regression (SVR)) on the calibration dataset containing the selected spectra and EPCs (separately and in combination), ii) RK approach; to develop PLSR-K, cubist-K, RF-K, and SVR-K models on the calibration dataset containing the selected spectra and EPCs (separately and in combination), and iii) Kriging approach; to develop OK models on the calibration dataset containing the measured and ML predicted SOC. Results obtained by separated and combined implementation of the spectra and EPCs were then compared to each other eour second Objective. Model inputs are the main sources of uncertainty (Arrouays et al., 2014), especially when these inputs are rasters produced by interpolation of point measurements. Hence, the SGS technique was used to quantify and map the uncertainty associated with the interpolated spectra as input to the SOC mapping models. 2. Materials and methods 2.1. Site description, soil sampling and analysis The study area consisted of three agricultural sites of Nov a Ves nad Popelkou (Nov a Ves), Údrnice (J cin) and Klu cov located in rural regions of the Czech Republic (Fig. 1). The sites mostly produce maize, potatoes, cereals and oilseed rape. According to the World Reference Base (WRB) of soils (IUSS Working Group WRB, 2014), Nov a Ves soil type is mainly Haplic Cambisol formed on the Permian-Carboniferous parent sandstone while Údrnice is comprised of Luvisols, Albic Luvisols, Luvic Chernozems on the Pleistocene loess units. Calcic Chernozem, Luvic Chernozem and Regosols are the dominant soil types in Klu cov developed on the Proterozoic and Paleozoic schist and granodiorite rocks. The sites have similar climate conditions with high rates of actual and potential degradation byerosion. The main terrain features of the sites include side valley, plateau, toe-slope and back-slope. Long-term tillage, no conservation, average 25 cm plough depth and 5e6 course rotation based on the Norfolk system are considered as the main land management practices in the study area (Gholizadeh et al., 2018). Two hundred and forty (240) soil samples were collected from the topsoil (0e20 cm) in June 2021. The position of each sampling point was recorded using a GeoXM (Trimble Inc., Sunnyvale, California, USA) global positioning system (GPS) with 1 m accuracy. All samples were air-dried, ground, sieved (<2 mm), and mixed thoroughly before the total oxidized carbon was measured using the WalkleyeBlack method. V. Khosravi, A. Gholizadeh, R. Kode sov a et al. International Soil and Water Conservation Research xxx (xxxx) xxx 2 2.2. Lab spectroscopy and data pre-processing ASD FieldSpec III Pro FR spectrometer (ASD Inc., Denver, Colorado, USA) was used to record the spectral reflectance of the samples across the VNIReSWIR range (350e2500 nm) by contacting with a high-intensity probe. Soil samples were placed in 5.5 cm diameter petri dishes forming 2 cm layers of soil to avoid beam reflectance from the bottom of the dish (Gholizadeh et al., 2018). Samples were leveled off to obtain a flat surface for maximum light reflection and a high signal to noise ratio. All spectral readings were from the center of the samples. The spectrometer was calibrated using a white Spectralon™(Lab-sphere, North Sutton, New Hampshire, USA) before the first scan and after every ten measurements. Before modeling, the raw spectra (Fig. 2 (a)) were subjected to several pre-processing scenarios including smoothing and noise removal, scatter corrections, derivative transformations and detrending. The most efficient spectral treatment was: i) removing the artificial noise between 350e449 nm and 2451e2500 nm, caused by the device, ii) transformation of reflectance to absorbance via log (1/R) in which R is the reflectance spectra, iii) Savitzky-Golay smoothing with a second-order polynomial fit and 11 smoothing points, and iv) first derivative (FD) transformation to remove the baseline offset (Fig. 2 (b)). Finally, the outliers were removed by applying the H-distance technique on the FD-spectra. Consequently, six samples (four from the Nov a Ves and two from Údrnice) were found not to be consistent with the other observations and removed from further processing. 2.3. Wavelength selection and creating rasters The GA technique aims to select the optimum variables in searching for the best solution from a population of candidates. The higher the solution quality of the regression model, known as fitness, the higher the reproduction probability. This study implemented PLSR on the spectra at GA-selected wavelengths because of Fig. 1. (a) Study sites in the Czech Republic and sampling locations in (b) Klu cov, (c) Nov a Ves nad Popelkou and (d) Údrnice. Fig. 2. The (a) raw and (b) first derivative spectra of the soil samples. V. Khosravi, A. Gholizadeh, R. Kode sov a et al. International Soil and Water Conservation Research xxx (xxxx) xxx 3 its capability to capture the maximum variance and provide the highest correlation between the spectra and the SOC content. ChemometricsWithR package was used for GA-PLSR with random initialization, simple threshold selection, and uniform type crossover (Wehrens, 2011). 500 iterations were done with 30e50 spectra and the default mutation probability of 1%. The same procedure was performed for each site and the most important common wavelengths were selected. The OK method was performed on the FD value at each selected wavelength to obtain a 5 m resolution pixel raster of that wavelength, with the value of every pixel/grid corresponding to 5 m by 5 m square in the study area. It was tried to fit the best theoretical variogram model to the FD values at each selected wavelength and perform an optimum FD estimation of unsampled areas with the lowest error and standard deviation. 2.4. Predictor variables Depending on the approach, various types of predictor variables were used individually or in combination to predict and map the SOC contents. For the ML and RK approaches, the GA-selected measured spectra, EPCs, and their combination were used to construct the prediction models on the calibration dataset. Two groups of EPCs were used in this study (Table A1): i) Sentinel-2Aderived spectral indices, and ii) terrain-based EPCs derived from Advanced Space-borne Thermal Emission and Reflection Radiometer Global Digital Elevation Model (ASTER GDEM) at 30 m spatial resolution (by using SAGA GIS (Conrad et al., 2015),). 2.5. Modeling approaches The classic KennardeStone method, which chooses samples based on a distance measure (Kennard &Stone, 1969), was used to split samples of each site into calibration (75%) and test (25%) datasets. Using 5-fold cross-validation, RK, and ML models were established on the calibration dataset consisting of EPCs and samples' FD value at the GA selected wavelengths, separately and in combination, as predictor variables, and measured SOC as the target variable. The obtained models were evaluated using the test dataset, which included the measured and interpolated spectra, EPCs, and their combination as input, and the measured SOC as the target variables (first and second approaches). Models developed based on the measured spectra were evaluated by applying them separately to the test dataset's measured and interpolated spectra. The models developed on the EPCs were tested on the EPCs data in the pixels, where the test samples were located. The models developed on the combination of the spectra and EPCs were applied to both EPC-measured spectra and EPCs-interpolated spectra in the test samples' locations. This could help to gain a better insight into the efficiency of using the interpolated spectra to predict the desired variable in unsampled locations. Finally, the OK models were developed on the measured SOC and the SOC predicted by the ML models on the calibration dataset (third approach). In the case of ML and RK, the results obtained by using the interpolated spectra and the interpolated spectra-EPCs were compared with those yielded by the measured spectra, EPCs, and the measured spectra-EPCs for the test dataset. This could help to determine how reliable the prediction of SOC in unsampled locations was, by applying the models on the interpolated spectra. Then, the ML and RK models constructed on soil spectra, EPCs, and their combination were compared to each other to understand if their synergy can improve the prediction accuracy. Various types of datasets used for the models’calibration and evaluation in each approach are presented in Table 1. 2.5.1. Machine learning The PLSR, cubist, RF, and SVR algorithms were used to develop SOC models in this study. The PLSR algorithm includes the characteristics of both principal component regression, and stepwise multiple linear regression (Martens &Naes, 1992). Cubist, as a regression tree-based algorithm, is an extension of Quinlan's M5 model tree (Quinlan, 1992). RF, developed by (Breiman, 2001), is a tree-based ensemble learning technique for both classification and regression applications. SVR, as a subcategory of support vector machines, is based on calculating a linear regression function in a multidimensional feature space, in which nonlinear functions are used for mapping the data (Awad &Khanna, 2015). 2.5.2. Kriging As a subclass of spatial estimation methods, kriging is a probabilistic tool that assumes a statistical model for the data. It is based on the semi-variogram (Eq. (1)) of the study area, which provides the spatial variation structure of the target variable in a quantitative form (Webster &Oliver, 2007). g ðhÞ¼ 1 2nX n i¼1 ðZðXiÞZðXiþhÞÞ2(1) where, g (h) is the average of semi-variances between all possible n pairs of samples separated by the lag distance vector of h, and Z(X i ) and Z(X i þh) are the values of a Zvariable at ith pair with distance of h. A reliable variogram model is needed for optimal explanation of the spatial complexity of the area under study and high accuracy kriging results. This can be achieved by separating samples by short to very large lag distances. The most common kriging method is OK, which is a linear unbiased estimator with an error mean equal to zero. The mean value of the target variable over the study area is not needed to be known in OK (Eq. (2)). ZðuÞ¼XNðuÞ j¼1 u jðuÞZuj;s:t:X NðuÞ j¼1 u j¼1 (2) where, Z(u) is the estimated value at point u, N(u) is the number of observed points employed for estimation at point u, u j (u) are the kriging weights with the sum equal to 1 assigned to the observed points. In this study, exponential, spherical, circular, stable, and Gaussian models were tested for experimental semi-variogram calculation based on cross-validation. Furthermore, the SOC content was log-transformed to comply with the normality assumption Table 1 Various types of datasets used for models’calibration and test in each approach. Approach Model Dataset Calibration (75%) Test (25%) ML PLSR Cubist RF SVR Measured spectra - Measured spectra -Interpolated spectra EPCs -EPCs Measure spectra þEPCs -Measure spectra þEPCs -Interpolated spectra þEPCs RK PLSR Cubist RF SVR Measured spectra -Measured spectra -Interpolated spectra EPCs -EPCs Measure spectra þEPCs -Measure spectra þEPCs -Interpolated spectra þEPCs Kriging OK Measured SOC Test samples' locations PLSR predicted SOC Cubist predicted SOC RF predicted SOC SVR predicted SOC V. Khosravi, A. Gholizadeh, R. Kode sov a et al. International Soil and Water Conservation Research xxx (xxxx) xxx 4 of kriging. 2.5.3. Regression kriging (RK) Introduced by Odeh et al. (1995), RK is one of the most popular hybrid spatial techniques for predicting soil properties. In RK, the experimental variogram of regression techniques prediction residuals (on the calibration dataset) is best fitted by an appropriate theoretical variogram. The SOC content at the test samples location, Z(u), is then estimated by the OK of residuals and summing up the kriged values, r(u), to the ML predictions on the test set, P(u) (Eq. (3)). ZðuÞ¼PðuÞþrðuÞ(3) 2.6. Digital soil mapping Further to quantitative evaluation of the calibrated models, the best of those obtained via the ML and RK approaches were applied to the interpolated spectra and EPCs, separately and in combination, to produce the SOC maps of the study area. The best OK model obtained in the third approach was also used to provide a SOC map to compare with the above-mentioned maps. 2.7. Accuracy assessment and uncertainty analysis The mean error (ME), root mean square error (RMSE), coefficient of determination (R 2 ), and Lin's concordance correlation coefficient (LCCC) were used as the evaluation criteria obtained by applying the calibrated models on the test set (Eqs. (4)e(7)). ME ¼Pn i¼1ðoipiÞ n(4) R2¼1Pn i¼1ðoipiÞ2 Pn i¼1ðoi m oÞ2(5) RMSE ¼ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 nX n i¼1 ðoipiÞ2 v u u t(6) LCCC ¼2 rs o s p s 2 oþ s 2 pþ m o m p2(7) where o i is the observed value, p i is the predicted value, n is the number of samples, m o and m p are the means for the observed and predicted values, s 2 o and s 2 p are the corresponding variances, and r is the Pearson's correlation coefficient between the two variables. In general, a robust model has high R 2 and LCCC and low ME and RMSE values. Analysis of uncertainty was performed using SGS by simulating 100 realizations (truths) of each selected wavelength raster into 5 m grids across each site and using the simulated realizations to feed 100 different uncertain ML prediction models. The standard deviation of the 100 predicted SOC values were determined and considered as each ML-models’uncertainty. The uncertainty associated with RK approach models were also obtained in a very similar way with 100 more realizations simulated on the prediction residuals of each of the 100 realizations of the ML models. The uncertainty distribution, which was equal in all samples, was not involved in any of the modeling processes of this study. A schematic view of the methodology is shown in Fig. 3. 3. Results 3.1. Statistical description Summary statistics of the SOC content of the collected samples (excluding outliers) are presented in Table 2. It can be seen that the SOC contents were low, below 2.59% as the maximum value. The values ranged between 0.88 and 2.59%, 0.50e1.87%, and 0.77e1.61% with the mean values of 1.38%, 1.04%, and 1.11% for Nov a Ves, Údrnice, and Klu cov, respectively. SOC skewed positively in all sites, with the highest skewness value of 1.31 for Nov a Ves, indicating the largest asymmetry in distribution of the target variable. The coefficient of variation (CV) indicated low to moderate variabilities for all sites (0.18 <CV <0.24), with the highest value for Údrnice, which can mainly be due to that it lies in a highland with complex relief that is naturally more heterogeneous. Samples of no site followed the normal distribution; hence logarithmic transformation was applied to the samples to guarantee the normality, where needed. 3.2. GA-selected wavelengths For each site, FD rasters at 30 top important GA-selected wavelengths were used as the predictor variables for our further processing. The OK-interpolated rasters of the selected spectral variables of Nov a Ves can be seen in Fig. 4. All rasters were created using spectra of the calibration dataset. The same resulting rasters for Údrnice and Klu cov can be seen in Fig. A1 and Fig. A2, in the supplementary materials. 3.3. Performance of models on interpolated spectra The results of applying the ML and RK models, on interpolated spectra (as test dataset) are presented in Tables 3 and 4. As shown in Table 3, most of the models developed under both ML and RK approaches, yielded relatively weak results when evaluated on the interpolated spectra. Considering the RK approach, the best result on the interpolated spectra was obtained by SVR-K with ME ¼0.03%, RMSE ¼0.12%, R 2 ¼0.69, and LCCC ¼0.76, for Údrnice, ME ¼0.05%, RMSE ¼0.13%, R 2 ¼0.62, and LCCC ¼0.72 for Klu cov, and ME ¼0.09%, RMSE ¼0.19%, R 2 ¼0.42, and LCCC ¼0.61 for Nov a Ves. Applying models on the measured spectra, yielded outputs with the same order, but higher accuracy ei.e., SVR-K with ME ¼0%, RMSE ¼0.09%, R 2 ¼0.77, and LCCC ¼0.85 for Údrnice, ME ¼0.01%, RMSE ¼0.10%, R 2 ¼0.75, and LCCC ¼0.78 for Klu cov, and ME ¼0.08%, RMSE ¼0.16%, R 2 ¼0.59, and LCCC ¼0.63 for Nov a Ves. Similarly, the performance of all models under the ML approach dropped, and SVR efor instanceeshowed lower prediction accuracy on the interpolated spectra (ME ¼0.03%, RMSE ¼0.15%, R 2 ¼0.65, and LCCC ¼0.73 for Údrnice, ME ¼0.06%, RMSE ¼0.17%, R 2 ¼0.53, and LCCC ¼0.64 for Klu cov, and ME ¼0.1%, RMSE ¼0.22%, R 2 ¼0.40, and LCCC ¼0.51 for Nov a Ves) than the measured spectra (ME ¼0%, RMSE ¼0.11%, R 2 ¼0.73, and LCCC ¼0.82 for Údrnice, ME ¼0.03%, RMSE ¼0.12%, R 2 ¼0.66, and LCCC ¼0.69 for Klu cov, and ME ¼0.10%, RMSE ¼0.19%, R 2 ¼0.56 and LCCC ¼0.61 for Nov a Ves). Similar results were obtained using the combination of spectra and EPCs as the test dataset (Table 4). Accordingly, the models under RK provided higher SOC prediction performance when applied to the combination of measured spectra and EPCs, compared to the combination of interpolated spectra and EPCs. Alike for the ML approach, better prediction results were obtained using the combination of measured spectra and EPCs than the interpolated spectra and EPCs. In addition, the results obtained on the interpolated spectra V. Khosravi, A. Gholizadeh, R. Kode sov a et al. International Soil and Water Conservation Research xxx (xxxx) xxx 5 separately and combined with EPCs were compared with the results of the third approach (presented in Table 5 and Fig. A3), the SOC estimations (at the test dataset) obtained by the OK models calibrated on the measured and ML-predicted SOC contents (at the calibration dataset). Accordingly, the results of OK on the measured SOC (ME ¼0.05, RMSE ¼0.12%, R 2 ¼0.71, and LCCC ¼0.79 for Údrnice, ME ¼0.03%, RMSE ¼0.14%, R 2 ¼0.65, and LCCC ¼0.77 for Klu cov, ME ¼0.08%, RMSE ¼0.21%, R 2 ¼0.42, and LCCC ¼0.60 for Nov a Ves) were better than those obtained by the ML and RK approaches on the interpolated spectra, alone. That was while the SOC predicted by the ML and RK approaches on the interpolated spectra were more accurate than all the SOC values estimated by the OK of ML-predicted SOC. In other words, using the interpolated spectra as ML and RK models input, yields more accurate SOC than interpolation of the RK and ML models’output SOC. It is also noteworthy to mention that the RK and ML approaches performed on combination of the interpolated spectra and EPCs had superior performance than all models under the third approach. 3.4. The effect of combining spectra and EPCs on developed models Comparing results obtained using the separated and combined samples spectra and EPCs (Tables 3 and 4), can highlight the better performance of models calibrated on the measured spectra and EPCs than those calibrated on any of them separately. For Nov a Ves, as an example, the best performance between all approaches and methods was attained using SVR-K on the combination of measured spectra and EPCs (R 2 ¼0.65, RMSE ¼0.15%), followed by SVR-ML on the same dataset (R 2 ¼0.62, RMSE ¼0.17%), SVR-K on the measured spectra (R 2 ¼0.59, RMSE ¼0.16%), and SVR-ML on the measured spectra (R 2 ¼0.56, RMSE ¼0.19%). In addition, the combination of interpolated spectra and EPCs yielded better predictions than using them individually (Tables 3 and 4). ML and RK on the interpolated spectra combined with the EPCs also provided higher accuracy (R 2 ¼0.51 and RMSE ¼0.17% for SVR-K and R 2 ¼0.49 and RMSE ¼0.19% for SVR-ML in Nov a Ves) than OK on the measured SOC (R 2 ¼0.42 and RMSE ¼0.21%). This further proved that data types had the highest impact on SOC prediction, in this study. By comparing the results obtained on all sites and types of datasets, following decreasing performance order can be observed: the ML and RK approaches calibrated and tested on the measured spectra combined with the EPCs >the ML and RK approaches calibrated and tested on the measured spectra >the ML and RK approaches calibrated on the measured spectra combined with the EPCs, and tested on the interpolated spectra combined with the EPCs >the OK-calibrated on the measured SOC >the ML and RK approaches calibrated and tested on the EPCs >the ML and RK approaches calibrated on the measured spectra and tested on the interpolated point spectra >the OK-calibrated on the ML-predicted SOC. 3.5. General comparison between approaches and algorithms Regardless of the sites and datasets, SVR presented the best accuracies among all algorithms, followed by RF. Between PLSR and Cubic, no method showed a clear superiority. However, by accepting that RMSE is a more important criterion than R 2 (Willmott, 1981), PLSR was more successful than cubist. The poorest SOC prediction was yielded using OK applied on the cubistpredicted SOC values with ME ¼0.14%, RMSE ¼0.25%, R 2 ¼0.21, and LCCC ¼0.29. 3.6. Spatial distribution maps Our main motivation for spectra interpolation was to employ it for digital mapping of SOC. As obtained previously (Tables 3 and 4), SVR eunder the RK (SVR-K) and ML (SVR) approachesecalibrated Fig. 3. Methodology flowchart of this study. Table 2 Summary statistics of SOC content (%) of the samples gathered from the study area. Site n Min Max Mean StD CV Skewness Nov a Ves 76 0.88 2.59 1.38 0.29 0.21 1.31 Údrnice 78 0.50 1.87 1.04 0.25 0.24 0.94 Klu cov 80 0.77 1.61 1.11 0.20 0.18 0.33 n: number of samples, Min: Minimum, Max: Maximum, StD: Standard deviation, CV: Coefficient of variation. V. Khosravi, A. Gholizadeh, R. Kode sov a et al. International Soil and Water Conservation Research xxx (xxxx) xxx 6 on the separated and combined spectra and EPCs, provided the highest accuracies for SOC prediction. These models were used to create the 5 m resolution SOC spatial distribution maps of all sites, along with the distribution map obtained by OK of the measured SOC (Fig. 5). As can be seen in all five maps of Nov a Ves, the lowest SOC contents are in the northern and southern parts being increased gradually toward the middle of the site, showing higher SOC contents of over 1.3%. For Údrnice, the southern and middle parts (upper portion) show relatively higher SOC contents in all the maps. The low SOC is also evident in the northern and middle parts (lower portion) of the site with contents lower than about 1%. In the case of Klu cov, the southern and middle parts show the highest SOC values, while the northern and eastern parts have relatively lower SOC contents of under about 1.1%. In general, comparing the approaches and datasets used, the spatial pattern of SOC was almost similar, though, some differences were observed in the range of predicted SOC. As an example, for all sites, the range of the predicted SOC using SVR (Fig. 5(b)) and SVR-K (Fig. 5d) on the interpolated spectra combined with EPCs, were more similar to the actual SOC limits, than the range of the predicted SOC using SVR (Fig. 5(a)) and SVR-K (Fig. 5(c)) on the interpolated spectra. These results were compatible with those reported in Tables 3 and 4. Finally, the maps obtained using SVR-K and OK (Fig. 5 (c), (d), (e)) are characterized by being smoother than the maps obtained using SVR (Fig. 5(a) and (b)), which might be due to the inherent smoothing effect of the OK method. 3.7. Uncertainty The average uncertainty associated with the SOC prediction models applied on the interpolated spectra is shown in Table 6. The Fig. 4. Rasters obtained by OK-based interpolation of the samples FD spectra at the selected wavelengths (Nov a Ves). V. Khosravi, A. Gholizadeh, R. Kode sov a et al. International Soil and Water Conservation Research xxx (xxxx) xxx 7 average uncertainty resulting from applying the SGS on the measured SOC is also added for comparison. As evidenced, the ML approach mostly showed lower uncertainty than RK, regardless of the data types, algorithms, and sites (Table 6). Compared to the kriging approach, the uncertainty obtained by SVR-ML was lower for both types of datasets, while, the uncertainty of the other algorithms, was more than the measured SOC. Models developed on the combination of spectra and EPCs released lower uncertainty than those solely calibrated on the spectra. Comparing the algorithms, SVR caused the lowest uncertainty in both approaches and datasets, followed by RF. No significant difference was observed between PLSR and cubist models’resulting uncertainties. Considering the sites, Klu cov had lower uncertainty than Nov a Ves and Údrnice. The uncertainty maps of SOC obtained by applying the best ML and RK models on the interpolated spectra and the combination of Table 3 Results of ML and RK approaches applied on measured and interpolated point spectra (test samples). Site Approach Model Measured point spectra Interpolated point spectra ME RMSE R 2 LCCC ME RMSE R 2 LCCC Nov a Ves ML PLSR 0.07 0.20 0.46 0.54 0.09 0.24 0.22 0.33 Cubist 0.12 0.22 0.43 0.49 0.14 0.23 0.26 0.37 RF 0.08 0.22 0.56 0.59 0.11 0.23 0.34 0.44 SVR 0.10 0.19 0.56 0.61 0.10 0.22 0.40 0.51 RK PLSR 0.09 0.21 0.45 0.58 0.13 0.24 0.25 0.54 Cubist 0.10 0.22 0.49 0.53 0.14 0.22 0.29 0.47 RF 0.06 0.19 0.58 0.62 0.09 0.22 0.34 0.59 SVR 0.08 0.16 0.59 0.63 0.09 0.19 0.42 0.61 Údrnice ML PLSR 0.06 0.12 0.77 0.84 0.08 0.15 0.35 0.53 Cubist 0.03 0.13 0.65 0.79 0.10 0.18 0.39 0.44 RF 0.04 0.14 0.59 0.68 0.06 0.16 0.47 0.56 SVR 0 0.11 0.73 0.82 0.03 0.15 0.65 0.73 RK PLSR 0.02 0.11 0.78 0.86 0.05 0.14 0.42 0.60 Cubist 0.04 0.14 0.67 0.81 0.09 0.16 0.49 0.57 RF 0.01 0.12 0.70 0.73 0.07 0.13 0.59 0.68 SVR 0 0.09 0.77 0.85 0.03 0.12 0.69 0.76 Klu cov ML PLSR 0 0.15 0.43 0.71 0.04 0.18 0.41 0.59 Cubist 0.05 0.12 0.57 0.56 0.09 0.18 0.51 0.48 RF 0.02 0.13 0.72 0.75 0.08 0.16 0.52 0.58 SVR 0.03 0.12 0.66 0.69 0.06 0.17 0.53 0.64 RK PLSR 0.05 0.12 0.59 0.73 0.04 0.16 0.55 0.62 Cubist 0.03 0.14 0.46 0.59 0.07 0.15 0.40 0.51 RF 0.05 0.13 0.74 0.78 0.08 0.14 0.52 0.65 SVR 0.01 0.10 0.75 0.78 0.05 0.13 0.62 0.72 Table 4 Results of ML and RK approaches applied on environmental predictor covariates (EPCs), measured point spectra and EPCs and interpolated point spectra and EPCs (test samples). Site Approach Model EPCs Measured point spectra & EPCs Interpolated spectra & EPCs RMSE R 2 RMSE R 2 RMSE R 2 Nov a Ves ML PLSR 0.23 0.24 0.18 0.44 0.22 0.36 Cubist 0.22 0.28 0.22 0.46 0.22 0.38 RF 0.23 0.35 0.20 0.59 0.21 0.42 SVR 0.20 0.43 0.17 0.62 0.19 0.49 RK PLSR 0.23 0.28 0.18 0.51 0.23 0.39 Cubist 0.21 0.32 0.19 0.47 0.21 0.41 RF 0.22 0.37 0.19 0.62 0.21 0.45 SVR 0.18 0.45 0.15 0.65 0.17 0.51 Údrnice ML PLSR 0.14 0.41 0.12 0.79 0.14 0.45 Cubist 0.18 0.42 0.12 0.68 0.16 0.43 RF 0.16 0.49 0.11 0.63 0.15 0.53 SVR 0.14 0.68 0.10 0.79 0.12 0.74 RK PLSR 0.12 0.44 0.11 0.81 0.12 0.49 Cubist 0.16 0.43 0.11 0.69 0.14 0.58 RF 0.13 0.54 0.10 0.71 0.12 0.62 SVR 0.11 0.71 0.09 0.87 0.11 0.75 Klu cov ML PLSR 0.17 0.52 0.12 0.59 0.15 0.55 Cubist 0.17 0.41 0.14 0.47 0.16 0.43 RF 0.16 0.56 0.13 0.75 0.14 0.58 SVR 0.15 0.59 0.12 0.68 0.13 0.61 RK PLSR 0.16 0.44 0.12 0.61 0.14 0.59 Cubist 0.15 0.59 0.13 0.50 0.15 0.46 RF 0.13 0.61 0.12 0.75 0.13 0.66 SVR 0.11 0.67 0.10 0.78 0.11 0.71 Table 5 Results of OK applied on measured SOC and ML predicted SOC (test samples). Site Metrics Measured SOC PLSR SOC Cubist SOC RF SOC SVR SOC Nov a Ves ME 0.08 0.16 0.14 0.13 0.11 RMSE 0.21 0.24 0.25 0.24 0.22 R 2 0.42 0.20 0.21 0.30 0.37 LCCC 0.60 0.32 0.29 0.41 0.46 Údrnice ME 0.05 0.10 0.11 0.08 0.05 RMSE 0.12 0.19 0.21 0.17 0.16 R 2 0.71 0.33 0.31 0.51 0.64 LCCC 0.79 0.50 0.42 0.63 0.67 Klu cov ME 0.03 0.07 0.1 0.09 0.06 RMSE 0.14 0.21 0.18 0.18 0.16 R 2 0.65 0.39 0.40 0.46 0.52 LCCC 0.77 0.47 0.43 0.49 0.60 V. Khosravi, A. Gholizadeh, R. Kode sov a et al. International Soil and Water Conservation Research xxx (xxxx) xxx 8 interpolated spectra and EPCs are presented in Fig. 6. The measured SOC uncertainty map (obtained by SGS under the kriging approach) is also added for comparison. The uncertainty distribution patterns of the sites were almost similar for all ML, RK, and Kriging approaches. However, the uncertainty ranges were variable among different approaches, as it is also evident in Table 6. For instance, for Nov a Ves, the range of uncertainty was higher for RK (Fig. 6(c) and (d): 0.17%e0.35%) than ML (Fig. 6(a) and b: 0.07%e0.18%) and OK (Fig. 6(e): 0.11%e0.18%). Likewise, for the other sites, the uncertainty ranges of the maps created using the RK approach (Fig. 6(c) and d: 0.12%e0.32% for Údrnice and 0.11%e0.28% for Klu cov) were higher than those obtained from ML and OK. Compatible with the results reported in Table 6, the lowest uncertainty variations were observed for SVRML (applied on the interpolated spectra combined with EPCs) followed by OK. The uncertainty maps illustrated in Fig. 6,presented higher spatial details than the prediction maps, which is mainly due to the higher variabilities in prediction confidence and SGS-related details of the uncertainty maps. The inherent smoothing effect of OK, which eliminates the fine-scale SOC variabilities in prediction maps can be considered as the other important reason. 4. Discussion 4.1. Interpolated spectra as input variables Focusing on the main objective of this study, applying the developed models on the interpolated spectra yielded relatively acceptable results. But how logical is incorporating rasters of the Fig. 5. The spatial distribution maps of SOC produced by (a) SVR-ML applied on the interpolated spectra (ML approach), (b) SVR-ML applied on the interpolated spectra combined with EPCs (ML approach), (c) SVR-K applied on the interpolated spectra (RK approach), (d) SVR-K applied on the interpolated spectra combined with EPCs (RK approach), and (e) OK of measured SOC (kriging approach). Table 6 Average uncertainty (%) of SOC prediction models under different approaches. Approach Model Lab spectra Lab spectra &EPCs Measured SOC Nov a Ves Údrnice Klu cov Nov a Ves Údrnice Klu cov Nov a Ves Údrnice Klu cov RK PLSR 14.13 15.26 11.02 12.11 10.84 9.36 eee Cubist 15.73 11.29 10.81 11.89 9.97 10.65 eee RF 14.19 10.86 10.76 10.25 9.36 8.29 eee SVR 13.51 10.91 9.35 11.48 8.91 7.17 eee ML PLSR 12.16 13.20 10.89 11.07 10.15 8.97 eee Cubist 11.62 12. 77 11.72 10.45 9.21 10.53 eee RF 11.87 9.74 8.92 10.23 8.73 6.97 eee SVR 10.76 7.81 7.65 9.70 6.62 5.49 eee Kriging SGS eeeeee10.18 7.07 6.91 V. Khosravi, A. Gholizadeh, R. Kode sov a et al. International Soil and Water Conservation Research xxx (xxxx) xxx 9