Adapting Prediction Models to Bare Soil Fractional Cover for Extending Topsoil Clay Content Mapping Based on AVIRIS-NG Hyperspectral Data

Abstract

The deployment of remote sensing platforms has facilitated the mapping of soil properties to a great extent. However, the accuracy of these soil property estimates is compromised by the presence of non-soil cover, which introduces interference with the acquired reflectance spectra over pixels. Therefore, current soil property estimation by remote sensing is limited to bare soil pixels, which are identified based on spectral indices of vegetation. Our study proposes a composite mapping approach to extend the soil properties mapping beyond bare soil pixels, associated with an uncertainty map. The proposed approach first classified the pixels based on their bare soil fractional cover by spectral unmixing. Then, a specific regression model was built and applied to each bare soil fractional cover class to estimate clay content. Finally, the clay content maps created for each bare soil fractional cover class were mosaicked to create a composite map of clay content estimations. A bootstrap procedure was used to estimate the standard deviation of clay content predictions per bare soil fractional cover dataset, which represented the uncertainty of estimations. This study used a hyperspectral image acquired by the Airborne Visible/Infrared Imaging Spectrometer-Next Generation (AVIRIS-NG) sensor over cultivated fields in South India. The proposed approach provided modest performances in prediction ($R_{val}^2$ ranging from 0.53 to 0.63) depending on the bare soil fractional cover class and showed a correct spatial pattern, regardless of the bare soil fraction classes. The model’s performance was observed to increase with the adoption of higher bare soil fractional cover thresholds. The mapped area ranged from 10.4% for pixels with bare soil fractional cover >0.7 to 52.7% for pixels with bare soil fractional cover >0.3. The approach thus extended the mapped surface by 42.4%, while maintaining acceptable prediction performances. Finally, the proposed approach could be adopted to extend the mapping capability of planned and current hyperspectral satellite missions.

1. Introduction

The achievement of sustainable land management requires a comprehensive understanding of earth’s soil at high spatio-temporal resolutions. The information on several soil properties such as moisture content, nutrient availability, soil organic carbon (SOC) content and water holding capacity are crucial for policymaking and land resource management for ensuring global food security [1]. Soil texture, which refers to the relative proportions of sand, silt, and clay particles in a soil mass, is widely considered as a significant soil characteristic. Soil texture impacts on plant growth, crop yield, soil water infiltration and retention capacities, and absorption of nutrients by plants [2] and its mapping is crucial over large areas for cropland management. However, conventional methods of soil analysis are expensive and time consuming and are consequently characterized by low spatial resolutions that fail to adapt the spatial variation exhibited by soil texture in a landscape.

To address this limitation, visible–near-infrared and short-wave infrared (VNIR–SWIR, 400 nm to 2500 nm) imaging spectroscopy has been extensively used for mapping several soil properties, including soil texture [3]. Several studies have estimated the clay content of the soil from its VNIR–SWIR reflectance spectra based on a diagnostic absorption signature near 2200 nm to 2300 nm induced by the vibrations associated with the OH bond and the OH-Al-OH bonds of kaolinite, montmorillonite, and illite minerals [4]. These studies typically employ continuum removal methods [5,6,7], multivariate regression models [8,9,10], or machine learning techniques [11,12,13] for estimating the clay content from VNIR–SWIR reflectance spectra.

A major limitation of soil spectroscopy from remote sensing platforms is that the mapping of soil properties, including clay content mapping, is limited to bare soil surfaces. The presence of non-soil cover affects the reflectance spectra acquired over the pixels and introduces a spectral mixing effect, which further affects the soil property estimated from the spectra [14,15,16]. This limitation is often encountered, particularly over cultivated areas that are typically covered by photosynthetic or non-photosynthetic vegetation during most of the year.

Some studies have attempted to address this issue by employing spectral unmixing to produce new spectra that are reasonably well representative of the pixel components. This approach, which takes advantage of data with high spectral resolution, was developed on hyperspectral VNIR–SWIR airborne imagery. Bartholomeus et al. [16] used a residual spectral unmixing approach on Airborne Hyperspectral Scanner-160 data to extract pure soil spectra over partially vegetated pixels and predict SOC content. Ouerghemmi et al. [17,18] used a blind source separation technique to extract pure soil spectra over partially vegetated pixels in a Hymap image and estimate their clay contents. Franceschini et al. [19] used a progressive spectral dataset selection based on unmixed soil fractions for soil property estimation using airborne AisaDUAL data.

With the availability of multi-temporal data from multispectral sensors such as Sentinel-2 and Landsat, several recent studies extracted bare soil pixels from multiple images over a common study site and combined them to generate a new composite image with a larger number of bare soil pixels [20,21,22,23,24]. These multi-temporal images could either be from the same sensor or from different sensors, over which an empirical selection of a spectral index threshold is applied to differentiate bare soil cover from non-soil cover [23]. Some spectral indices are commonly used to this extent, such as the Normalized Difference Vegetation Index (NDVI) for identifying green vegetation [21,22,25], Normalized Burn Ratio 2 (NBR2) for identifying dry vegetation [21,26,27], and the Bare Soil Index (BSI) for identifying bare soil [20,23]. Dematte et al. [21] proposed to use a multi temporal data mining approach to generate a Synthetic Soil Image (SYSI) by using the median of the spectral reflectance of bare soil pixels based on the Geospatial Soil Sensing System (GEOS3). This image provides a synthetic spectral reflectance of soils, and the approach has been adopted in several similar studies [13,28,29,30,31]. These approaches, creating bare soil composite images and synthetic imagery, enabled extended soil property mapping, including clay content, SOC, and cation exchange capacity [24,30,31]. Some studies have also explored the effect of factors such as the index used in bare soil selection [32], the length of the compositing period [33], and the pixelwise mosaicking approach used for generating the composite [34] on the prediction accuracy of estimated soil properties.

In this context, while the recent literature uses a temporal perspective to extend the soil properties mapping, our study proposes a new approach from a spatial perspective to (i) extend the soil properties mapping beyond bare soil pixels based on one hyperspectral VNIR–SWIR airborne image and (ii) provide an uncertainty map of the soil property estimations. This approach first classifies the pixels based on their bare soil fractional cover from the spectral unmixing process and then applies a specific regression model to each bare soil fractional cover class to estimate clay content along with the associated uncertainty in prediction. The selected study area is in the southern part of the Indian subcontinent, where the crop fields are typically small-sized, scattered, and farmed with diverse crop varieties.

2. Materials

2.1. Study Area

The study area covers 300 km² in Gundlupet Taluk of the Indian state of Karnataka (Figure 1a). A large part of the Berambadi catchment which belongs to the Kabini Critical Zone Observatory (Assimilation of Multi-satellite data at Berambadi watershed for Hydrology And land Surface (AMBHAS), Service d’Observation—Bassins Versants Expérimentaux Tropicaux (SO-BVET), Tomer et al. [35]; Sekhar et al. [36]) and is a part of the Critical Zone Observatories: Research and Application network (OZCAR; Gaillardet et al. [37]), falls within the study area. The climate of the region is tropical sub-humid, and it receives rainfall of 800 mm annually. The study area is covered by both cultivated fields (where multi-cropping and intercropping are practiced) and forested lands. The size of the crop fields over the cultivated part varies from 0.01 to 9.3 hectares, with an average size of 1.2 hectares [38]. These are farmed with a large variety of seasonal crops, and farmers opt for irrigation using drip techniques, sprinklers, or furrow channels, depending on the availability of water and type of crop [38]. The cropping system in the region is regulated by three seasons: (1) kharif or the south-western monsoon season, from June to September, when the majority of the fields are cultivated; (2) rabi or the north-eastern monsoon, from October to January, when most of the fields are irrigated; and (3) the hot and dry summer season, from February to May, when only a few of the fields are cultivated where irrigation is viable [38].

The terrain of the study area comprises uplands with gentle slopes and lowlands that are nearly level, in addition to hills and hillslopes. The soils of the study area are mostly developed on granitic gneiss bedrock [39]. The region is characterized by red soil over uplands and hillslopes (Ferralsols and Chromic Luvisols), black soil over valley bottoms (Vertisols and Vertic intergrades), and rocky/weathered soil [39].

2.2. Soil Sample Collection and Clay Content Measurements

A total of 272 soil samples were collected over the study area during the rabi cropping season, with 180 samples in November 2019 and 92 in January 2020 (Figure 1c). As clay content is a perennial property in topsoil horizons, it was assumed that the clay content analyzed from soil samples collected in 2019 and 2020 and the clay content over the field in 2016 (during the AVIRIS-NG campaign) remained unchanged as previously assumed in many studies [5]. During the period of data collection, some of the crop fields in the study area were cultivated, harvested, or left bare. The strategy adopted for soil sampling was to characterize the heterogeneity of soils, crops, and landforms within the study area using stratified sampling based on both soils, landuse and geological maps, and human expertise. Since red soils located over uplands and hillslopes make up the majority of the study area [40], more samples (195 samples) were taken over them than over black soils, which are located over valley bottoms (77 samples). In order to represent the main landforms of the region, the sampling locations were dispersed throughout nearly level lowlands, moderately sloping uplands, and valley bottoms. The soil samples were collected from a depth of 0 to 15 cm. The samples thus collected were perceived to characterize the cultivated soil horizon as the process of tilling the fields prior to planting homogenizes the topsoil [41]. Each sample was a mixture of five sub-samples collected at the center and the four vertices, respectively, of a 10 m × 10 m square. These samples were assumed to be representative of an area of 5 m × 5 m, as the variation in clay content remains low within small areas [42,43]. A Garmin GPS instrument was used to record the location of each of the sampled plots.

The clay content of the soil samples, i.e., the percentage of particles with diameter <0.002 mm as per the USDA textural classification, was measured using Robinson’s pipette method [44,45]. The samples were found to have clay contents that ranged from 5.8% to 62.2%, with a mean of 26.6%, a standard deviation of 13%, and a coefficient of variation of 48.8%. This soil dataset was previously used by George et al. [46].

2.3. AVIRIS-NG Hyperspectral Data

This work used a hyperspectral image acquired by the Airborne Visible/Infrared Imaging Spectrometer—Next Generation (AVIRIS-NG) imaging sensor over the study area on 10th of January 2016 from IST 12:00 to 13:45 h (Figure 1c). This acquisition was a part of the joint Indian Space Research Organization (ISRO)—National Aeronautics and Space Administration (NASA) hyperspectral science campaign [47] and corresponds to the image over the site ID 99, Muddur, Karnataka. The image was acquired under clear sky conditions. Since the study area had received rain 20 days prior to the date of image acquisition, the land surface can be considered dry at the time. The AVIRIS-NG data used in this study was obtained from the ISRO VEDAS portal (https://vedas.sac.gov.in/, accessed on 1 July 2021).

The AVIRIS-NG imaging spectrometer measures the reflected radiance in wavelengths in the 380–2510 nm range, at every 5 nm interval, thus resulting in a total of 425 spectral bands, with a ground sampling distance of 3.8 m for a flight altitude of 4.6 km and a high signal-to-noise ratio (SNR) > 2000 at 600 nm and > 1000 at 2200 nm [47]. The sensor had a field-of-view (FOV) of 34° and an instantaneous FOV (IFOV) of 1 milliradian [48,49]. The data had been ortho-rectified and atmospherically corrected using the ATmosphere REMoval algorithm (ATREM, Gao et al. [50]) to produce level 2 surface reflectance output, with values from 0 to 1 [47]. After eliminating water absorption bands and noisy bands, 260 spectral bands were retained for the modeling.

The images corresponding to nine flight lines were mosaicked to encompass the study area (Figure 1c). Water bodies in the study area were masked using an expert-calibrated threshold: pixels with a reflectance of less than 8% at 1663.67 nm (Band 258) were removed [43]. In addition, urban areas, forested areas towards the western end of the study area, and uncultivated hills were identified by visual inspection and field knowledge and were masked, and only the cultivated areas were kept for this work. This AVIRIS-NG data was previously described and used in George et al. [46].

3. Methods

The study proposed an approach to extend the soil properties mapping over non-soil pixels and provide an associated uncertainty map. The proposed approach first estimated the bare soil fractional cover of each pixel by spectral unmixing and then classified the pixels according to their bare soil fractional cover (Figure 2a). Specific regression models were then built for datasets with varying bare soil fractional cover thresholds to estimate the clay content associated with each pixel in the dataset (Figure 2b). A bootstrap procedure was used to estimate the standard deviation of clay content predictions per bare soil fractional cover dataset, which represented the uncertainty of estimations. Each regression model was then applied to the pixels belonging to the corresponding bare soil fractional cover class, providing a clay content map for the class (Figure 2c). Finally, the clay content maps created for each bare soil fractional cover class were mosaicked to create a composite map of clay content estimations, and the standard deviation of clay content predictions obtained through bootstrapping was used to create the associated uncertainty map (Figure 2d).

3.1. Spectral Unmixing for Bare Soil Fractional Cover Estimation

The spectral unmixing of this AVIRIS-NG image has been previously conducted by George et al. [46] (Figure 2a). Spectral unmixing is a process of deconvoluting the reflectance spectrum of a mixed pixel into its constituent pure component (also called ‘endmember’) spectra and the corresponding fractional covers of fractions (also called ‘abundances’) [51]. Based on the assumption that reflectance spectra of endmembers are linearly independent, the fully constrained least squares unmixing (FCLS) algorithm [52] was used to estimate the abundance of each endmember within each pixel [53]. The linear mixture model for a pixel with n spectral bands and p endmembers can be represented as:

$$y_i=M\alpha +ϵ$$

(1)

where $y_i$ is $n\times 1$ spectral signature of the pixel i, $M=\left[M_1, \dots , M_p\right]$ is an $n\times p$ matrix of endmember spectral reflectances, $\alpha$ is the $p\times 1$ abundance vector, and $ϵ$ is the measurement error [52]. The endmember abundances are then estimated by adopting the least squares technique to minimize the error $\Vert y-M\alpha \Vert$ as follows:

$$\alpha ={\left(M^{T}M\right)}^{-1}{M}^{T}y$$

(2)

Further, constraints were imposed for the estimated endmember abundances within a pixel to add up to unity and to be non-negative [52]:

$${\alpha }_j\ge 0 : \forall p$$

(3)

$${\sum }_1^p{\alpha }_j=1$$

(4)

Since the study area consisted of crop fields with actively growing vegetation or harvested fields that were left bare or had dry crop residue cover, George et al. [46] assumed that the AVIRIS-NG spectra could be modeled based on the three following endmembers: soil, photosynthetic vegetation (PV), and non-photosynthetic vegetation (NPV). The endmembers were derived directly from the AVIRIS-NG image as proposed by Plaza et al. [54], based on visual examination of the image, ground truth knowledge, and examining the reflectance spectra. Eight pure pixels belonging to each of the three classes (soil, PV, and NPV) were identified, and their respective average reflectances were considered the endmember spectra for unmixing. Finally, the map of the bare soil fractional cover, with values in pixels varying from 0 to 1, where 0 indicates no bare soil coverage and 1 indicates full bare soil coverage within a pixel, was produced over the study area [46]. The pixels in the bare soil fractional cover map were subsequently classified into nine classes C_p ranging from 0.3 to 1 (with p ranging from 1 to 9, Figure 2a), providing a map called SFClassMap (Supplementary Material Figure S1). The pixels with a bare soil fractional cover below 0.3 were not considered for further analysis as they corresponded to the lowest bare soil coverage, which implied the densest vegetation coverage (NPV or PV), and were subsequently masked in the SFClassMap.

3.2. PLSR Model

The Partial Least Squares Regression (PLSR) model was the quantitative statistical analysis technique used for estimating the clay content (response variable) from the AVIRIS-NG spectra (predictor variable). PLSR extracts latent predictor variables, which account for the maximum amount of variation in the response variables [55]. This technique is commonly used for soil property estimation through spectroscopy because of its ability to analyze data with noisy and collinear variables [56,57,58].

Before establishing the PLSR models, the AVIRIS-NG reflectance spectra were converted to “pseudo-absorbance” (log [1/reflectance]) to reduce nonlinearities [56], and noise reduction was achieved using the standard pretreatments: a second-degree Savitzky–Golay filter with a span of five was applied across the discontinued AVIRIS-NG spectra [59] and a mean centering and variance scaling were applied. Spectra with a Mahalanobis distance [60] greater than 3.5 were considered outliers and removed from the calibration set. The selection of the number of latent variables (LVs) was done such that the root mean squared error of prediction for the calibration dataset estimated through leave-one-out cross validation was minimized.

3.3. Dataset Preparation for PLSR Training and Validation

For estimating clay content from the AVIRIS-NG image, nine datasets, denoted as DB_p (with p ranging from 1 to 9), were created. Each dataset DB_p consisted of AVIRIS-NG pixels meeting two criteria:

(1) The corresponding clay contents measured in the laboratory were available;

(2) The bare soil fractional cover in the pixel was up to a specified threshold T_p. This threshold T_p ranged from >0.3 to >0.7 in increments of 0.05 (Figure 2b).

A first set of PLSR models, denoted as M_p, was trained from subsets of DB_p using an equal number of samples (N_cal) in all calibration datasets. These PLSR models M_p were then validated from subsets of DB_p using an equal number of samples (N_val) in all validation datasets (where N_cal and N_val are different). N_cal and N_val were limited by the most restrictive DB_p dataset, i.e., DB₉, which corresponded to the bare soil fractional cover greater than 0.7 and consisted of 46 samples. By adopting an 80–20% split for calibration–validation datasets [61,62], each M_p model was calibrated on 36 samples (N_cal) and validated on 10 samples (N_val), belonging to the corresponding dataset DB_p. This equal number of samples in both the calibration (N_cal = 36) and validation (N_val = 10) datasets among the models allowed a robust comparison of the effect of bare soil fractional cover on the model performance.

A second set of PLSR models, denoted as M_{all_p}, were trained using the maximum number of available samples in the DB_p dataset, keeping aside N_val samples for validation. Similarly, we maintained consistent sample sizes N_val (where N_val = 10) in validation datasets across all M_{all_p} models to facilitate performance comparisons.

3.4. Bootstrap Process

To provide a robust performance evaluation and standard deviation of predictions, a bootstrap procedure was adopted [63] (Figure 2b). Following previous soil spectroscopic studies that had adopted bootstrapping procedures [64,65], 100 M_p models were built, each one denoted M_{p_i} (with i ranging from 1 to 100) (Figure 2b). For each of the M_{p_i} models, the selection of 36 samples (N_cal) for calibration and 10 samples (N_val) for validation followed stratified random sampling based on the clay contents from the corresponding DB_p dataset.

Similarly, 100 M_{all_p} models were built, each one denoted by M_{all_p_i} (with i ranging from 1 to 100) (Figure 2b). For each of the M_{all_p_i} models, the selection of 10 samples (N_val) for validation followed stratified random sampling based on the clay contents from the corresponding DB_p dataset. The remaining samples in the DB_p dataset were used for calibrating the M_{all_p_i} models.

It is to be noted that in both configurations, the models M_{all_9} and M₉ are the same, as they are calibrated and validated on the same number of samples.

3.5. Performance Evaluation

The performances of each of the PLSR models were evaluated using the coefficient of determination of validation ($R_{val}^2$) and the root mean square error in prediction (RMSEP) (Figure 2b), which were calculated as follows:

$$R_{val}^2=1-\frac{{\sum }_{j=1}^n{\left({\tilde{y}}_j-y_j\right)}^2}{{\sum }_{j=1}^n{\left(y_j-\widehat{y}\right)}^2}$$

(5)

$$RMSEP=\sqrt{\frac{{\sum }_1^n{({\tilde{y}}_j-y_j)}^2}{n}}$$

(6)

where ${\tilde{y}}_j$ is the jth estimated clay content, $y_j$ is the jth measured clay content, $\widehat{y}$ is the average value of measured clay contents, and n is the number of samples. Further, the mean and standard deviation of the $R_{val}^2$ and RMSEP of the 100 M_{p_i} models from the bootstrap procedure were calculated. Similarly, the mean and standard deviation of the $R_{val}^2$ and RMSEP from the 100 M_{all_p_i} models from the bootstrap procedure were also calculated (Figure 2b).

3.6. Clay Content Composite Mapping

The model M_{all_p_i} trained using pixels with bare soil fractional cover greater than a threshold T_p, was applied to all pixels belonging to class C_p in the SFClassMap for estimating the corresponding clay contents (Figure 2c). Thus, 100 regression models M_{all_p_i} were applied to each pixel belonging to class C_p, resulting in 100 clay content estimations, whose mean value represented the final clay content estimated at the pixel (Figure 2c). A final clay content composite map named CompositeClayMap was then created by mosaicking the clay content estimations obtained per bare soil fractional cover class C_p (Figure 2d).

Likewise, the standard deviation of the 100 clay content estimations obtained from bootstrapping (M_{all_p_i}) was calculated, which represented the prediction uncertainty. The mosaicking of these standard deviation values generated the associated uncertainty map, named UncertaintyMap (Figure 2d).

4. Results

4.1. Bare Soil Fractional Cover Class Map

The SFClassMap generated from the bare soil fractional cover values estimated through FCLS unmixing covered 52.7% of the study area (Figure 3a; Supplementary Material Figure S1). The percentage of the study area covered by the bare soil fractional cover classes varied from 2.8% for class C₈ to 10.4% for class C₉ (

Table 1. Area covered by each of the bare soil fractional cover classes obtained by the FCLS unmixing.

Bare Soil Fractional Cover	Class (Cp)	Number of Available Soil Samples	Area (%)
0.30–0.35	C1	17	7.3
0.35–0.40	C2	18	7.4
0.40–0.45	C3	20	6.9
0.45–0.50	C4	14	5.9
0.50–0.55	C5	13	4.9
0.55–0.60	C6	14	3.9
0.60–0.65	C7	8	3.2
0.65–0.70	C8	17	2.8
0.70–1.00	C9	46	10.4
0.30–1.00		167	52.7

). The number of soil samples available for model ing per bare soil fractional cover class varied from 8 for C₇ (bare soil fractional cover from 0.6 to 0.65) to 46 for C₉ (bare soil fractional cover from 0.7 to 1) (Table 1). The number of soil samples N_DBp available for model ing per bare soil fractional cover threshold (T_p) varied from 46 for T₉ (threshold > 0.7) to 167 for T₁ (threshold > 0.3) (

Table 2. Details of the available soil samples associated with each bare soil fractional cover threshold T_p.

Dataset (DBp)	Bare Soil Fractional Cover Threshold (Tp)		Number of Available Samples (NDBp)	Minimum Clay Content (%)	Maximum Clay Content (%)	Mean Clay Content (%)	Standard Deviation of Clay Content (%)
DB1	T1	>0.30	167	5.81	59.61	26.21	12.72
DB2	T2	>0.35	150	5.83	59.61	25.70	12.53
DB3	T3	>0.40	132	5.83	59.61	26.10	12.61
DB4	T4	>0.45	112	5.83	59.61	26.70	12.50
DB5	T5	>0.50	98	5.83	59.61	27.12	12.43
DB6	T6	>0.55	85	6.96	58.91	26.60	11.47
DB7	T7	>0.60	71	6.96	49.04	26.01	11.16
DB8	T8	>0.65	63	6.96	49.04	25.48	10.99
DB9	T9	>0.70	46	6.96	47.30	24.91	9.88

). So while 272 soil samples were initially available (Section 2.2), only 167 were finally used, as the remaining 107 samples were collected over surfaces associated with bare soil fractional cover less than 0.3. Finally, the nine datasets DB_p (with p ranging from 1 to 9) had similar ranges and distributions of measured clay contents, regardless of the number of samples associated with them (Table 2).

4.2. Effect of Bare Soil Fractional Cover on Clay Content Estimation

The first set of regression models M_p built from the bootstrap procedure with the same number of calibration data (N_cal = 36) yielded modest to good performances (

Table 3. Performances of M_p models for clay content estimation.

Bare Soil Fractional Cover	Dataset (DBp)	PLSR Model (Mp)	Ncal	Nval	Validation (Mean ± Standard Deviation over 100 Iterations)
					${\mathit{R}}_{\mathit{v}\mathit{a}\mathit{l}}^2$	RMSEP (%)
>0.30	DB1	M1	36	10	0.46 ± 0.22	9.27 ± 2.00
>0.35	DB2	M2	36	10	0.47 ± 0.20	9.16 ± 2.23
>0.40	DB3	M3	36	10	0.49 ± 0.20	9.01 ± 1.90
>0.45	DB4	M4	36	10	0.49 ± 0.19	9.04 ± 1.95
>0.50	DB5	M5	36	10	0.50 ± 0.19	8.63 ± 1.69
>0.55	DB6	M6	36	10	0.50 ± 0.19	8.04 ± 1.90
>0.60	DB7	M7	36	10	0.52 ± 0.18	7.53 ± 1.51
>0.65	DB8	M8	36	10	0.59 ± 0.17	7.00 ± 1.41
>0.70	DB9	M9 *	36	10	0.63 ± 0.13	6.13 ± 1.09

* M₉ = M_{all_9}.

). The best performance was obtained for the model M₉, which considered the higher bare soil fractional cover (>0.7) (an $R_{val}^2$ of 0.63 and an RMSEP of 6.13% as the mean of 100 iterations), while the worst performance was obtained for the model M₁, which considered the lowest bare soil fractional cover (>0.3) (an $R_{val}^2$ of 0.46 and an RMSEP of 9.27% as the mean of 100 iterations) (Table 3). As expected, an increase in bare soil fractional cover in the data used for training the PLSR models involved an improvement in clay content estimation accuracy (Table 3). This increase in performance may be attributed to a better selection of bare soil pixels by the stricter threshold, which implies a better selection of soil spectra unaffected by vegetation.

The standard deviation of $R_{val}^2$ from the 100 iterations ranged from 0.13 for a bare soil fractional cover higher than 0.7 to 0.22 for a bare soil fractional cover higher than 0.3 (Table 3). The standard deviation of RMSEP from the 100 iterations ranged from 1.09% for a bare soil fractional cover exceeding 0.7 to 2.23% for a bare soil fractional cover higher than 0.35 (Table 3). Finally, an increase in the bare soil fractional cover threshold involved a decrease in the standard deviation of clay content predictions (Table 3).

The second set of regression models M_{all_p} was trained from the maximum number of available samples (N_{cal_all_p}) per bare soil fractional cover threshold, keeping aside N_val samples for validation (where N_val = 10). Therefore, the number of calibration samples (N_{cal_all_p}) ranged from 36 to 157 (

Table 4. Performances of M_{all_p} models for clay content estimation.

Bare Soil Fractional Cover	Dataset (DBp)	PLSR Model (Mall_p)	Ncal_all_p	Nval	Validation (Mean ± Standard Deviation over 100 Iterations)
					${\mathit{R}}_{\mathit{v}\mathit{a}\mathit{l}}^2$	RMSEP (%)
>0.30	DB1	Mall_1	157	10	0.53 ± 0.20	8.30 ± 1.96
>0.35	DB2	Mall_2	140	10	0.54 ± 0.21	8.20 ± 2.00
>0.40	DB3	Mall_3	122	10	0.54 ± 0.18	8.19 ± 1.78
>0.45	DB4	Mall_4	102	10	0.53 ± 0.19	8.19 ± 1.72
>0.50	DB5	Mall_5	88	10	0.53 ± 0.21	7.90 ± 1.86
>0.55	DB6	Mall_6	75	10	0.55 ± 0.17	7.26 ± 1.72
>0.60	DB7	Mall_7	61	10	0.59 ± 0.14	6.97 ± 1.29
>0.65	DB8	Mall_8	53	10	0.61 ± 0.16	6.50 ± 1.37
>0.70	DB9	Mall_9 *	36	10	0.63 ± 0.13	6.13 ± 1.09

* M_{all_9} = M₉.

). These models also yielded modest to good performances, with $R_{val}^2$ (as the mean of 100 iterations) from 0.53 (obtained for bare soil fractional cover >0.3, >0.45, and >0.5), to 0.63 (for bare soil fractional cover > 0.7) and RMSEP (as the mean of 100 iterations) from 6.13% (for bare soil fractional cover > 0.7) to 8.30% (for bare soil fractional cover > 0.3) (Table 4). As observed with the M_p models, the performance of M_{all_p} in clay content estimation improved with an increase in bare soil fractional cover.

For the M_{all_p} models, the standard deviation of $R_{val}^2$ from the 100 iterations ranged from 0.13 for bare soil fractional cover > 0.7 to 0.21 for bare soil fractional cover >0.35 and >0.5 (Table 4). The standard deviation RMSEP for the M_{all_p} models from the 100 iterations ranged from 1.09% for bare soil fractional cover > 0.7 to 2.00% for bare soil fractional cover > 0.35 (Table 4). Finally, an increase in bare soil fractional cover involved a decrease in the standard deviation of clay content estimation, as observed for the M_p models (Table 3 and Table 4).

The M_{all_p} models thus exhibited a similar trend in PLSR model performance for clay content estimation as the M_p models, i.e., an increase in $R_{val}^2$ and a decrease in RMSEP as the bare soil fractional cover increased. It was also observed that the M_{all_p} models exhibited better performances than the M_p models for the same p (Table 3 and Table 4). Thus, using a greater number of samples for training a PLSR model improved its performance in clay content estimation.

4.3. Clay Content Composite Map and Associated Uncertainties Map

To take advantage of the maximum number of soil samples in regression models and benefit from the best model performances, the M_{all_p_i} regression models were used for clay content mapping. The 100 regression models M_{all_p_i} were applied to generate 100 clay content estimations for each pixel in the AVIRIS-NG image, based on their bare soil fractional cover class C_p (Figure 2c). The mean of these 100 values was reported as the clay content at the pixel, and the mosaicking of these estimations produced the clay content composite map, CompositeClayMap (Supplementary Material Figure S2). This approach allowed mapping over 52.7% of the study area, whereas a single model M_{all_9} would have allowed mapping over only 10.4% of the study area (Table 1).

The spatial pattern of the CompositeClayMap was coherent with the soilscape trend of the region. High clay content values were estimated over the black soils corresponding to Vertisol and Vertic intergrades, found in the valley bottom around the riverbeds (with a mean clay content of around 40%), while lower clay contents were estimated over the red soils corresponding to Ferralsol and Chromic Luvisols, found in the uplands and hillslopes (with a mean clay content < 40%). The region towards the east of the Berambadi lake that is characterized by red soils was rightly identified with 10 to 30% of the estimated clay content (Figure 3b, pink ellipse). Similarly, the black soils present to the south of the lake were identified with higher clay contents of around 40 to 60% in the CompositeClayMap (Figure 3b, black ellipse). These observations were consistent with the pedological map and knowledge of the study area. The clay content composite map also compared well with the soil texture classification by Gomez et al. [66] and clay content mapping by George et al. [46] over the study area using Sentinel-2 and AVIRIS-NG images, respectively.

This coherent spatial pattern for clay content was observed from CompositeClayMap, regardless of the bare soil fractional cover classes (Figure 3a,b). Even the lower bare soil fractional cover classes (C₁ to C₆) with relatively higher RMSEP and lower $R_{val}^2$ were able to discriminate between the red and black soil in the study area. This implied that even over low bare soil fractional cover, the PLSR models were able to provide the correct spatial pattern of predicted clay values, which can be useful for supplementing the pedological knowledge at the regional scale.

The standard deviation of the predicted clay contents over 100 iterations, which represents the uncertainty of prediction in a pixel over the entire area, had a mean of 4.62% (Supplementary Material Figure S3). In these, 0.01% of the pixels spread across the study area without exhibiting any spatial pattern were associated with high uncertainties (>30%) and considered outliers and excluded from the analysis. Over the study area, pixels associated with high uncertainties (standard deviation from 15 to 30%) were located over the forested region to the north-west (SSupplementary Material Figure S1), where bare soil fractional cover was estimated to vary from 0.3 to 0.45, which could possibly be an overestimation of the bare soil fraction in this area. Over the cultivated area, both red and black soils (Figure 3b, pink and black ellipses, respectively) were associated with lower uncertainties (standard deviation < 10%), except over some specific agricultural fields that had high uncertainties (e.g., Figure 3c, purple ellipse), which indicated that clay estimations over these specific fields must be considered with caution.

The mean values of the uncertainty in prediction associated with each bare soil fractional cover class varied from 2.81% for C₉ to 5.78% for C₃, while the median varied from 2.54% for C₉ to 4.6% for C₃ and C₄ (Figure 4). Violin plots for each bare soil fractional cover class C_p were generated to compare the distribution of the uncertainties associated with the class (Figure 4). The higher bare soil fractional cover class C₉ was associated with a lower range of uncertainty as compared to the lower bare soil fractional cover classes (Figure 4). Finally, though there was a statistically significant difference between the uncertainties of bare soil fractional cover classes as determined by one-way ANOVA (F(8, 1.1 × 10⁷) = 188,787.14, p < 0.0001), there was no correlation between the bare soil fraction classes (Figure 3a) and the standard deviations (Figure 3c) of predictions over the entire mapped study area (R² of 0.08, p < 0.0001).

5. Discussion

5.1. Bare Soil Fractional Cover Estimation from Spectral Unmixing

The choice of a suitable collection of endmembers may have a significant impact on the result of spectral unmixing [54]. The endmembers used for unmixing can be selected either from a spectral library of known materials (acquired from a laboratory or field) or directly from the image [67]. In this study, the three endmembers (soil, PV, and NPV) were derived from the AVIRIS-NG image since the spectra of these image endmembers would have the same spatial and spectral resolutions and similar acquisition conditions as the unmixing input image [54].

Estimating the number of endmembers to be used for unmixing is an application dependent and subjective task [68]. Incorporating any of the available endmember variability reduction techniques, such as iterative mixture analysis, spectral feature selection, spectral transformations, etc., could account for the temporal and spatial variability between and among endmembers [69]. While the study area encompasses cultivated lands featuring a variety of crops at different growth stages, leading to intra-class spectral variability for vegetation endmembers (PV and NPV), we have simplified their representation by assigning each a single spectral signature. Similarly, despite the Berambadi catchment being characterized by both red and black soil, a singular soil endmember was chosen. It would be worth exploring the inclusion of two soil endmembers (black soil and red soil) and multiple PV and NPV endmembers in the unmixing process to better capture the intra-class spectral variability across this complex and heterogeneous study area.

5.2. Model Performance over the Highest Bare Soil Fractional Cover

The PLSR models M₉, calibrated using samples corresponding to the highest bare soil fractional cover class (>0.7), provided accurate predictions for topsoil clay content (Table 3 and Table 4). These model performances fell within the range reported in the literature dedicated to bare soil pixels (R² from 0.53 to 0.75). It slightly surpassed the performance for clay content estimation ($R_{val}^2$ of 0.55) achieved by Casa et al. [8] using PLSR modeling based on airborne hyperspectral MIVIS data. It also surpassed the performance for clay content estimation ($R_{val}^2$ from 0.39 to 0.57) achieved by Castaldi et al. [9] using PLSR modeling and resampled spectra from various hyperspectral satellite sensors, including EO-1 Hyperion, PRISMA, HyspIRI, and EnMAP, with the addition of noise and atmospheric effects. However, this performance was slightly lower compared to some studies employing the PLSR model and concentrating on bare soil pixels, such as Selige et al. [70] reporting an $R_{cv}^2$ of 0.71 from HyMap, Gomez et al. [43] reporting an $R_{val}^2$ of 0.75 from AisaDual data, and Nouri et al. [10] reporting an $R_{test}^2$ of 0.72 from HyMap data. Finally, it was in accordance with the performance obtained by Gomez et al. [5], reporting an $R_{cv}^2$ of 0.64 using HyMap hyperspectral data for clay content estimation.

The large range of performances of topsoil clay content prediction using hyperspectral airborne sensors reported in the literature could be attributed to several factors, such as the difference in the number of soil samples used to train the regression models [11], the soil surface roughness and moisture [23], the pedological context and soil diversity in each study area, the residual PV or NPV cover on the selected bare soil pixels [3], and the variability in spatial resolution [43].

5.3. Model Performances Depending on Bare Soil Fractional Cover

The performance in clay content estimations decreased slightly in terms of $R_{val}^2$ and RMSEP as the bare soil fractional cover decreased for both M_p and M_{all_p} models (Table 3 and Table 4). This result is consistent with the literature, as a better selection of bare soil pixels by the stricter threshold would ensure lesser spectral influence from non-soil components [23].

The M_{all_p} models exhibited satisfactory performances even at lower bare soil fractional cover, with $R_{val}^2$ consistently exceeding 0.5. This could be attributed to the utilization of the maximum number of soil samples by these models in their calibration. These results were consistent with those of Dvorakova et al. [71], where models predicting SOC showed only slight variations across a Cellulose Absorption Index (CAI) threshold ranging from 0.00 to 3.00. In that study, the R² values increased from a CAI threshold of 0.00 to reach the maximum at 0.75 (R² = 0.59), then decreased until a threshold of 1.25, and finally remained relatively stable for CAI thresholds between 1.50 and 3.00 (R² around 0.47) [71]. However, this stability in performance contrasts with the observations by Bartholomeus et al. [16], who conducted simulations on the impact of PV on SOC estimation over maize fields and reported significant inaccuracies beyond 5% of PV cover. Similarly, Ouerghemmi et al. [17] reported similar results for clay content estimation over vineyards, noting a significant decrease in accuracy beyond 15% of PV cover.

Finally, following the work of Dutta et al. [11] analyzing the impact of the number of samples in the calibration dataset on model performance or the work of Vaudour et al. [23] analyzing the impact of soil roughness and moisture on soil property predictions, future research could evaluate the impact of these factors in the proposed approach. It is also to be noted that while each model M_{all_p} was built using pixels of classes C_p to C₉ to obtain the best possible performance for clay content mode ling, it was applied only to pixels of the unique class C_p. Increasing the spatial sampling density to have more soil samples per bare soil fractional cover class would allow extensive studies on the proposed composite mapping approach by building and testing regression models on each of the bare soil fractional cover class C_p.

5.4. Clay Content Composite Mapping Approach

The accuracy of the spatialized clay predictions in this study with $R_{val}^2$, ranging from 0.53 to 0.63 depending on the bare soil fractional cover classes (Table 4), fell within the range of accuracy obtained from bare soil composite maps in the literature. For example, Gasmi et al. [72] obtained $R_{val}^2$ of 0.64 for clay prediction using a bare soil composite map based on reflectance median along multispectral Landsat-TM time-series images, while Sorenson et al. [31] obtained $R_{val}^2$ of 0.44 for clay prediction using a bare soil composite map based on Landsat-5 time-series images combined with legacy data.

Despite our modest prediction accuracies, the proposed clay content composite map showed a correct soil pattern, regardless of the bare soil fractional cover classes, where clay-rich areas were identified over black soils found in the valley bottoms around the riverbeds and low clay content areas were identified over red soils found in the uplands and hillslopes. While exercising caution in the interpretation of clay predictions at the local scale, such as the field level, our model demonstrates its utility in providing detailed and nuanced soil pattern recognition at the regional scale. This is particularly evident in its ability to discriminate between black and red soils, offering valuable insights into the heterogeneity of the soil landscape.

The standard deviation of clay predictions along the bootstrap process was considered an indicator of the prediction uncertainty (Figure 3c). The bootstrap, as proposed by Breiman [73], measures the uncertainty of predictions by generating several models based on different subsets that are simulated from a single dataset. Our approach to estimating and mapping the uncertainties of prediction follows the work of Brodský et al. [74], Gomez et al. [43], and more recently, Castaldi et al. [34]. Brodský et al. [74] proposed a stochastic simulation (100 realizations) to assess the uncertainty as the standard deviation of predictions, Gomez et al. [43] proposed a bootstrap procedure (999 iterations) to obtain the pixel-wise variance of predictions. Castaldi et al. [34] proposed to estimate the pixel-wise uncertainty as the ratio between the standard deviation and the mean of predictions. While both standard deviation and variance reflect the variability in a distribution, the standard deviation is expressed in the same unit as the original data (% of clay content, in this case), which can be useful for evaluating both the prediction and uncertainty maps. As these prediction uncertainties were uncorrelated to the bare soil fractional cover (Figure 3a,c), they cannot be attributed to the bare soil fractional cover. The prediction uncertainties could also be related to soil roughness and soil moisture, which may impact the reflectance spectra and, in turn, the predictions, as shown by Ben-Dor et al. [75], Wu et al. [76], and Denis et al. [77]. In the context of this study, the soil surface roughness and soil moisture might be highly variable as they depend on the type and intensity of the irrigation and tillage operations, respectively, and the study area is characterized by small sized crop fields and frequent crop rotation throughout the year [38]. So, as demonstrated by Vaudour et al. [23] and Urbina-Salazar et al. [33], future studies would have to consider the influence of several spectral perturbing factors and environmental parameters that would affect the soil property to be estimated, including soil surface characteristics such as soil moisture and roughness.

5.5. Comparison with Other Approaches to Extend Predictions

To overcome the issue of spectral mixing in imaging soil spectroscopy, several approaches have been proposed to extend the mapping capability over bare soil pixels. A first approach consists of (i) extracting the spectrum of “pure soils” from a mixed spectrum of soil and vegetation and then (ii) using it to predict the targeted soil properties [16,17,18]. A second approach consists of (i) creating a bare soil composite map from multi-temporal data, acquired from a single sensor such as Sentinel-2 [23,33] or Landsat [20,78] or multiple sensors [24,30] and then (ii) using it to predict the targeted soil properties [23,24]. These studies synthesize bare soil composites using the median or mean spectral reflectance of bare soil pixels identified using spectral indices like NDVI, BSI, and NBR2 [20,21,22,33,72]. While most of these studies used multispectral satellite data to create a bare soil composite map from multi-temporal data, Diek et al. [79] produced a bare soil composite map from multi-temporal hyperspectral data acquired from the Airborne Prism EXperiment (APEX) imaging sensor.

In contrast, our study used a single AVIRIS-NG image and thus benefit ed from the relatively finer spatial and spectral resolution inherent in airborne hyperspectral data. This finer spatial resolution of AVIRIS-NG data (4 m) is useful for mapping soil properties over small-sized fields with highly diverse vegetation types, whereas the spatial resolution of Sentinel-2 or Landsat data (10–30 m) might present limitations in capturing soil pixels in such fragmented landscapes. Additionally, the finer spectral resolution of the hyperspectral AVIRIS-NG data (5 nm) is useful in identifying dry vegetation (which is often a major component covering crop fields, especially after harvest), while the limited number of spectral bands and coarse spectral resolution of Sentinel-2 or Landsat may be insufficient [80].

Finally, airborne hyperspectral data have limited availability due to the expenses associated with the flight and the requirement of clear sky conditions for airborne platforms. The availability of data from new generation spaceborne hyperspectral sensors like the Italian PRecursore IperSpettrale della Missione Applicativa (PRISMA) [81], American Hyperspectral InfraRed Imager (HyspIRI) [82], Japanese Hyperspectral Imager Suite (HISUI) [83], and German Environmental Mapping and Analysis Program (EnMAP) [84] could provide multi-temporal hyperspectral observations. Combining our proposed approach with the ones based on bare soil composite images using multi-temporal data [23,71,79] could significantly increase soil property estimation coverage and help provide associated prediction uncertainties.

6. Conclusions

This study proposed a novel clay content composite mapping approach for extending topsoil clay content from an airborne AVIRIS-NG hyperspectral image over a heterogeneous agricultural study area. The approach classified the image pixels according to their bare soil fractional cover determined by spectral unmixing and assigned specific regression models to each bare soil fractional cover class to estimate the clay content along with the uncertainty in prediction. It was observed that the model performance increases with the adoption of higher bare soil fractional cover thresholds in terms of an increase in R² and a decrease in RMSEP values. Despite the modest prediction accuracy, the clay content composite map showed a correct spatial pattern regardless of the bare soil fractional cover classes. The potential mapped area in terms of clay content in this study ranged from 10.4% to 52.7% with a bare soil fraction threshold of >0.7 and >0.3, respectively, and the compositing approach allowed an extension of the mapped surface by 42.42%. The proposed approach could be adopted to extend the mapping capability of planned and current hyperspectral satellite missions that have relatively coarser spatial resolutions, which amplifies the spectral mixing problem. In such cases, a combination of the proposed approach and multi-temporal images could significantly increase both the extent of the mapped area and associated prediction performances. Further, additional research will be needed to mitigate the influence of factors that affect the spectral reflectance of soils, such as moisture and surface roughness, so as to improve the prediction performance.