Soil Salinity Mapping of Plowed Agriculture Lands Combining Radar Sentinel-1 and Optical Sentinel-2 with Topographic Data in Machine Learning Models

Abstract

This study assesses the relative performance of Sentinel-1 and -2 and their combination with topographic information for plow agricultural land soil salinity mapping. A learning database made of 255 soil samples’ electrical conductivity (EC) along with corresponding radar (R), optical (O), and topographic (T) information derived from Sentinel-2 (S2), Sentinel-1 (S1), and the SRTM digital elevation model, respectively, was used to train four machine learning models (Decision tree—DT, Random Forest—RF, Gradient Boosting—GB, Extreme Gradient Boosting—XGB). Each model was separately trained/validated for four scenarios based on four combinations of R, O, and T (R, O, R+O, R+O+T), with and without feature selection. The Recursive Feature Elimination with k-fold cross validation (RFEcv 10-fold) and the Variance Inflation Factor (VIF) were used for the feature selection process to minimize multicollinearity by selecting the most relevant features. The most reliable salinity estimates are obtained for the R+O+T scenario, considering the feature selection process, with R² of 0.73, 0.74, 0.75, and 0.76 for DT, GB, RF, and XGB, respectively. Conversely, models based on R information led to unreliable soil salinity estimates due to the saturation of the C-band signal in plowed lands.

1. Introduction

1.1. Soil Salinity: An Agriculture Threat

Soil salinization is a problem that affects agriculture worldwide [1], especially soil in dry sub-humid, semi-arid, and arid regions (taking into account the increasing degree of aridity) where salinity has reached alarming levels [2]. The presence of salts leads to a degradation of the soil’s physical properties by the dispersion of mineral colloids, as well as generating an osmotic potential in the soil, limiting the easy exchange of water and nutrients with the roots of the plants, which in turn results in a delay in the growth of the plants [3]. This process induces negative effects on crop productivity and quality [4], which can generate significant economic losses [5]. In the current context of a likely increase in temperature and decrease in precipitation in the south-central highlands [6], an increase in irrigation volumes is necessary to meet crop water requirements. However, irrigation without proper leaching and drainage leaves salt precipitates in the soil profile, while over-irrigation favors the gradual rise of the groundwater level, which in turn contributes to soil salt accumulation [7]. In this context, it is estimated that 50% of global cropland areas will be affected by salinization by 2050 [8], thus representing a major socioeconomic threat that can lead to population outmigration [4,9].

1.2. Soil Salinity Monitoring: Traditional vs. Remote Sensing Techniques

Conventional studies mapping soil properties such as salinity have the advantage of allowing a comprehensible representation of complex soil formation processes through conceptual models [10] and are not restricted by the quantity of soil types that can be mapped with a single conceptual model [11]. However, they require field campaigns that can be (i) slow, (ii) expensive, and (iii) limited in space and time. In this context, salinity data cannot be provided in a timely manner or at the agriculture plot scale to make decisions about planting and/or to understand the influence of agricultural practices on salinity variability [12,13]. To overcome these limitations, digital soil mapping based on remote sensing data emerges as a promising technique to update soil information in a cost-effective manner and with a regional perspective [14,15]. More specifically, soil salinity mapping has been substantially improved across the last decade through the integration of remote sensing data into machine learning (ML) models [2,16]. Freely available, with a global coverage and a high frequency revisit time (5–10 days), Sentinel-2 (S2) optical and/or Sentinel-1 (S1) radar images are generally used as explanatory features in the ML models.

Recent studies have made it possible to map soil salinity on a global scale by combining soil information and remote sensing, which can then be used to estimate the area affected by salinity worldwide, covering about one billion hectares [1,17]. Many studies carried out in China [18], Algeria [19], Iran [20], and Bolivia [21] have reported on the potential of Sentinel-2 optical satellites images (S2) to retrieve soil salinity estimates with a relatively good degree of confidence. However, these techniques are not suitable for regions or time periods affected by clouds. In the presence of clouds, unlike radar sensors, optical sensors (such as S2) cannot observe the ground surface or, therefore, soil salinity content. In this context, several studies carried out in Iran [22], Vietnam [23], India [24,25], and Bolivia [26] reported a good degree of confidence for soil salinity estimates based solely on Sentinel-1 radar images (S1), whereas some authors took advantages of the S1 and S2 complementarity to improve soil salinity estimates [27]. Despite the significant advances in soil salinity mapping brought by the above-mentioned studies, those works are limited to undisturbed natural soils. However, agriculture is developed on plowed soils with different roughness and compaction features than the ones of natural undisturbed soils. As roughness and compaction features affect S2 reflectance and S1 polarization [28,29,30], soil salinity models calibrated across natural unperturbed land cannot be applied to plow soil. In this context, some authors focused on agriculture plots to assess the potential of S2 alone [31,32,33] and in combination with S1 [34,35] to estimate soil salinity in cultivated land.

1.3. Which Machine-Learning Model Performs the Best?

Independently of the type of soil (i.e., natural or agricultural), the choice of the machine learning model could affect the reliability of the soil salinity estimates. In natural undisturbed soils, the sensitivity of soil salinity estimates to different machine learning models, such as Artificial Neural Network (ANN), decision tree (DT), Random Forest (RF), Support Vector Machine (SVM), Convolutional Neural Network (CNN), Classification and Regression Trees (CART), Partial Least Squares Regression (PLSR), Multi Linear Regression (MLR), and k-Nearest Neighbor (kNN) models, has been extensively reported [18,19,36,37,38,39]. In most cases, more reliable soil salinity estimates are obtained when considering the RF or SVM models [21]. For agricultural lands, only one study has reported on the sensitivity of soil salinity estimates to different machine learning models [35]. Carried out in Egypt, this study compared the soil estimates obtained with three machine learning models (i.e., RF, SVM, and Back Propagation Neural Network (BPNN)). The results show that the most reliable soil salinity estimates are obtained with the RF model, as this model is more suitable for mapping heterogeneous regions, resulting in a more accurate representation of the spatial distribution of areas affected by salinity.

1.4. Study Objectives

Based on the previously established state of the art, this study aims to assess the potential of S1 and S2 images (i.e., separately and as a whole) for soil salinity estimates in plowed areas. Several machine learning models (i.e., RF, DT, Gradient Boosting, and XGBoost) are considered to expand the current limited knowledge on the sensitivity of agricultural soil salinity estimates to different machine learning models. With 900 km² of its area at risk of salinization [40], the Bolivian Altiplano is selected as the study site.

2. Materials

2.1. Study Area

The study was carried out in the central and southern Altiplano of Bolivia (Figure 1). The region is part of the TDPS endorheic system (Titicaca, Desaguadero, Poopó, Salar) located between the eastern and western Andean Cordilleras at a mean elevation of 4000 m m.a.s.l. with a very flat topography (i.e., mean slope of 4.7°) [41,42]. The region is semi-arid with (i) a mean precipitation of less than 600 mm∙year⁻¹ concentrated in the austral summer [43], (ii) a mean temperature ranging from 4 °C to 8 °C [44], and (iii) a high evapotranspiration rate estimated at 1700 mm∙year⁻¹ around the Poopó Lake [45,46]. Agriculture is a main source of economic income that increasingly relies on irrigation practices to cope with particularly adverse weather conditions [47]. The study area encompasses three main types of soil: Haplic Xerosols (Xh-northern part), Eutric Fluvisols (Je-surroundings of Poopó Lake), and Haplic Yermosols (Yh-southern part) [18]. According to the guidelines of the WRB System (World Reference Base for Soil), their approximate equivalents could be placed as Xerosols = Association of Leptosols − Durisols − Regosols, and Yermosols = Association of Regosols − Solonetz − Solonchaks [48]. The combination of a semi-arid climatic context and irrigation in this endorheic system favors the accumulation of salts on the soil surface, resulting in extreme saline soil conditions [7,21].

2.2. Reference Soil Salinity Data

A total of 255 soil samples (their location based on accessibility and the presence of plowed plots) were collected during the pre-sowing period (April to October), which correspond to the dry period with almost no precipitation events that could affect soil salinity (Figure 1c–e). The sample collections were made during 5 field campaigns lead in 2022 and 2023 (

Table 1. Soil sampling dates and image acquisition dates S1 and S2.

Mission	Date	Number of Samples	Date S1	Date S2
1	17 to 19 July 2022	17	18, 30 July 2022	17, 22, 27 July 2022
2	10 to 12 August 2022	18	11, 23 August 2022	11, 16, 21 August 2022
3	7 to 9 April 2023	57	8, 20 April 2023	8, 13, 18 April 2023
4	10 to 12 July 2023	79	13, 25 July 2023	12, 17, 22 July 2023
5	2 to 4 October 2023	84	5, 17 October 2023	5, 10, 15 October 2023
	Total:	255

).

Each soil sample corresponds to a composite sample made of 5 subsamples taken at an average depth of 10 cm at the center and four corners of each sampling site (approximately 10 m × 10 m in size) to take into account the potential soil salinity heterogeneity that may occur at the Sentinel images’ spatial resolution (10 m × 10 m).

The soil electrical conductivity (EC) of soil samples was measured in the soil laboratory of the Faculty of Agronomy of the Universidad Mayor de San Andrés (UMSA), as a proxy of soil salinity. To prepare the samples, they were first air-dried, pulverized, homogenized, and then sieved with a 2 mm mesh. The method used to measure the soil EC involved creating a suspension of soil and water using a 1:5 ratio method [49]. In this process, a soil sample weighing 5 ± 0.1 g was placed in a polypropylene centrifuge tube and 25 ± 0.5 mL of deionized water was added. The samples were then agitated on a stirring table for 60 min and allowed to decant for an additional 30 min. The EC (25 °C) was then measured in the supernatant without disturbing the sediment layer, using the Potentiometric Method, with a calibrated multi-parameter device fitted with a conductivity electrode.

The collected samples correspond to Haplic Yermosol (n = 11), Haplic Xerosol (n = 109) and Eutric Fluvisol (n = 135) (Figure 1f,g). Considering the FAO soil salinity classification [1], most of the soil samples analyzed do not present high levels of salinity (Figure 1f). Experienced farmers favor sowing in low salinity areas rather than high salinity areas to guarantee decent crop production. Soil samples with high salinity are mainly located in Eutric Fluvisol soils because of the presence of a nearby water table [40,50]. Indeed, regional groundwater is more saline in nature due to factors such as the dissolution of sediments from ancient Quaternary lakes, the recharge of the Desaguadero River, natural contamination due to the erosion of rocks with high concentrations of salts, and contamination due to centuries of mining activity [51,52]. Consequently, groundwater becomes more saline near Poopó Lake and Desaguadero River [52], as identified by the collected samples (Figure 1f).

2.3. Sentinel-2 Images and Pre-Processing

Sentinel-2 images (S2) provide reflectance in 12 bands with wavelengths ranging from the visible spectrum to the infrared region. For this study, B1 and B9 bands were not considered due to their low spatial resolution (60 m) and because their applications are mainly associated with water vapor and atmospheric correction (aerosols and water vapor).

S2 images are available prior (level 1C) and after (level 2A) the application of an atmospheric correction. To ensure time and space consistency in S2 reflectance values, the atmospherically corrected level 2A images are used. The “QA60” quality band was used to identify the pixels affected by cloud cover. All pixels identified as “cloudy” were systematically flagged as missing values. To reduce the effect of noise that could occur locally and obtain an average condition of the sampled site, composite images were elaborated for each field campaign [21,53]. The composite images are obtained considering the mean value of the three S2 images obtained during and just after the sampling date (Table 1). Images before the sampling date were not considered to avoid using images recorded when the agricultural plot may not yet have been plowed. S2 pre-processing (i.e., download and composite) was carried out using GEE [54] in the free Google Colab cloud platform (retrieved from https://colab.research.google.com/, accessed on 12 June 2023).

2.4. Sentinel-1 Images and Pre-Processing

The Sentinel-1 images (S1) provide polarization values in VV (Vertical/Vertical) and VH (Vertical/Horizontal) in interferometry mode with a central frequency of 5.405 GHz. S1 images are available as a Single Look Complex product (SLC), providing both phase and amplitude information, and as the post-processed Ground Range Detected product (GRD), providing direct surface characteristics in the form of intensity images. Due to its dual polarization antenna that could discriminate subtle changes in soil salinity content, the GRD product has become a valuable resource for mapping soils affected by salinity and therefore is being considered for the present study [55]. S1 images are available in both ascending and descending orbit. To ensure consistency in S1 observations in space and time, only the descending orbit is used.

S1 images were pre-processed according to four successive steps: (i) noise removal, (ii) radiometric calibration, (iii) topographic correction, and (iv) focal mean filter to reduce speckle effect. As for S2 images, composite images made of the two closest S1 post samples are considered to reduce the effect of noise that could occur locally (Table 1). S1 pre-processing (i.e., download, pre-process, and composite) was carried out using a similar way as in S2, via GEE [54] in the free Google Colab cloud platform (retrieved from https://colab.research.google.com/, accessed on 12 June 2023).

2.5. Digital Elevation Model

The Shuttle Radar Topographic Mission Digital Elevation Model (SRTM-DEM) was used in this study because it was previously identified as the most reliable DEM for the Altiplano region [41]. Data to produce the SRTM-DEM were collected over 11 days in February 2000 using dual spaceborne imaging radar (SIR-C) and dual X-band synthetic aperture radar (X-SAR). The SRTM-DEM has a spatial resolution at about 30 m and spans 60°N–56°S.

2.6. Machine Learning Models

Four supervised machine learning (ML) models based on decision trees are considered: (i) Random Forest (RF) [56], (ii) Extreme Gradient Boosting (XGB) [57], (iii) decision tree (DT) [58], and (iv) Gradient Boosting (GB) [59]. Decision trees are non-parametric supervised predictive models used for both regression and classification problems. They represent decisions and their possible consequences in a hierarchical manner similar to the structure of a tree. In this way, simple decision rules are inferred from the characteristics of the dataset to create a model that predicts the value of a target variable (i.e., soil salinity) [60].

In the DT model, a single decision tree is constructed by recursively dividing the dataset into smaller subsets based on the features that provide the greatest information gain or impurity reduction [58,60]. In the RF model, multiple decision trees are built using different subsets of the learning database (a technique called bagging). The final prediction corresponds to the average of the predictions of all trees (ensemble) [56]. In the GB model, decision trees are sequentially built so that the new tree attempts to correct the previous tree’s errors by minimizing a specific loss function. In this process, each tree is built on the residuals (errors) of the previous tree [59]. Finally, the XGB model is an optimized version of GB, which applies parallel computing techniques to speed up the training and regularization process to avoid overfitting [59].

The DT, RF, and GB models were implemented using the scikit-learn library [61], while XGB was implemented using the xgboost library, which is compatible with scikit-learn functions [57]. The default hyper-parameters assigned by the scikit-learn library in Python were considered for all models (100 decision trees excluding DT (n_estimator), and considering all features for each decision tree (max_features)). It is worth mentioning that these hyper-parameters provided good results in ML models [61].

3. Methods

3.1. Elaboration of the Learning Database

The datasets (i.e., S1, S2, and SRTM) were homogenized to a 10 m spatial resolution using the nearest-neighbor method. Then, S1 polarization (VV and VH), S2 spectral reflectance (B2, B3, B4, B5, B6, B7, B8, B8a, B11 and B12), and SRTM-DEM elevation values were extracted from the pixels where a soil sample was collected. At this stage, the learning database includes 255 EC observations with corresponding S1 polarization and S2 spectral reflectance values and SRTM-DEM elevation. Environmental features, including 24 spectral indices, 5 soil texture indices, and 12 topographic indices, were computed from S2 reflectance, S1 polarization, and SRTM-DEM elevation, respectively, and added to the learning database (

Table 2. Sentinel-2 features.

N	Index/Acronym	Formula * or Description	Reference
	Sentinel-2 bands	B2, B3, B4, B5, B6, B7, B8, B8A, B11, B12
1	Normalized Differential Vegetation Index (NDVI)	$\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}={\displaystyle \frac{\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{R}\mathrm{e}\mathrm{d}}{\mathrm{N}\mathrm{I}\mathrm{R}+\mathrm{R}\mathrm{e}\mathrm{d}}}$	[66]
2	Soil-Adjusted VegetationIndex (SAVI)	$\mathrm{S}\mathrm{A}\mathrm{V}\mathrm{I}={\displaystyle \frac{\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{R}\mathrm{e}\mathrm{d}}{\mathrm{N}\mathrm{I}\mathrm{R}+\mathrm{R}\mathrm{e}\mathrm{d}+0.5}}\times 1.5$	[67]
3	Enhanced Vegetation Index (EVI)	$\mathrm{E}\mathrm{V}\mathrm{I}=2.5\times {\displaystyle \frac{\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{R}\mathrm{e}\mathrm{d}}{\mathrm{N}\mathrm{I}\mathrm{R}+6\times \mathrm{R}\mathrm{e}\mathrm{d}-7.5\times \mathrm{B}\mathrm{l}\mathrm{u}\mathrm{e}+1}}$	[68]
4	Normalized Difference Moisture Index (NDMI)	$\mathrm{N}\mathrm{D}\mathrm{M}\mathrm{I}={\displaystyle \frac{\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1}{\mathrm{N}\mathrm{I}\mathrm{R}+\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1}}$	[69]
5	Moisture Stress Index (MSI1)	$\mathrm{M}\mathrm{S}\mathrm{I}1={\displaystyle \frac{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1}{\mathrm{N}\mathrm{I}\mathrm{R}}}$	[70]
6	Moisture Stress Index (MSI2)	$\mathrm{M}\mathrm{S}\mathrm{I}2={\displaystyle \frac{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1-\mathrm{N}\mathrm{I}\mathrm{R}}{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1+\mathrm{N}\mathrm{I}\mathrm{R}}}\times {\displaystyle \frac{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1-\mathrm{R}\mathrm{e}\mathrm{d}}{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1+\mathrm{R}\mathrm{e}\mathrm{d}}}$	Modified from [69]
7	Bare Soil Index (BSI)	$\mathrm{B}\mathrm{S}\mathrm{I}={\displaystyle \frac{(\mathrm{R}\mathrm{e}\mathrm{d}+\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1)-(\mathrm{N}\mathrm{I}\mathrm{R}+\mathrm{B}\mathrm{l}\mathrm{u}\mathrm{e})}{(\mathrm{R}\mathrm{e}\mathrm{d}+\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1)+(\mathrm{N}\mathrm{I}\mathrm{R}+\mathrm{B}\mathrm{l}\mathrm{u}\mathrm{e})}}$	[71]
8	Normalized Burn Ratio (NBR)	$\mathrm{N}\mathrm{B}\mathrm{R}={\displaystyle \frac{\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}2}{\mathrm{N}\mathrm{I}\mathrm{R}+\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}2}}$	[72]
9	Normalized Burn Ratio 2 (NBR2)	$\mathrm{N}\mathrm{B}\mathrm{R}2={\displaystyle \frac{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1-\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}2}{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1+\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}2}}$	[73]
10	Normalized Difference Water Index (NDWI)	$\mathrm{N}\mathrm{D}\mathrm{W}\mathrm{I}={\displaystyle \frac{\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{G}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{n}}{\mathrm{N}\mathrm{I}\mathrm{R}+\mathrm{G}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{n}}}$	[74]
11	Normalized Salinity Index (NSI1)	$\mathrm{N}\mathrm{S}\mathrm{I}1={\displaystyle \frac{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1-\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}2}{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1+\mathrm{N}\mathrm{I}\mathrm{R}}}$	[1]
12	Modified Normalized Salinity Index (NSI2)	$\mathrm{N}\mathrm{S}\mathrm{I}2={\displaystyle \frac{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1-\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}2}{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1-\mathrm{N}\mathrm{I}\mathrm{R}}}$	Modified from [1]
13	Salinity Index 1 (SI1)	$\mathrm{S}\mathrm{I}1=\sqrt{\mathrm{R}\mathrm{e}\mathrm{d}\times \mathrm{G}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{n}}$	[75]
14	Salinity Index 2 (SI2)	$\mathrm{S}\mathrm{I}2=\sqrt{\mathrm{B}\mathrm{l}\mathrm{u}\mathrm{e}\times \mathrm{R}\mathrm{e}\mathrm{d}}$	[14]
15	Salinity Index 3 (SI3)	$\mathrm{S}\mathrm{I}3=\sqrt{{\mathrm{R}\mathrm{e}\mathrm{d}}^2\times {\mathrm{G}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{n}}^2}$	[14]
16	Salinity Index 4 (SI4)	$\mathrm{S}\mathrm{I}4={\displaystyle \frac{\mathrm{N}\mathrm{I}\mathrm{R}\times \mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1-{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1}^2}{\mathrm{N}\mathrm{I}\mathrm{R}}}$	[1]
17	Salinity Index 5 (SI5)	$\mathrm{S}\mathrm{I}5={\displaystyle \frac{\mathrm{B}\mathrm{l}\mathrm{u}\mathrm{e}}{\mathrm{R}\mathrm{e}\mathrm{d}}}$	[14]
18	Salinity Index 6 (SI6)	$\mathrm{S}\mathrm{I}6={\displaystyle \frac{\mathrm{R}\mathrm{e}\mathrm{d}\times \mathrm{N}\mathrm{I}\mathrm{R}}{\mathrm{G}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{n}}}$	[14]
19	Vegetation Soil Salinity Index (VSSI)	$\mathrm{V}\mathrm{S}\mathrm{S}\mathrm{I}=2\times \mathrm{G}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{n}-5\times (\mathrm{R}\mathrm{e}\mathrm{d}+\mathrm{N}\mathrm{I}\mathrm{R})$	[76]
20	Normalized Difference Salinity Index (NSI)	$\mathrm{N}\mathrm{S}\mathrm{I}={\displaystyle \frac{\mathrm{R}\mathrm{e}\mathrm{d}-\mathrm{N}\mathrm{I}\mathrm{R}}{\mathrm{R}\mathrm{e}\mathrm{d}+\mathrm{N}\mathrm{I}\mathrm{R}}}$	[77]
21	Salinity Ratio (SR)	$\mathrm{S}\mathrm{R}={\displaystyle \frac{\mathrm{G}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{n}-\mathrm{R}\mathrm{e}\mathrm{d}}{\mathrm{B}\mathrm{l}\mathrm{u}\mathrm{e}+\mathrm{R}\mathrm{e}\mathrm{d}}}$	[76]
22	Canopy Response Salinity Index (CRSI)	$\mathrm{C}\mathrm{R}\mathrm{S}\mathrm{I}=\sqrt{{\displaystyle \frac{\mathrm{N}\mathrm{I}\mathrm{R}\times \mathrm{R}\mathrm{e}\mathrm{d}-\mathrm{G}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{n}\times \mathrm{B}\mathrm{l}\mathrm{u}\mathrm{e}}{\mathrm{N}\mathrm{I}\mathrm{R}\times \mathrm{R}\mathrm{e}\mathrm{d}+\mathrm{G}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{n}\times \mathrm{B}\mathrm{l}\mathrm{u}\mathrm{e}}}}$	[1]
23	Brightness Index (BI1)	$\mathrm{B}\mathrm{I}1=\sqrt{{\mathrm{R}\mathrm{e}\mathrm{d}}^2+{\mathrm{N}\mathrm{I}\mathrm{R}}^2}$	[78]
24	Brightness Index (BI2)	$\mathrm{B}\mathrm{I}2=\sqrt{{{\mathrm{G}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{n}}^2+\mathrm{R}\mathrm{e}\mathrm{d}}^2+{\mathrm{N}\mathrm{I}\mathrm{R}}^2}$	[75]

* Where Blue = B2; Green = B3; Red = B4; NIR = B8; SWIR1 = B11; SWIR2 = B12.

,

Table 3. Sentinel-1 features.

N	Index/Acronym	Formula * or Description	Reference
	Sentinel-1 polarization (dB)	VV, VH
1	Soil texture index 1 (IT1)	$\mathrm{I}\mathrm{T}1=\mathrm{V}\mathrm{V}+\mathrm{V}\mathrm{H}$	[25]
2	Soil texture index 1 (IT2)	$\mathrm{I}\mathrm{T}2={\mathrm{V}\mathrm{V}}^2+\mathrm{V}\mathrm{H}$
3	Soil texture index 1 (IT3)	$\mathrm{I}\mathrm{T}3={\mathrm{V}\mathrm{H}}^2-\mathrm{V}\mathrm{V}$
4	Soil texture index 1 (IT4)	$\mathrm{I}\mathrm{T}4={\mathrm{V}\mathrm{V}}^2+{\mathrm{V}\mathrm{H}}^2$
5	Soil texture index 1 (IT5)	$\mathrm{I}\mathrm{T}5={\displaystyle \frac{{\mathrm{V}\mathrm{V}}^2+{\mathrm{V}\mathrm{H}}^2}{\mathrm{V}\mathrm{H}}}$

* Where VV = Vertical transmission/Vertical reception; VH = Vertical transmission/horizontal reception.

and

Table 4. Topographic features.

N	Index/Acronym	Formula *	Reference
	DEM	$\mathrm{z}=\mathrm{f}(\mathrm{x},\mathrm{y})$	[79]
1	Slope	$\mathrm{G}=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\sqrt{{\mathrm{p}}^2+{\mathrm{q}}^2}$	[62]
2	Aspect	$\begin{array}{c}\mathrm{A}=-90\left[1-\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left(\mathrm{q}\right)\right]\left(1-\left|\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left(\mathrm{p}\right)\right|\right)+180\left[1+\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left(\mathrm{p}\right)\right]\\ -{\displaystyle \frac{180}{\mathsf{\pi }}}\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left(\mathrm{p}\right)\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{o}\mathrm{s}\left({\displaystyle \frac{-\mathrm{q}}{\sqrt{{\mathrm{p}}^2+{\mathrm{q}}^2}}}\right)\end{array}$$\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left(\mathrm{x}\right)=\left\{\begin{array}{c}1\mathrm{for}\mathrm{x}>0\\ 0\mathrm{for}\mathrm{x}=0\\ -1\mathrm{for}\mathrm{x}<0\end{array}\right.$
3	Hillshade	${\mathrm{H}}_{\mathrm{L}}=255(\cos \left(\mathrm{Z}\mathrm{e}\right)\times \cos \left(\mathrm{G}\right)+\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{Z}\mathrm{e})\times \mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{G})\times \mathrm{c}\mathrm{o}\mathrm{s}(\mathrm{A}\mathrm{z}-\mathrm{A})$
4	Northness	${\mathrm{A}}_{\mathrm{N}}=\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{A}$
5	Eastness	${\mathrm{A}}_{\mathrm{E}}=\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{A}$
6	Horizontal Curvature	${\mathrm{k}}_{\mathrm{h}}={\displaystyle \frac{{\mathrm{q}}^2\mathrm{r}-2\mathrm{p}\mathrm{q}\mathrm{s}+{\mathrm{p}}^2\mathrm{t}}{({\mathrm{p}}^2+{\mathrm{q}}^2)\sqrt{1+{\mathrm{p}}^2+{\mathrm{q}}^2}}}$
7	Vertical Curvature	${\mathrm{k}}_{\mathrm{v}}={\displaystyle \frac{{\mathrm{p}}^2\mathrm{r}-2\mathrm{p}\mathrm{q}\mathrm{s}+{\mathrm{q}}^2\mathrm{t}}{({\mathrm{p}}^2+{\mathrm{q}}^2)\sqrt{{\left(1+{\mathrm{p}}^2+{\mathrm{q}}^2\right)}^3}}}$
8	Mean Curvature	$\mathrm{H}={\displaystyle \frac{(1+{\mathrm{q}}^2)\mathrm{r}-2\mathrm{p}\mathrm{q}\mathrm{s}+(1+{\mathrm{p}}^2)\mathrm{t}}{2\sqrt{{\left(1+{\mathrm{p}}^2+{\mathrm{q}}^2\right)}^3}}}$
9	Gaussian Curvature	$\mathrm{K}={\displaystyle \frac{\mathrm{r}\mathrm{t}-{\mathrm{s}}^2}{{\left(1+{\mathrm{p}}^2+{\mathrm{q}}^2\right)}^2}}$
10	Minimal Curvature	${\mathrm{K}}_{\mathrm{m}\mathrm{i}\mathrm{n}}=\mathrm{H}-\sqrt{\left({\mathrm{H}}^2-\mathrm{K}\right)}$
11	Maximal Curvature	${\mathrm{K}}_{\mathrm{m}\mathrm{a}\mathrm{x}}=\mathrm{H}+\sqrt{\left({\mathrm{H}}^2-\mathrm{K}\right)}$
12	Shape Index	$\mathrm{S}\mathrm{I}={\displaystyle \frac{2}{\mathsf{\pi }}}\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}{\displaystyle \frac{\mathrm{H}}{\left({\mathrm{H}}^2-\mathrm{K}\right)}}$

* Where z = elevation (meters); Ze = Zenit (rad); Az = Azimuth (rad); r, t, s, p y q are the partial derivatives obtained by operating on linear distances calculated from the Haversine formula [62].

). The topographic indices were computed using the “tagee” library [62]. Those indices were selected as they have shown strong correlation with soil salinity [1,14,25,63,64,65].

3.2. Machine Learning Modelling Set-Up

To identify the potential of S1 and S2 features for soil salinity prediction, S1 and S2 features are considered explanatory features (i) separately to assess their respective potential and (ii) together to assess their complementarity. Finally, S1, S2, and topographic features are assessed together to highlight the potential benefit brought by topographic features. In this way, four scenarios are considered for each of the four ML models (i.e., RF, GT, GB, and XGB) (Figure 2).

Multi-collinearity is a well-known problem occurring in ML to decrease models’ robustness due to redundancy in selected independent features. To minimize such effects, the Recursive Feature Elimination Cross Validation algorithm (RFEcv) is used to find the subset of independent features with the best score [80]. As a model-specific feature selection method (i.e., wrapping method), the RFEcv is run independently for each model. In this approach, the model is run iteratively, removing one redundant feature at a time until the model performance decreases the least. In this process, the k-fold cross validation method (10-fold) with RMSE as the objective function is used. Then, the Variance Inflation Factor (VIF) is applied to the subset of independent features selected by the RFEcv to reduce the multi-collinearity even further. In this process, only the independent features with a VIF inferior to 10 are selected [81].

To assess the potential benefits of the feature selection process, each machine learning model (i.e., RF, GT, GB, and XGB) was trained and validated with (i) the optimized independent feature subset (i.e., obtained after the successive application of RFEcv and VIF) and (ii) with all independent features (i.e., without features selection). In this process, 70% and 30% of the learning database were used for the training and validation steps, respectively. The models’ reliability was evaluated with the coefficient of determination (R2, Equation (1)) and the Root Mean Square Error (RMSE, Equation (2)).

$${\mathrm{R}}^2={\displaystyle \frac{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{P}}_{\mathrm{i}}-{\mathrm{O}}_{\mathrm{i}}\right)}^2}{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{O}}_{\mathrm{i}}-{\mathrm{O}}_{\mathrm{i}}\right)}^2}}$$

(1)

$$\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt{{\displaystyle \frac{1}{\mathrm{n}}}{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{P}}_{\mathrm{i}}-{\mathrm{O}}_{\mathrm{i}}\right)}^2}$$

(2)

where n represents the number of samples and P_i and O_i represent the predicted and observed values of salinity content at the site i, respectively.

4. Results

4.1. Features Selection and Importance

Figure 3 shows the feature importance of the selected independent variables for each ML model and scenario. For scenario-1 (R), only three out of the five texture indices (ITs) were selected for modeling (IT2, IT3, and IT4). According to their equations (Table 4), all ITs are very similar, and therefore, their combination is not expected to bring additional benefits, which can explain why only a few remain after the feature selection process.

For scenario-2 (O), the process of feature selection significantly decreases the quantity of independent variables. The RF and GB models only select six out of twenty-six independent features, while the XGB and DT models selected only five out of twenty-six independent features. As for scenario-1, the indices corresponding to scenario-2 were similar in terms of the considered bands. Half of the indices was based on a combination of NIR and/or Red and/or SWIR bands (Table 2). The MSI2 index is the only one based on the three SWIR, NIR, and Red bands and thus was the most important feature. Therefore, MSI2 was seen as an effective way to synthetize the information brought by these bands. The NDWI is the second most important for all models. This index includes other bands’ information (i.e., Green) than the Red, NIR, and SWIR bands. No indices based on the Blue or Vegetation Red Edges (B5, B6, B7, and B8a) bands were selected, which shows that soil salinity content is not sensitive to the spectral range of these bands.

For scenario-3, all selected S2 indices (excepting SI2 for XGB model) are based on Red, NIR, SWIR, and Green bands, as for scenario-2. When combining S1 and S2 products, the most important indices come from S2 (i.e., MSI2 and NDWI), with the indices’ importance three times higher than those of the selected S1 features. This shows that in plowed lands, soil salt content is more sensitive to S2 than S1.

For scenario-4, no S1 indices remain, which confirms the low sensitivity of soil salt contents to S1 (as highlighted for scenario-3). The selected features are optical (MSI2 and NDWI) and topographic (Elevation and Slope). Elevation and slope were related to the identification of sites where water accumulates, evaporates, and leaves salts in the soil. Furthermore, lower areas are more exposed to groundwater rise and evaporation that will concentrate soil salts. These processes involve soil humidity dynamics that can explain the high contribution of MSI2, NDWI, and NDVI indices. Indeed, these indices are based on the NIR band (Table 2), which allows for the identification of soil moisture and was previously used across the considered region to separate dry soil from water bodies [46].

In summary, the selected variables show that humidity and drainage pattern play an important role in the soil salinization process.

4.2. Soil Salinity Estimates Prediction

Table 5. Performance of ML models with (SR) and without (TT) feature selection.

Model	Metric	Scenario-1 (R)		Scenario-2 (O)		Scenario-3 (R+O)		Scenario-4 (R+O+T)
		TT	SR	TT	SR	TT	SR	TT	SR
RF	R2	0.21	0.19	0.57	0.62	0.58	0.63	0.66	0.75 **
	RMSE (µS∙cm−1)	5646	2004	4143	3617	3776	3550	3820	2230 **
XGB	R2	0.03	0.07	0.64	0.69	0.63	0.66	0.76	0.76 **
	RMSE (µS∙cm−1)	7151	5247	3652	3336	3708	3773	3042	2890 **
DT	R2	0.07	−0.02	0.54	0.66	0.55	0.56	0.64	0.73 **
	RMSE (µS∙cm−1)	6638	6643	4832	3252	4610	3398	4206	2658
GB	R2	0.16	0.10	0.56	0.62	0.64	0.73	0.70	0.74 **
	RMSE (µS∙cm−1)	4513	6212	3897	3604	3507	3026	3421	3027 **

Where TT = All features; SR = selected features with RFEcv and VIF. ** Performance obtained without radar variables (R).

shows the metrics obtained in the validation stage for the models with and without feature selection.

The feature selection considerably improves the models’ reliability: for most scenarios, feature selection (RFEcv + VIF) increases (decreases) R² (RMSE). For the XGB model and scenario-1 (scenario-3), the improvement is limited to RMSE (R²). For example, when considering scenario-4, RMSE is reduced by 42%, 5%, 37%, and 12% (RF, XGB, DT, and GB models, respectively) by using feature selection.

All models provided much less reliable soil salinity estimates for scenario-1 than for scenario-2. Indeed, improvements in R² values of 226%, 886%, 3200%, and 520% were observed for RF, XGB, DT, and GB for scenario-2 (O) with respect to scenario-1 (R) (after the feature selection). This shows that estimates of the soil salinity of plowed agricultural lands are more sensitive to optical than radar indices, confirming the lower contribution of radar data highlighted in Figure 3c.

Finally, all models provide the most reliable soil salinity estimates for scenario-4 (O+R+T), indicating that the inclusion (exclusion) of topography (radar) indices lead to significant improvements in the model’s reliability (Figure 3d). An increase of 19%, 15%, 30%, and 1% is observed for RF, XGB, DT, and GB models for R² when comparing scenario-3 and -4 (Table 5).

4.3. Soil Salinity Mapping

Figure 4 shows soil salinity maps derived from the models based on the most reliable scenario for October 2023 (i.e., scenario-4).

Although all models produced similar metrics in the validation stage (Table 5), substantial differences are observed in the salinity spatial distribution predicted by the models. Taking as example an area near Lake Poopó (Figure 4), the XGB model yields more homogenous soil salinity mapping with very high salt concentrations (i.e., EC from 8000 to 15,000 µS/cm, Figure 4). On the contrary, the DT shows abrupt changes from high and extreme soil salinities in most of the area to no saline condition (i.e., <750 µS/cm) in the northwestern corner.

RF and GB models report a more balanced distribution with dominant high to very high salinity conditions. Both models agree with the higher salinity content in the southwestern part. However, the soil samples have a better agreement with the RF estimates than with the GB ones.

5. Discussion

Topography plays a fundamental role in soil formation and is therefore used as a key feature to estimate different soil properties, including salinity [1]. Actually, topography provides valuable information to identify regions where water is prone to accumulate (and therefore evaporate) and leave salt precipitates in the soils. Moreover, topography allows for the identification of depression zones which are prone to be closer to the groundwater table and are therefore more exposed to groundwater upraise which leaves salt precipitates in the soil [7]. In this context, numerous studies have already reported on the benefits of including topography indices as an independent feature in ML models for soil salinity mapping [27,82,83]. These findings are consistent with the results obtained in the current investigation. Indeed, a significant improvement in the soil salinity estimates is observed once topographic indices are included as independent features (i.e., scenario-4). By identifying topography as a key factor controlling soil salinity distribution, more reliable models could be generated with a more precise DEM. A study led in the United Kingdom shows that the use of a precise DEM derived from drones considerably improves the ML reliability for soil organic carbon predictions in plow lands [84]. As micro-topography inside plowed areas has already been identified as a key factor in soil salinity spatial distribution [85,86], such a DEM should considerably improve the model’s reliability.

Additionally, regions where water accumulates (and therefore where salts are present) are prone to higher soil humidity, which can be detected with spectral indices based on SWIR and NIR bands [87,88]. In this context, the MSI2 and NDWI indices are among the most relevant features for all scenarios and models. At an order of magnitude lower in relevance relative to the previous optical indices, the NDVI was also identified for most models, with this index being more relevant when monitoring saline soils with crops [89]. Numerous studies reported that S1 gathered valuable information on soil humidity estimates [90,91,92,93]. Based on these studies, it was expected that S1-derived features should improve soil salinity modeling. However, for all models, scenario-1 (only S1) reaches the lowest performance and S1 features’ contribution is insignificant to scenario-3 (S1 + S2) and -4 (S1 + S2 + T). This is due to the fact that the S1 C-band radar signal saturates in the case of plowed soils because plowed land roughness generally surpasses C-band wavelength (approximately 5.5 cm) [94,95].

Previous studies reporting on S1 sensitivity to soil humidity (and therefore soil salinity) were limited to non-disturbed soils. Indeed, more reliable soil salinity estimates were obtained in non-disturbed soil (low roughness) with S1 variables than with S2 features [26]. These contradictory observations regarding radar backscatter are sensitive to the dielectric constant [23] in non-disturbed and plowed soils, drawing attention to the necessity of providing knowledge of local soil conditions in order to identify the model limitations. Along this line, it is worth mentioning that S1 and S2 are not only sensitive to soil salinity, but also to soil organic carbon [96], nitrogen [97,98,99,100], texture [101,102], and pH [103,104,105]. In this regard, a model trained in a specific region may not apply to another region with different soil parameters (i.e., organic matter, nitrogen and texture). Thus, the proposed approach could be enhanced by training a specific model for every group of soils with a range of similar properties (chemical, physical, and biological) and that are spatially close (i.e., soil associations).

Finally, the model was developed over a specific time period (2022–2023), but trained for the main Altiplano soil type (Figure 1) where agriculture take place. As the soil type is not expected to change over time, the proposed model could be used to retrieve soil salinity maps along the Sentinel (-1 and -2) observation period (i.e., 2015-present). Therefore, this model could be used to assess soil salinity trends, and it links with agriculture practices (i.e., extensive quinoa cultivation in the Altiplano known as quinoa boom [106,107]). Such information could highlight the limit of current agricultural practice regarding its impact on soil salinity which might decrease soil productivity, therefore threatening both food security and socio-economic conditions. Such knowledge will promote investigations on regenerative agriculture effects [108] to define sustainable agriculture practices.

6. Conclusions

This study assesses the potential of Sentinel-2 and -1 (S2 and S1) products alone, together, and combined with topographic information for the soil salinity mapping of agricultural plowed lands. For this task, four machine learning (ML) models were used (Decision Tree—DT, random forest—RF, Gradient Boosting—GB, Extreme Gradient Boosting—XGB). The main results can be summarized as follows:

• Plowed land soil salinity estimates cannot be retrieved from S1 as the C-band wavelength is inferior to plowed land roughness, thus leading to the saturation of the signal.
• Accurate soil salinity estimates can be retrieved from S2 features alone (R² > 0.6), but the addition of topography features considerably improve the model’s reliability (R² > 0.7).
• The most relevant independent variables are the Moisture Saturation Index 2, Normalized Difference Water Index, and Elevation. These variables are related to the identification of water accumulation spots where salt is also expected to accumulate.
• For all ML models, the feature selection process significantly improves the model’s reliability, up to 13.6% in the best case.
• The four ML models reach similar global performance. However, high discrepancy is observed when using the models for the mapping of soil salinity in plowed lands in the studied region.

This study highlights the possibility of model soil salinity on plowed land in the Bolivian Altiplano, integrating data from remote sensing and machine learning algorithms. This information could be retrieved before the sowing period in order to provide information to minimize soil salinity risks with appropriate drainage measures. Furthermore, the higher spatial resolution of the elevation/relief information than that provided by the global elevation model (i.e., SRTM) should enhance the micro-topography representation on the agriculture plot scale. Such information would help to identify water accumulation zones inside the agriculture plot and therefore improve soil salinity modeling.

Finally, this study highlights remote sensing potential to assess soil salinization sensitivity to different factors (i.e., climate, irrigation, soil management). This knowledge could support the formulation of modern policies towards the sustainable management of agricultural lands affected by salinity under an integrative and preventive approach to minimize soil degradation and/or favor their rehabilitation.