Xianglin Zhang

Songchao Chen
email: chensongchao@zju.edu.cn

Jie Xue

Nan Wang

Yi Xiao

Qianqian Chen

Yongsheng Hong

Yin Zhou

Hongfen Teng

Bifeng Hu

Zhiqing Zhuo

Wenjun Ji

Yuanfang Huang

Yuxuan Gou

Anne C Richer-De-Forges

Dominique Arrouays

Zhou Shi

Morgan Cristine

Improving model parsimony and accuracy by modified greedy feature selection in digital soil mapping

Keywords:

In the context of increasing soil degradation worldwide, spatially explicit soil information is urgently needed to support decision-making for sustaining limited soil resources. Digital soil mapping (DSM) has been proven as an efficient way to deliver soil information from local to global scales. The number of environmental covariates used for DSM has rapidly increased due to the growing volume of remote sensing data, therefore variable selection is necessary to deal with multicollinearity and improve model parsimony. Compared with Boruta, recursive feature elimination (RFE), and variance inflation factor (VIF) analysis, we proposed the use of modified greedy feature selection (MGFS), for DSM regression. For this purpose, using quantile regression forest, 402 soil samples and 392 environmental covariates were used to map the spatial distribution of soil organic carbon density (SOCD) in Northeast and North China. The result showed that MGFS selected the most parsimonious model with only 9 covariates (e.g., brightness index, mean annual temperature), much lower than RFE (22 covariates), VIF (30 covariates), and Boruta (76 covariates). The repeated validation (50 times) showed that the MGFS derived model performed better (R 2 of 0.60, LCCC of 0.74, RMSE of 13.80 t ha -1 ) than these using full covariates, Boruta, RFE and VIF (R 2 of 0.48-0.57, LCCC of 0.64-0.72, RMSE of 14.24-15.79 t ha -1 ). Despite the similar performance of the uncertainty estimate (PICP), the model using MGFS and RFE had the lowest global uncertainty (0.86) as indicated by the uncertainty index. In addition, MGFS had the best computation efficiency when considering the steps of variable selection and map prediction. Given these advantages over Boruta, RFE and VIF, MGFS has a high potential in fine-resolution soil mapping practices, especially for these studies at a broad scale involving heavy computation on millions or billions of pixels.

Introduction

Soil is one of the Earth's most essential and finite resources. It enables life on Earth by delivering crucial ecosystem services, including the provision of food, fibre and fuel, water purification, contaminant reduction, nutrient cycling, carbon sequestration, climate and flood regulation, and biodiversity conservation [START_REF] Mcbratney | The dimensions of soil security[END_REF][START_REF] Adhikari | Linking soils to ecosystem services-A global review[END_REF][START_REF] Baveye | Soil "ecosystem" services and natural capital: critical appraisal of research on uncertain ground[END_REF][START_REF] Pereira | Soil ecosystem services, sustainability, valuation and management[END_REF]. Under the tremendous pressure of population growth, economic development and climate change, global soils are continuously degraded (e. g., erosion, salinization, fertility decline). [START_REF] Fao | Status of the world's soil resources: main report[END_REF] reported that 33 % of the Earth's soils are already degraded, and over 90 % could become degraded by 2050. For sustaining soil resources for the next generation, there is an urgent demand to improve management practices which require spatially explicit soil information to support decision-making. However, many conventional soil maps do not provide relevant information for the current global challenges, and most of the data for creating the soil maps is outdated [START_REF] Sanchez | Digital soil map of the world[END_REF][START_REF] Arrouays | GlobalSoilMap: Toward a fine-resolution global grid of soil properties[END_REF][START_REF] Arrouays | Digital soil mapping across the globe[END_REF]. A digital revolution has taken place in soil mapping in the last 20 years, which resulted in the concept of digital soil mapping (DSM, McBratney et al., 2003).

DSM has emerged as a fast-growing sub-discipline of soil science by integrating soil surveys, geostatistics, geographic information system, remote sensing and machine learning [START_REF] Minasny | Digital soil mapping: A brief history and some lessons[END_REF]. DSM has been widely used in delivering fine-resolution maps of soil information across scales with a focus switching from "primary" soil properties (e.g., soil organic carbon, particle size fractions, pH, soil depth) to "secondary" soil properties (e.g., available water capacity, carbon sequestration potential, carbon vulnerability) in the last decade [START_REF] Padarian | Predicting and mapping the soil available water capacity of Australian wheatbelt[END_REF][START_REF] Hengl | SoilGrids250m: Global gridded soil information based on machine learning[END_REF][START_REF] Chen | Fine resolution map of top-and subsoil carbon sequestration potential in France[END_REF][START_REF] Román Dobarco | Uncertainty assessment of GlobalSoilMap soil available water capacity products: A French case study[END_REF][START_REF] Viscarra Rossel | Continental-scale soil carbon composition and vulnerability modulated by regional environmental controls[END_REF][START_REF] Liu | Mapping high resolution National Soil Information Grids of China[END_REF].

The number of environmental covariates used for DSM has rapidly increased due to the growing volume of remote sensing data (Chen et al., 2022a). Therefore, variable selection is increasingly used before fitting the final predictive model. It has several advantages, such as (1) calibrating the predictive model faster;

(2) reducing model complexity; (3) increasing model performance; (4) avoiding multicollinearity; (5) producing the map faster [START_REF] Wadoux | Machine learning for digital soil mapping: Applications, challenges and suggested solutions[END_REF]. Currently, there are mainly-two strategies for variable selection in DSM [START_REF] Wadoux | Machine learning for digital soil mapping: Applications, challenges and suggested solutions[END_REF]: (1) covariate reduction as a pre-processing step (before calibrating a machine learning model), such as selecting the most relevant covariates by Pearson's correlation coefficient between soil properties and covariates, discarding the covariates which are highly correlated with other covariates, or keeping the first several components using principle component analysis [START_REF] Mosleh | The effectiveness of digital soil mapping to predict soil properties over low-relief areas[END_REF][START_REF] Hamzehpour | Exploring the driving forces and digital mapping of soil organic carbon using remote sensing and soil texture[END_REF][START_REF] Poggio | SoilGrids 2.0: producing soil information for the globe with quantified spatial uncertainty[END_REF][START_REF] Taghizadeh-Mehrjardi | Enhancing the accuracy of machine learning models using the super learner technique in digital soil mapping[END_REF]; (2) wrapper methods which rely on inference made by calibrating a machine learning model to assess covariate importance [START_REF] Xiong | Holistic environmental soil-landscape modeling of soil organic carbon[END_REF][START_REF] Brungard | Machine learning for predicting soil classes in three semi-arid landscapes[END_REF][START_REF] Nussbaum | Evaluation of digital soil mapping approaches with large sets of environmental covariates[END_REF][START_REF] Amiri | Assessment of the importance of gully erosion effective factors using Boruta algorithm and its spatial modeling and mapping using three machine learning algorithms[END_REF][START_REF] Keskin | Digital mapping of soil carbon fractions with machine learning[END_REF][START_REF] Poggio | SoilGrids 2.0: producing soil information for the globe with quantified spatial uncertainty[END_REF]. In addition to the limitation of not being applicable to categorical variables, most of the first strategy only accounts for the linear relationship between soil properties and covariates or even omits their correlations, so it can potentially neglect these covariates non-linearly correlated to the soil properties and thus decrease model performance [START_REF] Camera | A high resolution map of soil types and physical properties for Cyprus: A digital soil mapping optimization[END_REF][START_REF] Zeraatpisheh | Digital mapping of soil properties using multiple machine learning in a semi-arid region, central Iran[END_REF]. In addition, though machine learning models can be affected by multicollinearity [START_REF] Strobl | Conditional variable importance for random forests[END_REF][START_REF] Drobnič | On the interpretability of machine learning models and experimental feature selection in case of multicollinear data[END_REF], they are still more robust than traditional multiple linear regression [START_REF] Dormann | Collinearity: a review of methods to deal with it and a simulation study evaluating their performance[END_REF]; therefore, the second strategy is more appropriate for machine learning based DSM studies. Among these wrapper methods for variable selection, the most popular ones are Boruta [START_REF] Xiong | Holistic environmental soil-landscape modeling of soil organic carbon[END_REF][START_REF] Amiri | Assessment of the importance of gully erosion effective factors using Boruta algorithm and its spatial modeling and mapping using three machine learning algorithms[END_REF][START_REF] Keskin | Digital mapping of soil carbon fractions with machine learning[END_REF][START_REF] Rasaei | Rescue and renewal of legacy soil resource inventories in Iran as an input to digital soil mapping[END_REF][START_REF] Xu | Digital soil mapping of soil total nitrogen based on Landsat 8, Sentinel 2, and WorldView-2 images in smallholder farms in Yellow River Basin, China[END_REF][START_REF] Zeraatpisheh | Improving the spatial prediction of soil organic carbon X[END_REF] and recursive feature elimination (RFE, [START_REF] Brungard | Machine learning for predicting soil classes in three semi-arid landscapes[END_REF][START_REF] Nussbaum | Evaluation of digital soil mapping approaches with large sets of environmental covariates[END_REF][START_REF] Gomes | Modelling and mapping soil organic carbon stocks in Brazil[END_REF][START_REF] Chen | Digital mapping of the soil thickness of loess deposits over a calcareous bedrock in central France[END_REF][START_REF] Poggio | SoilGrids 2.0: producing soil information for the globe with quantified spatial uncertainty[END_REF][START_REF] Yang | The effectiveness of digital soil mapping with temporal variables in modeling soil organic carbon changes[END_REF]. Previous studies commonly chose either Boruta or RFE in variable selection, and only a few studies compared RFE and Boruta in DSM practices (Chen et al., 2022b;[START_REF] Luo | Regional soil organic matter mapping models based on the optimal time window, feature selection algorithm and Google Earth Engine[END_REF]. In addition, all these relevant studies only focused on the impact of variable selection on the model performance while ignoring its effect on the uncertainty estimate. Therefore, it remains unclear which variable selection method is superior in model performance and uncertainty estimates. From our previous experience with Boruta and RFE, we also found that they can potentially ignore some useful variables so as to decrease model performance [START_REF] Xiao | Improving pedotransfer functions for predicting soil mineral associated organic carbon by ensemble machine learning[END_REF]. Therefore, based on the recently proposed greedy feature selection (GFS, [START_REF] Drobnič | On the interpretability of machine learning models and experimental feature selection in case of multicollinear data[END_REF], a modified variable selection algorithm is expected to solve this problem for mapping soil properties.

Addressing the current knowledge gap and limitation, the objectives of this study are twofold: (1) apply a modified GFS algorithm to DSM regression problem; (2) compare the modified GFS algorithm to commonly used variable selection methods on model performance, uncertainty estimate and computation efficiency in soil organic carbon density (SOCD) mapping.

Materials and methods

Soil data

The study area is located in Northeast and North China, covering a total area of 56.30 × 10 4 km 2 in which 37.67 × 10 4 km 2 are used as cropland (Fig. 1). It has a temperate continental monsoon climate with mean annual precipitation (MAP) between 400 and 1200 mm and mean annual temperature (MAT) between -1.1 and 16 • C. According to the Genetic Soil Classification of China, the main soil types are Black (Phaeozems), Chernozem and Brown Earth (Luvisols) in Northeast Plain, and Fluvo-aquic (Eutric Cambisols), Shajiang Black (Vertisols), and Cinnamon (Eutric Luvisols) in Huai-Hai Plain [START_REF] Zhuo | Soil organic carbon storage, distribution, and influencing factors at different depths in the dryland farming regions of Northeast and North China[END_REF]. The single cropping system (i.e., spring corn, spring wheat, soybean) dominates in the Northeast Plain, and the double cropping system (winter wheat and summer corn) is the main farming system in the Huai-Hai Plain. Accounting for around 30 % of the national grain production, the study area is a crucial agricultural region in China.

A total of 402 sampling sites were selected in the cropland by stratified random sampling representing the major soil types, clay content and agricultural system. For each sampling site, three soil cores were collected at four depth intervals (0-10, 10-20, 20-30, and 30-40 cm) using an undisturbed soil sampler (diameter of 5.1 cm and height of 10 cm) between April and May in 2017. All the composite soil samples were air-dried and sieved to < 2 mm, and SOC content (g kg -1 ) was determined by the dichromate oxidation-external heating method [START_REF] Bao | Soil agro-chemistrical analysis[END_REF]. Soil bulk density (BD, g cm -3 ) was determined by the mean of three replicates of undisturbed soil cores (100 cm 3 in volume). More details about sampling design and laboratory analysis can be found in [START_REF] Zhuo | Soil organic carbon storage, distribution, and influencing factors at different depths in the dryland farming regions of Northeast and North China[END_REF]. In this study, we focused on the topsoil (0-30 cm) by integrating the soil information from the first three depth intervals (0-10, 10-20, 20-30 cm) using the weighted average method. SOCD (t ha -1 ) was calculated by the equation below:

SOCD = SOC × BD × Depth × 10 (1)
where Depth is the soil depth (cm) which is 30 cm. In this study, coarse fragments did not present in the topsoil (0-30 cm) of cropland, so it was not considered in this equation.

Environmental covariates

Following the Scorpan framework [START_REF] Mcbratney | On digital soil mapping[END_REF], 392 environmental covariates relevant to SOCD were investigated and listed in Table 1. The environmental covariates were derived from multiple sources: Landsat 5 Thematic Mapper (TM), Landsat 7 Enhanced Thematic Mapper (ETM + ), Landsat 8 Operational Land Imager (OLI), Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM), products from Moderate Resolution Imaging Spectroradiometer (MODIS), and legacy digital soil maps.

All available Landsat 5, 7 and 8 surface reflectance images at 30 m resolution covering the entire study area in the Collection 2 Tier1 Level 2 database were collected from 2003 to 2017 (https://www.usgs.gov/). These images received radiometric, geometric and atmospheric corrections with Landsat Ecosystem Disturbance Adaptive Processing System (Landsat 5 TM and Landsat 7 ETM+) and Land Surface Reflectance Code (Landsat 8 OLI). To ensure the image quality, we filtered the image collection with the criteria of cloud cover < 20 %, root mean square error (RMSE) of geometric residuals measured on the ground control points < 10 m, and the best image quality level. Afterwards, the CFmask algorithm was used to identify and remove the dilated cloud, cirrus, cloud and cloud shadow [START_REF] Zhu | Improvement and expansion of the Fmask algorithm: cloud, cloud shadow, and snow detection for Landsats 4-7, 8, and Sentinel 2 images[END_REF][START_REF] Foga | Cloud detection algorithm comparison and validation for operational Landsat data products[END_REF].

Considering the difference in the spectral characteristics of Landsat TM, ETM+ and OLI, we harmonized blue, green, red, near-infrared, shortwave infrared 1, and shortwave infrared 2 in Landsat TM, ETM+ to OLI with the ordinary least square regression coefficient [START_REF] Roy | Characterization of Landsat-7 to Landsat-8 reflective wavelength and normalized difference vegetation index continuity[END_REF]. The products of MODIS from 2003 to 2017 were also used as environmental covariates.

The Geospatial Soil Sensing System (GEOS3) algorithm was used to detect the bare soil pixels in all available Landsat images from 2003 to 2017 [START_REF] Demattê | Geospatial Soil Sensing System (GEOS3): A powerful data mining procedure to retrieve soil spectral reflectance from satellite images[END_REF], and then a median composite image of all bare soil pixels was used to calculate 11 soil indices such as Brightness Index (BI), Hue Index (HI) and Saturation Index (SI). Additionally, legacy digital soil maps such as soil erosion (SE) and soil particle size fractions were included [START_REF] Hengl | SoilGrids250m: Global gridded soil information based on machine learning[END_REF][START_REF] Teng | Modelling and mapping soil erosion potential in China[END_REF].

The climate covariates mainly originated from WorldClim2, including 7 climatological and 16 bioclimatic (BIO) variables at 1 km resolution for the period of 1970-2000 [START_REF] Fick | WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas[END_REF]. Moreover, the average of Day-time and Night-time Land Surface Temperature (LSTD, LSTN) at 1 km resolution were calculated using MOD11A1 [START_REF] Wan | MODIS/Terra Land Surface Temperature/Emissivity Daily L3 Global 1km SIN Grid V061[END_REF]. As for Normalized Difference Snow Index (NDSI), we used Landsat (2003Landsat ( -2017) ) to calculate the intra-annual and intra-quarterly mean and standard deviation (10 variables), mean of the intra-annual and intra-quarterly maximum and minimum (10 variables), and averaged intra-quarterly standard deviation (1 variable), resulting in 21 variables to better represent its temporal change (Zhou et al., 2019a). As for Potential Evapotranspiration (PET) and Evapotranspiration (ET), we used MOD16A2 at 500 m resolution to calculate the mean and standard deviation of the intra-annual and intra-quarterly aggregation (10 variables), and averaged standard deviation of intra-quarterly aggregation in 2003-2017 (1 variable) [START_REF] Running | MOD16A2 MODIS/Terra Net Evapotranspiration 8-Day L4 Global 500m SIN Grid V006[END_REF].

A total of 253 covariates related to organisms were calculated by Landsat and MODIS products [START_REF] Loveland | Landsat: Building a strong future[END_REF][START_REF] Myneni | MODIS/Terra+Aqua Leaf Area Index/FPAR 4-Day L4 Global 500m SIN Grid V061[END_REF][START_REF] Running | MOD17A2H MODIS/Terra Gross Primary Productivity 8-Day L4 Global 500m SIN Grid V006[END_REF]. Leaf Area Index (LAI), Enhanced Vegetation Index (EVI), Normalized Difference Vegetation Index (NDVI) and eight other indices were calculated in the same way as NDSI (21 variables for each). The calculation of Gross Primary Production (GPP) and Net Primary Production (NPP) from MOD17A2H was consistent with PET and ET (11 variables for each) [START_REF] Running | MOD17A2H MODIS/Terra Gross Primary Productivity 8-Day L4 Global 500m SIN Grid V006[END_REF].

The DEM at 90 m resolution was derived from SRTM of National Aeronautics and Space Administration (NASA) [START_REF] Andy | Hole-filled SRTM for the globe Version 4, available from the CGIAR-CSI SRTM 90m Database[END_REF]. Based on the elevation, we used SAGA GIS to derive 45 relief derivatives [START_REF] Conrad | System for Automated Geoscientific Analyses (SAGA) v. 2.1.4[END_REF], such as Analytical Hillshading (AH), Slope (SLO), Channel Network Base Level (CNBL), Plan Curvature (PLC), Multiresolution Index of Valley Bottom Flatness (MRVBF) and Topographic Wetness Index (TWI).

Based on the latitude and longitude, oblique geographic coordinates (OGC) were calculated at 15 • , 30 , 150 • and 165 • angles, with the following equation [START_REF] Møller | Oblique geographic coordinates as covariates for digital soil mapping[END_REF]:

OGC = ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ Lat 2 + Lon 2 √ × cos ( α -tan -1 (Lat/Lon) ) (2)
where Lon and Lat are the longitude and latitude, α is the angle (in degrees) of the titled axis relative to the x axis.

All environmental covariates were resampled to 250 m resolution using the bilinear method for spatial modelling and map prediction. The collection and pre-processing of remote sensing data were performed in Google Earth Engine [START_REF] Gorelick | Google Earth Engine: Planetary-scale geospatial analysis for everyone[END_REF]. All covariates were unified into the CGCS WGS 1984 geographic coordinate system, and predicted maps were projected into the Albers Conic Equal Area coordinate system.

Variable selection methods

Boruta

The Boruta algorithm is one of the widely used variable selection methods in DSM [START_REF] Xiong | Holistic environmental soil-landscape modeling of soil organic carbon[END_REF][START_REF] Xu | Incorporation of satellite remote sensing pan-sharpened imagery into digital soil prediction and mapping models to characterize soil property variability in small agricultural fields[END_REF][START_REF] Keskin | Digital mapping of soil carbon fractions with machine learning[END_REF][START_REF] Rasaei | Rescue and renewal of legacy soil resource inventories in Iran as an input to digital soil mapping[END_REF]. Proposed by [START_REF] Kursa | Feature selection with the Boruta package[END_REF], Boruta first duplicates the data set and then shuffles its predictors in each column (which are called shadow predictors). Afterwards, it trains a Random Forest (RF) model using the original and shuffled data sets combined and evaluates the variable importance (Z score) for each predictor. In each iteration, it checks whether a real predictor has higher importance (RMSE) than the best of its shadow predictors and marks the predictor as either confirmed (important) or rejected (unimportant). Finally, it stops when all the predictors are confirmed or rejected. The Boruta algorithm was performed in R package "Boruta" [START_REF] Kursa | Feature selection with the Boruta package[END_REF].

Recursive feature elimination

The RFE algorithm is commonly used to select the most relevant predictors for machine learning methods [START_REF] Gomes | Modelling and mapping soil organic carbon stocks in Brazil[END_REF][START_REF] Chen | Digital mapping of the soil thickness of loess deposits over a calcareous bedrock in central France[END_REF][START_REF] Poggio | SoilGrids 2.0: producing soil information for the globe with quantified spatial uncertainty[END_REF][START_REF] Hounkpatin | Assessment of the soil fertility status in Benin (West Africa)-Digital soil mapping using machine learning[END_REF]. Based on the backward selection, RFE works as follows: (1) fit a model using all the n predictors, calculate model performance by k-fold cross-validation (RMSE) and the variable importance; (2) remove the least important predictor from the pool, refit the model, assess model performance and remove the least important predictor again; (3) repeat the second procedure down to a pool from n to 1 with a step of 1; (4) determine the optimal number of predictors by taking the model with the best performance (RMSE). The RFE was implemented in R package "caret" using RF as an internal model [START_REF] Kuhn | caret: Classification and Regression Training. R package version 6[END_REF].

Variance inflation factor analysis

The VIF analysis is a commonly used method to evaluate the multi-collinearity among predictors and select the optimal subset of predictors that are not correlated [START_REF] Curto | The corrected vif (cvif)[END_REF]. During the VIF process, each predictor is regressed against the remaining predictors by a linear model using ordinary least squares regression. After the determination of the coefficient of determination (R 2 ) from the fitted linear model, the VIF is calculated by equation 3:

VIF i = 1 1 -R 2 i (3)
where R 2 i is the R 2 for predictor i. According to previous studies, predictors with VIF > 10 were determined to be highly correlated with other predictors and were therefore removed from the final DSM model [START_REF] Odhiambo | Spatial prediction and mapping of soil pH across a tropical afro-montane landscape[END_REF][START_REF] Taghizadeh-Mehrjardi | Enhancing the accuracy of machine learning models using the super learner technique in digital soil mapping[END_REF]. The R package "car" was used to perform VIF analysis [START_REF] Fox | An R Companion to Applied Regression, Third Edition[END_REF].

Modified greedy feature selection

Based on the GFS algorithm [START_REF] Drobnič | On the interpretability of machine learning models and experimental feature selection in case of multicollinear data[END_REF], a modified greedy feature selection (MGFS), is proposed for DSM regression problem. Opposite to RFE, MGFS adopts a forward selection strategy which includes steps (Fig. 2): (1) fit a model using all the n predictors, and calculate the variable importance; (2) select the most important predictor (only one) to fit an initial model, and calculate the model performance (RMSE) by k-fold cross-validation (note that there is only one predictor in the pool); (3) fit a list of models using 2 predictors (the combinations of predictor(s) in the pool and one of the remaining predictors), calculate their model performance, and record the model with the best performance; (4) update the pool by taking the predictors from the best model in the previous step; (5) repeat steps 3 and 4 by increasing the number of predictors from 3 to n. The predictors in the model with the best performance (RMSE) are selected for the final model. It is possible to set an early stop when the model performance starts to decrease for a large number of predictors (>50). Since Boruta and RFE both used RF, the same internal model was used in MGFS for a fair comparison in this study. The differences between original GFS and MGFS are in two aspects:

(1) GFS adds an additional constraint (least-trees-used criterion) when adding predictors to the model. This additional constraint likely has little impact since a tie would be required between 2 predictors, which is highly unlikely, and therefore it is not adopted in MGFS;

(2) Inspired by RFE, MGFS started with the most important predictor based the fitted model using all the predictors, while GFS tested all the single predictor in the first iteration. Therefore, MGFS should be slightly more efficient than GFS. In summary, VIF and Boruta are fundamentally different from RFE and MGFS. VIF operates independently from the dependent predictor and deals primarily with multicollinearity, and Boruta is based on individual predictor importance with regards to shadow predictors. Solely based on the variable importance, RFE and MGFS are equivalent.

Calibration model

Quantile regression forest (QRF) has seen increased use in DSM studies, especially at a broad scale, because it provides calculation of uncertainty and good model performance [START_REF] Meinshausen | Quantile regression forests[END_REF][START_REF] Vaysse | Using quantile regression forest to estimate uncertainty of digital soil mapping products[END_REF][START_REF] Loiseau | Satellite data integration for soil clay content modelling at a national scale[END_REF][START_REF] Kasraei | Quantile regression as a generic approach for estimating uncertainty of digital soil maps produced from machine-learning[END_REF][START_REF] Poggio | SoilGrids 2.0: producing soil information for the globe with quantified spatial uncertainty[END_REF].

We define X and Y as the predictor variables and target variables, QRF generates a large number of trees (b) using bootstrapping (random sampling with replacement) from p training samples (X i ,Y i ), i = 1,…, p. A random subset of the predictor variables is then used to select split-point for each node of the bootstrap tree. For a new sample N = X n , its prediction for each bootstrap tree is the conditional mean estimate of Y, which can be formulated as below:

X = ∑ p i=1 w i Y i (4)
where w i is the weight for the sample (X i ,Y i ) in the same leaf within the same bootstrap tree.

The mean prediction of b bootstrap trees is used to represent the final prediction of the new sample N. Using the weighted samples, QRF can also derive a conditional distribution from which the probability of Y being lower than a given percentile can be determined and thus calculate the prediction intervals. We refer to [START_REF] Meinshausen | Quantile regression forests[END_REF] for more details relevant to the calculation of conditional distribution. The variables to possibly split at each node (mtry) were optimized by 5-fold cross-validation with grid searching, and other parameters were set to default values as suggested by [START_REF] Kuhn | caret: Classification and Regression Training. R package version 6[END_REF]. The "caret" [START_REF] Kuhn | caret: Classification and Regression Training. R package version 6[END_REF] and "quantregForest" (Meinshausen, 2017) packages were used for optimizing and running QRF in R (R Core Team, 2021). Predictions at mean, 5th and 95th percentiles can be derived from the fitted QRF model directly [START_REF] Vaysse | Using quantile regression forest to estimate uncertainty of digital soil mapping products[END_REF][START_REF] Poggio | SoilGrids 2.0: producing soil information for the globe with quantified spatial uncertainty[END_REF].

Evaluation of model performance and computation efficiency

The whole dataset (402 samples) was randomly split into calibration (281 samples) and validation (121 samples) sets at a ratio of 70 % to 30 % with 50 repeats. The R 2 , Lin's concordance correlation coefficient (LCCC), and RMSE were used to evaluate model accuracy on the validation set. Higher R 2 and LCCC close to 1 and a lower RMSE close to 0 mean better model accuracy.

R 2 = 1 - ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ ̅ ∑ n i (ŷ i -y i ) 2 √ ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ ∑ n i (y i -y) 2 √ (5) RMSE = ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ ̅ ∑ n i (ŷ i -y i ) 2 n √ √ √ √ √ (6) LCCC = 2ρσ y σ ŷ (y -ŷ) 2 + σ 2 y + σ 2 ŷ (7)
where y i and ŷi are observation and prediction for sample i, y and ŷ are the mean of all the observations and predictions, n is the number of samples, ρ is the correlation coefficient, σ y and σ ŷ are the variances of all the observations and predictions. The 90th prediction interval coverage probability (PICP90), calculated by the proportion of observations that fall within the 90 % prediction intervals (PIs), was used to determine the performance of the uncertainty estimate on the validation set [START_REF] Malone | Comparing regression-based digital soil mapping and multiple-point geostatistics for the spatial extrapolation of soil data[END_REF]. A PICP90 close to 0.9 (or 90 %) means the uncertainty estimate is efficacious. To determine the magnitude of uncertainty, we calculated the uncertainty index (UI) by equation 8 [START_REF] Viscarra Rossel | Baseline map of organic carbon in Australian soil to support national carbon accounting and monitoring under climate change[END_REF]Zhou et al., 2019b). A greater UI indicates higher model uncertainty.

UI = 1 n ∑ n i=1 ( Q95 i -Q5 i Mean ) (8)
where Mean, Q95 i and Q5 i are the prediction at mean, 95th and 5th percentiles for sample i in the validation set.

The computation efficiency was determined by the computation time using a DELL Precision 3650 Tower Workstation (Intel i9-10900 K CPU with 10 cores and 64 GB RAM). Here we recorded the computation time in two steps, variable selection and map prediction, which are the most time-consuming. The total number of pixels to be predicted was about 8.7 million for the cropland in the study area, and parallel computation was not adopted in map prediction.

Results

Statistical summary of soil properties

Table 2 shows the statistical summary of SOCD in all, calibration and validation datasets. For all dataset, SOCD ranged from 10.41 to 149.04 t ha -1 . The median and mean of SOCD were 41.29 and 46.26 t ha -1 at 0-30 cm. According to the coefficient of variation (47.65 %), SOCD had high heterogeneity. The skewness of 1.61 and kurtosis of 6.36 indicated that SOCD had a slightly positive and leptokurtic distribution in all dataset. The calibration and validation for SOCD had similar distributions to all dataset, showing that the repeated random data split (50 times) was reasonable for proper model calibration and validation.

Number of covariates selected by variable selection methods

Table 3 indicates the number of final covariates selected by the four variable selection methods. The optimal covariates were determined by Z score for Boruta, VIF score (<10) for VIF, and best model performance (RMSE) for RFE and MGFS (Figures S1,S2,S3 and S4). Among 392 environmental covariates, Boruta retained 76 covariates, within which 3, 26, 36, 6 and 5 covariates were related to soil, climate, organisms, relief and position. VIF kept 30 covariates linked to soil (3), climate (4), organisms (21), and relief (2). RFE selected 22 covariates, including 3, 5, 10, and 4 covariates relevant to soil, climate, organisms, and position, respectively. Selecting the minimum number of covariates, MGFS only kept 9 covariates for the final modelling, including BI (Brightness Index, soil), Tavg (Mean annual temperature, climate), UH (Upslope height, relief), OGC105 (Oblique geographic coordinate at 105 • , position), NDVI3Min (Minimum NDVI at 3rd quarter, organisms), LAI2Min (Minimum LAI at 2nd quarter, organisms), NDRI1Avg (Average NDRI at 1st quarter, organisms), GARI3Max (Maximum GARI at 3rd quarter, organisms) and GPP3Avg (Average GPP at 3rd quarter, organisms) ranked descending by importance.

Model performance of four variable selection methods

The model performance using all 392 covariates and selected covariates by the four variable selection methods is presented in Fig. 3. VIF had the lowest model performance with R 2 of 0.48, LCCC of 0.64, and RMSE of 15.79 t ha -1 . The R 2 , LCCC and RMSE for the model using all the covariates (named Full hereafter) were 0.54, 0.70 and 14.79 t ha -1 in SOCD prediction. Compared to the Full model, RFE performed slightly better with R 2 of 0.57, LCCC of 0.72, and RMSE of 14.23 t ha -1 , while Boruta had a lower model performance with lower R 2 (0.55), LCCC (0.71) and greater RMSE (14.66 t ha -1 ) . Among all the models, MGFS performed best with R 2 of 0.60, LCCC of 0.74 and RMSE of 13.80 t ha -1 . Regarding the PICP, all the models were around 0.90 (0.88-0.92), indicating a similar quantification ability of 90 % PIs. The UI indicated that MGFS had the lowest global uncertainty (0.86) while the VIF model had the highest global uncertainty (0.98).

Computation efficiency of three variable selection methods

Table 4 shows the computation efficiency of variable selection and map prediction using the Full model and four variable selection methods. Since the Full model did not involve variable selection, the total computation time was 364 min which was equal to the time of map prediction. In the step of variable selection, Boruta (17 min) was the most efficient method, followed by MGFS (67 min), VIF (100 min) and RFE (150 min). In the step of map prediction, MGFS (34 min) performed much more efficient than RFE (51 min), VIF (59 min) and Boruta (102 min). Summing up two procedures, MGFS was the most efficient model (101 min), followed by Boruta (119 min), VIF (159 min), RFE (201 min), and the Full model (364 min).

Table 2

Descriptive statistics for SOCD (t/ha) in the all, calibration and validation datasets. The statistics for calibration and validation datasets are the mean of 50 replicates. Number of samples (No), minimum (Min), first quantile (Q1), third quantile (Q3), maximum (Max), standard deviation (SD), coefficient of variation (%CV), skewness (Skew) and kurtosis (Kurt). Fig. 3. Boxplots of QRF model performance with 50 repeats using all the predictors (Full), Boruta, recursive feature elimination (RFE), variance inflation factor (VIF), and modified greedy feature selection (MGFS).

X. Zhang et al.

Geoderma 432 ( 2023) 116383 8

Spatial distribution of SOCD and associated uncertainty

The spatial distribution of SOCD and 90 % PIs are shown in Fig. 4 and Figure S5. A general increasing gradient was observed in all the SOCD maps from the south-western part (20-50 t ha -1 ) to the north-eastern part (>70 t ha -1 ), and the difference among the four SOCD maps was relatively small. The lower and upper limits (5 % and 95 % PI) of 90 % PIs had a similar spatial pattern, while the lower limits were general-ly<50 t ha -1 and upper limits were commonly>57 t ha -1 over the whole study area.

Fig. 5 presents the spatial prediction of UI. The general spatial patterns of UI using five models were similar and the area with a lower density of calibration samples had greater UI (Fig. 1). It was also evident from the frequency plots that MGFS had fewer pixels with UI > 1.6 and more pixels with UI < 0.7, indicating a low global uncertainty.

Discussion

Comparison of variable selection methods

Our results showed that MGFS selected the most parsimonious model with only 9 covariates from a total of 392 covariates compared to Boruta (76 variables), VIF (30 variables), and RFE (22 variables) (Table 3). Among 9 covariates selected by MGFS, BI was the most important one sensitive to the brightness of soils. [START_REF] Liu | A soil colour map of China[END_REF] noted that soil colour is dominated by SOC in Northeast and North China due to the cold continental semi-humid temperature regime. Therefore, being a proxy of soil colour, BI can provide crucial information in SOCD modelling. The good correlation between BI and SOC (or SOCD) was also confirmed by previous studies [START_REF] Gholizadeh | Soil organic carbon and texture retrieving and mapping using proximal, airborne and Sentinel-2 spectral imaging[END_REF][START_REF] Ayala Izurieta | Multi-predictor mapping of soil organic carbon in the alpine tundra: a case study for the central Ecuadorian páramo[END_REF][START_REF] Wang | Prediction of Soil Organic Carbon under Different Land Use Types Using Sentinel-1/-2 Data in a Small Watershed[END_REF]. Mean annual temperature was the second important covariate in SOCD modelling. Previous studies have demonstrated the importance of temperature in SOC variation because temperature predominantly affects the microbial decomposition of SOC [START_REF] Wiesmeier | Soil organic carbon storage as a key function of soils-A review of drivers and indicators at various scales[END_REF]. A low temperature leads to low SOC decomposition, so numerous studies indicated an increase in SOC with decreasing temperature (e.g., [START_REF] Koven | Higher climatological temperature sensitivity of soil carbon in cold than warm climates[END_REF][START_REF] Chen | National estimation of soil organic carbon storage potential for arable soils: A data-driven approach coupled with carbon-landscape zones[END_REF]. The latitude of the study area expands from 30 • N to 48 • N, leading to an apparent temperature gradient from the south-western part (16 • C) to the northeastern part (-1.1 • C) that influences SOC decomposition. The upslope height (relief factor) is relevant to the potential energy entering the Fig. 4. Spatial distribution of SOCD produced by QRF using modified greedy feature selection (MGFS) and its associated 90% prediction intervals.

X. Zhang et al.

surface and indirectly affects the carbon input into the soil [START_REF] Thompson | Digital Soil Mapping: Interactions with and Applications for Hydropedology[END_REF]. The importance of OGC at 105 • indicates a spatial trend of SOCD at 105 • relative to the × axis [START_REF] Møller | Oblique geographic coordinates as covariates for digital soil mapping[END_REF]. The remaining 5 covariates (NDVI3Min, LAI2Min, NDRI1Avg, GARI3Max, GPP3Avg) were related to the organisms factor. It is interesting to note that all these covariates were calculated within a given quarter. This finding suggests that the remote sensing derived indicators calculated within a shorter time window (e.g., month, quarter) can capture more helpful vegetation phenology information related to soil properties than annual ones. This finding is in line with previous results from [START_REF] Hengl | SoilGrids250m: Global gridded soil information based on machine learning[END_REF], [START_REF] Keskin | Digital mapping of soil carbon fractions with machine learning[END_REF] and [START_REF] Kempen | Mapping topsoil organic carbon concentrations and stocks for Tanzania[END_REF].

Being the most parsimonious model, MGFS also performed best among four variable selection methods with R 2 of 0.60, LCCC of 0.74, and RMSE of 13.80 t ha -1 (Fig. 3, Table 3). Compared to the Full model (R 2 of 0.54, LCCC of 0.70, RMSE of 14.49 t ha -1 ), MGFS increased 11.1 % and 5.7 % in R 2 and LCCC, and decreased 6.8 % in RMSE. RFE and Boruta performed slightly better than the Full model while VIF had the lowest model performance. The unsatisfactory result of VIF is rooted in the fact that VIF only accounted for the correlation between predictors and neglected their relationships to the target variable. In addition, VIF eliminates multicollinearity but it also removed part of the useful information between correlated variables, leading to a lower model performance. The implication of this result is that highly correlated variables are not mandatory to be removed before fitting a machine learning model (i.e. QRF) because part of their uncorrelated information can provide useful information for modelling. Our result demonstrates that variable selection methods can efficiently decrease the number of covariates, but they do not ensure improving model performance simultaneously (Chen et al., 2022b;[START_REF] Luo | Regional soil organic matter mapping models based on the optimal time window, feature selection algorithm and Google Earth Engine[END_REF]. Therefore, digital soil mappers should compare the Full model with the parsimonious model derived from the selected variable selection method and then decide whether the latter is preferable. There was no considerable difference in 90 % prediction intervals since all the PICP were around 0.9 (0.88-0.92), indicating all these uncertainty estimates were efficacious. The UI indicated that MGFS and RFE had the lowest global uncertainty (0.86) when compared to Boruta (0.88), Full model (0.91) and VIF (0.98). These results indicated that we need to include UI for reporting the magnitude of uncertainty in DSM studies since it is beyond the capacity of PICP (e.g., Viscarra Rossel et al., 2014;Zhou et al., 2019b;[START_REF] Poggio | SoilGrids 2.0: producing soil information for the globe with quantified spatial uncertainty[END_REF].

From the aspect of computation efficiency, MGFS was the most efficient variable selection method (101 min) when considering the variable selection (392 variables, 67 min) and map prediction (8.7 million pixels, 34 min) together (Table 4). When improving the map resolution from 250 m to 30 m, it can be expected that the superiority of MGFS in computation efficiency would be amplified. As noted by Chen et al. (2022a), benefiting from rapidly growing remote sensing data and new soil observations obtained from multiple platforms, the DSM studies are moving forward to a broader spatial extent with a finer spatial resolution. For fine-resolution DSM at a broad scale, map prediction is the most time-consuming step. Take SoilGrids 2.0 for example, the total computation time (RFE, model training and map prediction) for a single soil property was 1500 CPU hours of which map prediction accounted for about two-thirds [START_REF] Poggio | SoilGrids 2.0: producing soil information for the globe with quantified spatial uncertainty[END_REF]. MGFS can efficiently promote computation efficiency by constructing a parsimonious model with a fair model performance, so it has a high potential for broad-scale DSM studies at a fine resolution.

Comparison SOCD product with SoilGrids 2.0

We compared our SOCD product with SoilGrids 2.0 mapped for the globe at 250 m resolution in Fig. 6 [START_REF] Poggio | SoilGrids 2.0: producing soil information for the globe with quantified spatial uncertainty[END_REF]. After evaluating the map accuracy using the same validation sets, we found that SoilGrids 2.0 had R 2 of 0.24, LCCC of 0.60 and RMSE of 18.84 t ha -1 , which indicated a poor accuracy of SOCD estimation compared to our product (R 2 of 0.60, LCCC of 0.74, RMSE of 13.80 t ha -1 ). The spatial distribution maps showed that our SOCD estimate was generally higher than SoilGrids 2.0, especially in the northern part and the centre of the southwestern part as illustrated in the difference map between our SOCD product and SoilGrids 2.0. The frequency plot indicated that>88 % of the region was underestimated by SoilGrids 2.0. The estimated SOC stock was 2.89 Pg for the study area in our product while the SOC stock estimation from SoilGrids 2.0 product was 24.2 % lower (2.19 Pg). The poor SOCD estimate of SoilGrids 2.0 resulted from the fact that the soil database (>140,000 locations) used for producing SoilGrids 2.0 only had<80 sampling sites in the study area, leading to a limited weight in the spatial predictive model to represent the unique pedo-climatic condition for the study area. Therefore, when using the global soil maps for regional or national studies, it is necessary to correct them by the regional or national soil observations, especially for these regions with limited soil observations in producing the global soil maps. For this purpose, previous studies have proved that model averaging was a potential solution for producing more accurate "local" digital soil property maps by coupling "global" soil maps and "local" soil observations [START_REF] Malone | Using model averaging to combine soil property rasters from legacy soil maps and from point data[END_REF][START_REF] Clifford | Combining two soil property rasters using an adaptive gating approach[END_REF]Guo, 2015, Caubet et al., 2019;[START_REF] Chen | Model averaging for mapping topsoil organic carbon in France[END_REF].

Limitations and perspectives

One limitation of MGFS is that it may take a long time to select the optimal variables for a large number of covariates (e.g., n > 50) as it needs to run the model (n 2 + n)/2 times. Therefore, the early stop is suggested in this study (see details in section 2.3.3) when using MGFS for a large set of covariates. It should be noted that the use of early stop may result in a local optimum which can be slightly different from the global optimum. From our personal experience, the local optimum was quite close to the global optimum but it still needs to be confirmed by other datasets [START_REF] Xiao | Improving pedotransfer functions for predicting soil mineral associated organic carbon by ensemble machine learning[END_REF].

From theory, MGFS method can work on all the machine learning methods, such as RF, support vector machine, gradient boosted machine, XGBoost. In this study, we only evaluated its performance on RF, so its applicability remains to be tested for other machine learning methods in future studies.

Conclusions

In this study, we proposed the use of a modified greedy feature selection (MGFS), to improve model parsimony and keep model accuracy simutanouly for digital soil mapping (DSM) regression. Taking soil organic carbon density (SOCD) mapping of cropland in Northeast and North China as a case study, we evaluated the model performance, uncertainty estimate and computation efficiency on MGFS together with most currently used variable selection methods Boruta, recursive feature elimination (RFE), and variance inflation factor (VIF). Our results showed that MGFS selected the most parsimonious model with only 9 covariates from a list of 392 covariates while RFE, VIF and Boruta retained 22, 30 and 76 covariates, respectively. The repeated validation (50 times) indicated that the model using MGFS selected covariates performed much better (R 2 of 0.60, LCCC of 0.74, and RMSE of 13.80 t ha -1 ) in SOCD mapping than VIF, Boruta and RFE. When compared to the model using all 392 covariates (R 2 of 0.54, LCCC of 0.70, RMSE of 14.79 t ha -1 ), MGFS had an increase of R 2 and LCCC at 11.1 %, 5.7 % and decrease of RMSE at 6.8 %. In addition, MGFS had the lowest global uncertainty (UI of 0.86) and the highest computation efficiency among all the models. Accordingly, we concluded that MGFS was superior to RFE, Boruta, VIF, and Full model in terms of model parsimony, model performance, map uncertainty and computation efficiency; therefore, it has a high potential in fine-resolution DSM studies at broad scales. Theoretically, MGFS can work on all machine learning methods; however, further tests on other datasets and machine learning models are still needed to evaluate its robustness.

Fig. 1 .

 1 Fig. 1. The location of study area (a) and the distribution of soil sampling sites (b).

Fig. 2 .

 2 Fig. 2. Diagram of modified greedy feature selection (MGFS).

Fig. 5 .

 5 Fig. 5. Uncertainty index of SOCD produced by QRF using all the predictors(Full) (a), Boruta (b), recursive feature elimination (RFE) (c), variance inflation factor (VIF) (d) and modified greedy feature selection (MGFS) (e).

Fig. 6 .

 6 Fig. 6. Spatial distribution of SOCD produced by modified greedy feature selection (MGFS) (a) and SoilGrids 2.0 (b), and difference between MGFS and SoilGrids 2.0 (c).

X

  .Zhang et al.

Table 1

 1 List of environmental covariates used in this study. The number in the columns of Scorpan factor and Covariate indicates the number of variables within each group.

						Table 1 (continued )
						Scorpan	Covariate	Abbreviation	Scale	Reference
						factor
	Scorpan factor	Covariate	Abbreviation	Scale	Reference					Vescovo and Gianelle (2008)
	Soil (15) Climate (68) Organisms (2 5 3)	Soil Erosion (1) Clay, Silt and Sand Content (3) Brightness Index (1) Carbonate Index (1) Coloration Index (1) Ferrous Minerals (1) Gypsum index (1) Hue Index, Saturation Index (1) Iron Oxide (1) Reflectance Absorption Index (1) Redness Index (1) Stress Related (1) Bioclimatic Variables (Mean Diurnal Range, Isothermality, Temperature Seasonality, Wind Speed, Vapor Pressure, Solar Radiation, Average Temperature (Average, Maximum, Minimum) (3) Day-time and Night-time Land Surface Potential Evapotranspiration Normalized Difference Snow Index (21) b Canopy Index (21) b (11) a Evapotranspiration (11) a and Temperature (2) Precipitation (4) Coldest Quarter) (16) Precipitation of of Warmest Quarter, Quarter, Precipitation Precipitation of Driest Wettest Quarter, Precipitation of Seasonality, Month, Precipitation Precipitation of Driest Wettest Month, Precipitation of Annual Precipitation, Coldest Quarter, Mean Temperature of Warmest Quarter, Temperature of Quarter, Mean Temperature of Driest Wettest Quarter, Mean Temperature Of Range, Mean Temperature Annual	SE Clay, Silt, Sand BI CarI Coll FM GI HI, SI IO RAI RI SR Tavg, Tmax, Tmin NDSI CANI PET, ET LSTD, LSTN Wind, VP, Sol, Prec BIO02, 03, 04, 07, 08, 09, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19	1000 m 250 m 30 m 30 m 30 m 30 m 30 m 30 m 30 m 30 m 30 m 30 m 1000 m 30 m 30 m 500 m 1000 m 1000 m 1000 m	Teng et al. (2019) Hengl et al. (2017) Abbas and Khan (2007) Boettinger et al. (2008) Escadafal (1994) Imbroane et al. (2007) Boettinger et al. (2008) Mandal (2016) Hewson et al. (2001) Ghaemi et al. (2013) Madeira et al. (1997) Foody et al. (2001) (1994) Riggs et al. et al. (2017) Running (2021) Wan et al. (2017) Hijmans Fick and (2017) Hijmans Fick and Fick and Hijmans (2017)	Relief (46)	Differenced Vegetation Index (21) b Enhanced Vegetation Index (21) b Green Atmospherically Resistant Vegetation Index (21) b Normalized Difference Red/Green Redness Index (21) b , Ratio Vegetation Index (21) b Normalized Difference Vegetation Index (21) b Normalized Difference Water Index (21) b Soil Adjusted Vegetation Index (21) b Gross (11) a and Net Primary Production (11) a Fraction of Diurnal Anisotropic Terrain Ruggedness, Valley Depth, Vector Upslope Height,, Wetness Index, Index, Topographic Topographic Position Power Index, Flatness, Stream of Valley Bottom Multiresolution Index of Ridge Top Flatness, Multiresolution Index Ruggedness Number, Position, Melton Length, Relative Slope Length Factor, Slope Upslope), Slope Plan, Profile, Maximum, Minimum, Longitudinal, Local Upslope, Local Downslope, Downslope, Local, Classification, (Cross-Sectional, Gradient, Curvature Geomorphons, Generalized Surface, Path Length, Accumulation, Flow Index, Flow Features, Convergence Morphometric Landform, Depressions, Clusters, Distance, Closed Channel Network Network Base Level, Area, Channel Area, Total Catchment Slope, Catchment Photosynthetically Active Radiation (21) b , Leaf Area Index (21) b Elevation, Analytical Hillshading, Aspect,	DVI EVI GARI NDRI, RVI NDVI NDWI SAVI GPP, NPP FPAR, LAI DAH, WE, SI MAH, VD, VTR, TPI, TWI, UH, MRVBF, SPI, MRN, MRRTF, LSF, SL, RSP, PLC, PRC, UC, LTC, MAC, MIC, LC, LDC, LUC, CSC, CVC, DSC, GS, GEO, GRA, CGI, FA, FPL, CLU, LAN, MF, CNBL, CND, CD, ELE, AH, ASP, SLO, CA, TCA,	30 m 30 m 30 m 30 m 30 m 30 m 30 m 500 m 500 m 90 m	Richardson and Wiegand (1977) Huete et al. (2002) Gitelson et al. (1996) Bannari et al. (1995) Tucker (1979) Gao (1996) Huete (1988) Running et al. (2015) Myneni et al. (2021) Jarvis et al. (2008)

(continued on next page) X.

Zhang et al.

Table 1

 1 (continued )

	Scorpan	Covariate	Abbreviation	Scale	Reference
	factor
		Heating, Wind
		Exposition and Shelter
		Index (46)
	Position	Oblique geographic	OGC15, 30, 45,	30 m	Møller et al.
	(10)	coordinate at 15 • , 30 • ,	60, 75, 105,		(2020)
		45 • , 60 • , 75 • , 105 • ,	120, 135, 150,
		120 • , 135 • , 150 • , 165 •	165
		(10)
	a Mean and standard deviation of the intra-annual and intra-quarterly ag-
	gregation in 2003-2017 (10 variables), and averaged standard deviation of
	intra-quarterly aggregation in 2003-2017 (1 variable). b Intra-annual and intra-quarterly mean and standard deviation in 2003-2017
	(10 variables), mean of the intra-annual and intra-quarterly maximum and
	minimum in 2003-2017 (10 variables), and averaged intra-quarterly standard
	deviation in 2003-2017 (1 variable).

Table 3

 3 Number of covariates after variable selection.

	Category		Number of covariates
		RFE	Boruta	VIF	MGFS
	Whole	22	76	30	9
	Soil	3	3	3	1
	Climate	5	26	4	1
	Organisms	10	36	21	5
	Relief	0	6	2	1
	Position	4	5	0	1

Table 4

 4 Computation efficiency of feature selection and map prediction (in minutes). The most efficient method is shown in bold font for each procedure.

	Method	Feature selection	Map prediction	All procedures
	Full	/	364	364
	Boruta	17	102	119
	RFE	150	51	201
	VIF	100	59	159
	MGFS	67	34	101

Acknowledgements

This study was supported by the National Natural Science Foundation of China (No. 41930754, 42001047, 41901055) and Ten-thousand Talents Plan of Zhejiang Province, China (2019R52004). Dominique Arrouays (DA) and Anne C. Richer-de-Forges (ARdF) are members of the research consortium GLADSOILMAP supported by LE STUDIUM Institute of Advanced Research Studies, France. ARdF is coordinator and DA is member of the French "Centre d'Expertise Scientifique Cartographie Numérique des sols" granted by the CNES-TOSCA program. We also would like to acknowledge two anonymous reviewers for their constructive suggestions which improved our manuscript a lot.

Data availability

The authors do not have permission to share data.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi. org/10.1016/j.geoderma.2023.116383.