Spatial Variability in Soil Attributes and Multispectral Indices in a Forage Cactus Field Irrigated with Wastewater in the Brazilian Semiarid Region

Abstract

Multispectral images obtained from Unmanned Aerial Vehicles (UAVs) have become strategic tools in precision agriculture, particularly for analyzing spatial variability in soil attributes. This study aimed to evaluate the spatial distribution of soil electrical (EC) and total organic carbon (TOC) in irrigated forage cactus areas in the Brazilian semiarid region, using field measurements and UAV-based multispectral imagery. The study was conducted in a communal agricultural settlement located in the Mimoso Alluvial Valley (MAV), where EC and TOC were measured at 96 points, and seven biophysical indices were derived from UAV multispectral imagery. Geostatistical models, including cokriging with spectral indices (NDVI, EVI, GDVI, SAVI, and NDSI), were applied to map soil attributes at different spatial scales. Cokriging improved the spatial prediction of EC and TOC by reducing uncertainty and increasing mapping accuracy. The standard deviation of EC decreased from 1.39 (kriging) to 0.67 (cokriging with EVI), and for TOC from 15.55 to 8.78 (cokriging with NDVI and NDSI), reflecting a 43.5% reduction in uncertainty. The indices, EVI, NDVI, and NDSI, showed strong potential in representing and enhancing the spatial variability in soil attributes. NDVI and NDSI were particularly effective at finer grid resolutions, supporting more efficient irrigation strategies and sustainable agricultural practices.

1. Introduction

Remote sensing (RS) is an essential tool in precision agriculture, with RGB and multispectral images widely used to estimate agricultural productivity and adjust management strategies, especially in intercropping systems [1,2]. The evolution of RS techniques has brought increasing relevance to Unmanned Aerial Vehicles (UAVs) equipped with multispectral cameras, which provide greater spatial resolution and data acquisition frequency compared with satellite images [3,4,5,6]. Silva et al. [7] highlighted the potential of biophysical indices derived from RGB images captured by UAVs for monitoring the vegetative state of forage cactus, as well as for detecting bare soils and salinized areas, factors that contribute to more sustainable agricultural practices [8,9].

Biophysical indices derived from UAV images offer robust support for agricultural monitoring and environmental diagnostics. Among the vegetation indices, the Normalized Difference Vegetation Index (NDVI) is one of the most widely used due to its sensitivity to photosynthetically active biomass, exploiting the contrast between red and near-infrared reflectance [10]. The Soil-Adjusted Vegetation Index (SAVI) and the Enhanced Vegetation Index (EVI) are also commonly applied to minimize the influence of atmospheric and soil background conditions, respectively, and these adjustments are particularly relevant in canopies with sparse or heterogeneous vegetation, such as those found in forage cactus systems [11,12,13]. In addition, the Modified Green–Red Vegetation Index (MGRVI) and the Visible Atmospheric Resistance Index (VARI), which use only the visible spectral bands, are effective under variable illumination and in conditions where near-infrared (NIR) data are limited [14].

To assess salinity, the Normalized Difference Salinity Index (NDSI) was adopted for its ability to highlight spectral differences caused by salt crusts or saline moisture on the soil surface, using the red and NIR bands [15]. Studies such as [16] have demonstrated strong correlations between NDSI and soil salinity in irrigated areas. The Generalized Difference Vegetation Index (GDVI) was also considered due to its flexibility in different spectral configurations and its accuracy in estimating vegetation density [17]. Collectively, these indices provide a comprehensive spectral framework that enhances spatial modeling capabilities when integrated with geostatistical techniques.

Recent studies reinforce the value of RGB-based indices in estimating parameters such as leaf area index (LAI) and chlorophyll content via machine learning [18,19,20]. Santos et al. [18] and [12,21] demonstrated the relevance of these indices for productive and environmental monitoring, while [2] showed that multispectral images, when including bands beyond the visible spectrum, outperform RGB images in the classification of forage cactus using artificial intelligence algorithms. According to [12], the integration of biophysical indices from visible and non-visible spectral bands improves the mapping of field conditions. Feng et al. [22] also reported a higher model accuracy with multispectral images compared with RGB-only data.

These indices allow the indirect estimation of soil attributes through regression models with high coefficients of determination [8] and assist in spatialization by interpolation methods, such as kriging. Among these methods, cokriging stands out for incorporating secondary variables (e.g., spectral indices) correlated with the primary variable, thus reducing the uncertainty of the estimate. Studies have shown that cokriging produces smaller prediction errors than conventional techniques, such as ordinary kriging and inverse distance weighting (IDW) [23,24,25].

At the same time, the Brazilian semiarid region presents unique environmental challenges due to its unpredictable precipitation patterns, shallow and compactable soils, and frequent salinization [26]. Marengo et al. [27] describes it as the wettest semiarid region in the world, with annual precipitation of up to 800 mm, but with high evapotranspiration (~2000 mm/year), which restricts groundwater recharge and compromises water security. The limited sanitation infrastructure, present in only 243 of the 1135 municipalities [28], contributes to the risks of contamination and reduced water quality [29]. Initiatives that use treated domestic wastewater for irrigation have emerged as a sustainable alternative. However, their reach is still limited, especially in rural areas [30]. Although positive results have been observed in cactus pear irrigated with treated wastewater [31], there are still concerns about environmental risks, such as microbial contamination and, especially, soil salinization [32]. On the other hand, the controlled reuse of wastewater can increase soil organic matter content [33].

Forage cactus is a strategic plant for semiarid regions, characterized by intense radiation, water scarcity, and high temperatures, affecting the availability of forage food [34]. The productivity of this crop is strongly correlated with the physical–chemical quality of the soil, with organic carbon content (TOC) being one of the main indicators of fertility in semiarid environments [35]. TOC contributes to greater water retention, aggregate stability, and nutrient cycling, both of which are essential for the development of crops under a water deficit. These benefits are enhanced by the use of forage species adapted to semiarid ecosystems, which generally accumulate more carbon when compared with agricultural crops due to the high density of their fasciculate root system and their ability to use small amounts of water available in the soil [36].

Soil electrical conductivity (EC), which is often high in areas irrigated with reclaimed water, acts as a salinization agent, limiting water absorption by agricultural crops, affecting their productivity. In this context, the spatial mapping of these attributes, especially with the support of remote sensing techniques, allows more effective interventions in irrigation management and the use of conservation practices, providing increased productive resilience of forage cactus. Cheng et al. [37] found that remote sensing can be used to map crop water productivity, leading to significant water savings while maintaining or increasing crop yields. Pizarro et al. [38] demonstrated that multispectral indices derived from UAVs can improve the prediction of soil properties, contributing to the more efficient and effective management of agricultural crops.

Despite the advances in RS and geotechnologies, there is still a scientific gap regarding the integration of RGB and multispectral indices derived from UAVs as auxiliary variables in cokriging models aimed at estimating soil attributes in areas irrigated with reused water. This approach remains little explored in forage cactus cultivation, especially in the Brazilian semiarid region. Therefore, this study aimed to evaluate the effectiveness of biophysical indices in reducing the uncertainties associated with the geostatistical mapping of organic matter and electrical conductivity in a forage cactus production area irrigated with reused water, in the semiarid region of Pernambuco, Brazil. In addition, we sought to investigate the sensitivity of cokriging to the spatial resolution of secondary variables, comparing different interpolation scales. This research offers a practical and scientific contribution to precision agriculture by proposing an innovative, relatively low-cost approach based on the integration of spectral data from UAVs to improve soil monitoring and increase the efficiency of wastewater irrigation in regions of high environmental vulnerability.

2. Material and Methods

2.1. Study Area

This study was carried out in an experimental agricultural area located in the Nossa Senhora do Rosário settlement, within the Mimoso representative basin, in the municipality of Pesqueira-PE. The Mimoso basin is part of the Alto Ipanema basin, situated in the semiarid region of the State of Pernambuco (Figure 1), and it serves as a tributary to the São Francisco River in its middle lower section. The climate of the region is classified as BSh (hot semiarid) according to the Köppen classification. The average annual reference potential evapotranspiration is 2.000 mm, and the average annual temperature is 23 °C. The region receives an average annual precipitation of 630 mm [39].

The soil in the study area is classified as a Regolithic Neosol, according to the Brazilian Soil Classification System [40]. The physical and chemical characteristics of the soil, based on saturated paste extracts, are shown in Table 1.

In the experimental plot, forage cactus (Opuntia stricta (Haw.) Haw.) was transplanted on 15 May 2022 and spaced 1.0 m between rows and 0.2 m between plants. Irrigation was applied using a localized drip system supplied with treated domestic wastewater, following the reuse system described by Mayer et al. [30]. The irrigation schedule was set at three-day intervals and based on 70% of the reference evapotranspiration (ET_o), calculated using the Penman–Monteith method [41].

2.2. Soil Salinity and Soil Total Organic Matter

Soil attribute measurements were carried out using a regular grid of 96 sampling points, spaced 2 × 1 m apart, at a depth of 0–0.1 m. Soil electrical conductivity (EC) was determined using the saturated paste extraction method, as described by Teixeira et al. [42]. Soil carbon and organic matter dynamics were analyzed following the method of Yeomans and Bremner [43], based on carbon oxidation with potassium dichromate (K₂Cr₂O₇) and subsequent titration with ammonium ferrous sulfate (0.4 mol L⁻¹), enabling the estimation of total organic carbon (TOC).

2.3. Multispectral Monitoring with Unmanned Aerial Vehicle (UAV)

Multispectral images were acquired using a DJI Phantom 4 Multispectral drone, which offers a spatial resolution of 0.5 m. The drone is equipped with an RGB camera and a set of multispectral cameras, which captures data across five spectral bands, blue (450 nm), green (560 nm), red (650 nm), red edge (730 nm), and near-infrared (840 nm), each with 2 MP resolution and global shutter, mounted on a stabilized 3-axis gimbal [9]. The flight was conducted in October 2023, following technical specifications and Brazilian Civil Aviation Agency (ANAC) regulations. Flight parameters included a speed of 2 m/s and an altitude of 40 m. The weather conditions on the day of capture were characterized by a maximum temperature of 32.37 °C, a minimum temperature of 20.57 °C, and the occurrence of 0.6 mm of precipitation. Cumulative rainfall in the preceding periods was 4.9 mm in the last 30 days, 35.9 mm over 60 days, and 77.0 mm over the previous 90 days.

An automatic survey mission was executed with 75% front and lateral image overlap across all study dates, achieving a Ground Sample Distance (GSD) of ~2 cm² for high-resolution analysis. Reflectance correction was performed using a Calibrated Reflectance Panel (CRP) to minimize lighting variability and enable accurate vegetation index calculations. All images met Metashape quality standards (scores > 0.75) after preliminary assessment. A set of biophysical indices were selected that explore both the visible spectrum bands (RGB) and the near-infrared (NIR) bands, as shown in Table 2.

The index obtained from visible spectrum bands (RGB) will be referred to as RGB indices, while those using bands both within and beyond the visible spectrum (RGB, RE, and NIR) will be referred to as multispectral indices. To calculate and extract index values at field sampling points, the Raster Calculator in QGIS^® software version 3.34.10 [50] was used. Spectral data from multispectral images were extracted based on three different sampling grids (Figure 2) to assess the influence of sampling density on correlations with primary soil variables, using a single image. This higher number of sampling points could reduce the uncertainties in soil variable map estimations. The extraction followed the original grid (grid 1) used in field sampling, with a spacing of 2 × 1 m, totaling 96 points. Regarding grid 2 (1 × 1 m), 179 points were extracted, and for grid 3 (0.5 × 0.5 m) 636 points were considered. The higher-resolution grids were generated through the spatial up-scaling of the UAV-derived data, without additional field soil sampling, to assess the influence of extraction resolution on the correlation analysis and geostatistical predictions.

2.4. Statistical Analysis

The data was analyzed using descriptive statistics, including the mean, minimum and maximum values, first and third quartiles, standard deviation, and coefficient of variation (CV). The CV was classified according to the criteria proposed by Warrick and Nielsen [51]: low variability (CV ≤ 12%), moderate variability (12% < CV ≤ 60%), and high variability (CV > 60%). Data normality was assessed using the Kolmogorov–Smirnov test at a 0.05 significance level, and logarithmic transformations were applied when necessary. Trend analysis was performed using linear models based on the coefficient of determination, following the methodology of Andrade et al. [52]. Outliers were identified and removed according to Hoaglin et al. [53], considering values below the lower limit (LL) or above the upper limit (UL).

Additionally, Spearman’s rank correlation and Principal Component Analysis (PCA) were performed to identify potential linear relationships between multispectral indices and soil attributes. PCA was conducted using a custom script in R software (version 4.4.2, 2022), with the FactoMineR package [54].

Spatial dependence was assessed using classical semivariograms, constructed from semivariance estimates according to Equation (1) [55]. The fitting of theoretical semivariogram models was evaluated using the Geo-EAS software Version 1.2.1 (Geostatistical Environmental Assessment Software, Environmental Monitoring Systems Laboratory, Las Vegas, NV, USA) [56] and GS+^® version 7.0 (Gamma Design Software, Plainwell, MI, USA) [57].

(1)$\widehat{\mathsf{\gamma }}\left(\mathrm{h}\right)={\displaystyle \frac{1}{2\mathrm{N}\left(\mathrm{h}\right)}}{\sum }_{\mathrm{i}=1}^{\mathrm{N}\left(\mathrm{h}\right)}{\left[\mathrm{Z}\left({\mathrm{X}}_{\mathrm{i}}+\mathrm{h}\right)-\mathrm{Z}\left({\mathrm{X}}_{\mathrm{i}}\right)\right]}^2$

where γ(h) = estimated value of the semivariance of the experimental data; Z(X_i + h) and Z(X_i) = observed values of the regionalized variable; N(h) = number of pairs of measured values, separated by a distance h; and X_i and X_i + h refer to two data points (locations where the variable Z is measured) that are a lag distance, h, apart.

Based on the experimental semivariogram, the data were fitted to theoretical models, testing spherical, exponential, and Gaussian models. The mathematical fitting enabled the estimation of key parameters: nugget effect (C₀), range of spatial dependence (a), and sill (C₀ + C₁). The degree of spatial dependence (DSD) was classified according to Cambardella et al. [58], where the degree of dependence is considered strong, moderate, or weak when the ratio of nugget effect (C₀) to sill (C₀ + C₁) is less than 25%, between 25% and 75%, and greater than 75%, respectively.

After fitting the semivariograms, model validation was performed using the leave-one-out cross-validation method [59], where each observed value is temporarily removed and estimated via kriging, then replaced with the predicted value. The distribution of standardized errors is then analyzed, which should present a mean close to zero and a standard deviation near one [60]. Following semivariogram validation, kriging was performed using the fitted models in Surfer^®, version 13.6.618 [61]. In addition to contour maps, standard deviation maps were also generated.

The spatial correlation between field variables and multispectral indices was assessed using a cross-semivariogram methodology. This correlation between distinct variables, known as inter-regionalization, was described by Goovaerts [60] and can be quantified using the cross-semivariogram or cross-covariance, as presented in Equation (2). Standard deviation values were extracted to analyze the frequency distribution and assess the consistency of the correlations obtained with the spectral indices.

(2)$\widehat{{\mathsf{\gamma }}_{12}}\left(\mathrm{h}\right)={\displaystyle \frac{1}{2\mathrm{N}\left(\mathrm{h}\right)}}{\sum }_{\mathrm{i}=1}^{\mathrm{N}\left(\mathrm{h}\right)}\left[{\mathrm{Z}}_1\left({\mathrm{X}}_{\mathrm{i}}\right)-{\mathrm{Z}}_1\left({\mathrm{X}}_{\mathrm{i}}+\mathrm{h}\right)\right][{\mathrm{Z}}_2({\mathrm{X}}_{\mathrm{i}})-{\mathrm{Z}}_2({\mathrm{X}}_{\mathrm{i}}+\mathrm{h})]$

where $\widehat{{\mathsf{\gamma }}_{12}}\left(\mathrm{h}\right)$ is the estimated value of the cross semivariogram of the experimental data; Z(X_i + h) and Z(X_i) are the observed values of the regionalized variable; Z₁ and Z₂ are spatially correlated variables; N(h) is the number of pairs of measured values separated by a distance, h; and Xi and X_i + h refer to two data points (locations where the variable Z is measured) that are a lag distance, h, apart.

3. Results and Discussion

3.1. Descriptive Analysis of Soil Attributes and Spectral Indices

Table 3 presents the descriptive statistics for soil electrical conductivity (EC), organic matter (OM), and the indices derived from UAV imagery, using three different sampling grid resolutions. The mean organic matter content was 33.831 g kg⁻¹, while the mean electrical conductivity was 1.83 dS m⁻¹. Based on the coefficients of variation (CVs), EC was classified as highly variable (85.4%), whereas organic matter showed moderate variability (27.6%). The high variability observed in EC may be related to localized factors such as variations in irrigation management with reclaimed water and topographic differences across the field [62].

Regarding the biophysical indices, the Enhanced Vegetation Index (EVI) and the Generalized Difference Vegetation Index (GDVI) showed the highest mean values (0.426 and 0.398, respectively). In contrast, the Normalized Difference Salinity Index (NDSI) and the Visible Atmospherically Resistant Index (VARI) exhibited negative values, which is a common pattern in areas with low vegetation cover or exposed soils [12,63]. The variability in the indices ranged from low to moderate, with coefficients of variation (CV) between 17.5% (NDVI) and 22% (EVI), indicating good consistency in the data, even when using a lower-resolution sampling grid.

Considering grid 2, a slight rise in the mean values of EVI and GDVI was observed. The data range remained stable, with minimum and maximum values similar to those of the previous grid, and the CVs remained comparable. This suggests that increasing the sampling density did not significantly affect the variability in the extracted indices. The lower (LL) and upper (UL) limits confirmed that most data were concentrated within a statistically acceptable confidence interval, with few outliers. In the highest-resolution grid (grid 3), a further increase in the mean values of EVI (0.440) and GDVI (0.419) was recorded, suggesting improved spatial detail in detecting variability in soil spectral parameters. However, the standard deviation values remained similar to those of the previous grids—for example, EVI showed standard deviations of 0.094 (grid 1), 0.084 (grid 2), and 0.092 (grid 3); for GDVI, the values were 0.072 (grid 1), 0.073 (grid 2), and 0.079 (grid 3)—indicating that increasing the sampling density did not significantly reduce data dispersion or uncertainty.

3.2. Correlations Between Soil Attributes and Spectral Indices

Table 4 presents the Spearman correlations between soil attributes and spectral indices across the three sampling resolutions (grids 1, 2, and 3), allowing the consistency of biophysical and physicochemical variables to be assessed under varying sampling densities. Overall, strong positive correlations were observed among the spectral indices EVI, NDVI, GDVI, and SAVI, particularly in grids 2 and 3. GDVI and SAVI exhibited near-perfect correlations across all grids (≈0.999), suggesting that both capture similar spectral patterns associated with green biomass. NDVI and GDVI also showed high correlations (>0.95), reinforcing their robustness in vegetation detection in irrigated areas. Conversely, NDSI showed a strong negative correlation with vegetation indices such as GDVI and NDVI (<–0.95), indicating that increased salinity tends to be associated with reduced vegetation cover [64].

Soil electrical conductivity (EC) showed moderate correlations with the spectral indices of NDVI, GDVI, EVI, and SAVI, with values of 0.209, 0.172, 0.17, and 0.162, respectively. The NDSI exhibited a negative correlation (−0.16), suggesting an indirect relationship between soil salinity and the presence of forage cactus in the study area. In contrast, organic matter (OM) showed weaker correlations. Overall, OM displayed weak or nearly null correlations with the spectral indices, corroborating the findings of Belenok et al. [65], who reported a low linear correlation between the OM and NDVI values derived from Landsat 5 and 8 imageries, based on Pearson’s test. The correlation patterns were consistent across all three sampling grids, indicating that changes in sampling density did not significantly affect the strength of the relationships among variables.

Principal Component Analysis (PCA) was applied to the dataset, excluding outliers, to identify correlation patterns between soil attributes and spectral indices. The resulting biplots for the 2 × 1 m, 1 × 1 m, and 0.5 × 0.5 m grids are presented in Figure 3. Dimension 1 (Dim1) explained between 40.6% and 42.2% of the variance in all three grids, while Dimension 2 (Dim2) accounted for 27.5% to 30.6%, together representing over 70% of the total variability explained by the first two components.

The biplots reveal that EC and OM contributed minimally to the variance explained by the first two principal components. However, in the 2 × 1 m grid, soil electrical conductivity showed a stronger correlation with the spectral indices NDVI, SAVI, GDVI, and EVI, suggesting that these indices have the potential for estimating soil salinity patterns in areas irrigated with reclaimed water. A similar behavior was observed in the higher resolution grids (1 × 1 m and 0.5 × 0.5 m), where EC maintained a certain degree of proximity to the spectral vectors, albeit with reduced intensity. The consistency of this pattern across different grids reinforces the robustness of PCA in capturing the structural relationships between spectral data and soil chemical attributes. Nevertheless, this result contrasts with the findings of Silva et al. [8], who reported inverse correlations between NDVI, SAVI, GDVI, and EVI and soil EC in irrigated areas of the semiarid region. This divergence can be attributed to the fact that the present study focused solely on the spectral response of the soil, whereas Silva et al. [8] evaluated the spectral response of the plant. Variations in soil type, the degree of salinization, irrigation management, soil cover, and image acquisition period also directly affect spectral responses and their relationship with soil properties.

3.3. Geostatistical Modeling of Soil Attributes and Spectral Indices

Table 5 and Figure 4 present the parameters of the fitted semivariogram models for the physicochemical soil variables (OM and EC) and the spectral indices (MGRVI, NDSI, NDVI, VARI, EVI, GDVI, and SAVI). All variables exhibited spatial dependence, with fitted exponential or Gaussian models and high coefficients of determination (R²), indicating good model performance.

Organic matter (OM) exhibited a strong spatial dependence (DSD = 11.15%), while electrical conductivity (EC) showed a moderate dependence (47.98%). The range values varied from 1.08 to 3.36 m, with EC presenting the highest range, likely due to its greater spatial heterogeneity, often influenced by irrigation practices and salt accumulation on the soil surface [33]. The spatial dependence observed for organic matter is similar to that reported by Fu et al. [66], who evaluated soil organic matter at depths of 0–20 cm and 20–40 cm and found models with DSD values of 10% and 13%, respectively.

Among the spectral indices, EVI presented the highest coefficient of determination (R² = 0.993), followed by NDSI and NDVI (0.987), GDVI (0.937), and SAVI (0.930), indicating an excellent fit of the theoretical model to the experimental semivariances. In contrast, VARI showed the lowest R² (0.671), likely due to its higher sensitivity to atmospheric interference and illumination variability, which is common in indices based solely on visible bands [23]. In terms of spatial dependence, VARI exhibited a strong spatial dependence (DSD < 25%), while the other indices—MGRVI, NDSI, NDVI, EVI, GDVI, and SAVI—demonstrated a moderate spatial dependence (DSD between 38.38% and 48.03%). The observed range values varied between 3.95 m (VARI) and 56.25 m (MGRVI), reflecting different spatial extents.

Table 6 and Figure 4 present the parameters of the semivariograms adjusted for the spectral indices in both grids 2 (1 × 1 m) and 3 (0.5 × 0.5). The adjustments of the semivariograms of the spectral indices showed different patterns of spatial structure, with a predominance of the Gaussian and exponential models.

The spectral indices MGRVI, NDSI, NDVI, EVI, and GDVI were best represented by the Gaussian model, with high coefficients of determination (R²) ranging from 0.986 (GDVI) to 0.850 (NDVI and NDSI), indicating an excellent fit of the data to the modeled spatial structure. On the other hand, the SAVI and VARI indices showed a better fit to the exponential model. The VARI stood out for having the highest range (R = 23.04 m) along with a good R² (0.886), while EVI (R = 16.47 m) and MGRVI (R = 12.30 m) also exhibited ranges above 10 m, indicating significant spatial continuity. Regarding the degree of spatial dependence (DSD), SAVI showed a strong dependence (13.67%), while the other indices exhibited a moderate dependence, with MGRVI (49.93%), NDVI (49.76%), and NDSI (49.76%) standing out.

In grid 3 (0.5 × 0.5 m), the model fitting patterns remained consistent, particularly for indices based on near-infrared (NIR) bands, such as NDVI, GDVI, SAVI, and EVI, reinforcing the robustness of these indices in representing the soil’s spectral response. GDVI showed the highest R² (0.928), followed by SAVI (0.927), and both NDVI and NDSI (0.910), indicating their strong capacity to represent the spatial structure of the data even under a higher sampling density. Although VARI (0.515) and MGRVI (0.556) presented lower R² values, all indices passed the leave-one-out cross-validation test, reinforcing the accuracy of the spatial models. Furthermore, the DSD values in grid 3 followed patterns similar to the previous grids: VARI (11.59%) and MGRVI (13.90%) maintained a strong spatial dependence, while EVI, GDVI, SAVI, NDVI, and NDSI showed a moderate dependence, ranging from 49.78% to 57.52%. The spatial range varied from 2.64 m (VARI) to 10.62 m (EVI), indicating that the spatial dependence structure was preserved even at higher spatial resolution.

The parameters of the cross-semivariogram models between organic matter (OM) and the spectral indices are presented in Table 7 and Figure 5. The Gaussian models were predominant, with EVI showing the highest spatial dependence (DSD = 26.36%), followed by SAVI with a very low DSD (0.17%) and a range of 16.09 m, indicating a well-defined spatial structure. NDVI and NDSI presented similar low DSD values (0.18%), suggesting a strong spatial correlation with OM. GDVI had a range of 12.66 m and DSD of 8.33%, indicating a moderate spatial dependence. MGRVI and VARI did not yield valid semivariogram models at this scale, indicating limited applicability in coarser sampling grids.

The cross-semivariogram adjustments between electrical conductivity (EC) and the spectral indices revealed different levels of spatial dependence, as shown in Table 8 and Figure 6. For this soil attribute, the SAVI exhibited the highest degree of spatial dependence (DSD = 41.91%) with a range of 11.58 m, followed by NDVI and NDSI (DSD = 34.41%), which both fitted with exponential models and showing a broader spatial range (21.09 m). EVI also showed a significant spatial relationship with EC (DSD = 29.85%) and a range of 12.97 m, indicating its potential for mapping soil salinity variations over larger scales. GDVI exhibited the lowest range (7.39 m) and a DSD of 23.24%, suggesting a less continuous spatial structure. MGRVI and VARI did not yield valid model fits at this scale.

3.4. Uncertainty Analysis and Mapping Performance

The frequency distributions of the standard deviations associated with ordinary kriging and cokriging between organic matter and biophysical indices are presented in Figure 7. In grid 1, with the lowest resolution (2 × 1 m), the ordinary kriging of organic matter (OM) yielded a mean standard deviation of 15.55. Cokriging with EVI reduced this value to 9.317, while GDVI, NDVI, and NDSI each reduced it to 13.19. In grid 2 (1 × 1 m), cokriging with EVI, NDSI, and NDVI produced the lowest deviations (9.409, 8.981, and 8.981, respectively), indicating improvements in prediction accuracy. In grid 3 (0.5 × 0.5 m), the lowest deviations were observed with NDSI and NDVI (both at 8.783), followed by EVI (9.317) and MGRVI (15.073).

The estimation of soil organic matter (OM) based on the EVI, GDVI, NDSI, NDVI, and SAVI indices resulted in a reduction in the mean standard deviation in all cases evaluated, except for the GDVI of the 0.5 × 0.5 m grid, where an increase in uncertainty was observed. This improvement in estimates is supported by Gou et al. [67], where the authors highlighted the effectiveness of indices such as NDVI in predicting organic carbon in small properties using UAV-derived images. Overall, it was observed that some indices, such as NDVI and NDSI, performed better as the extraction grid was refined, with successive reductions in mean standard deviations: 13.19, 8.918, and 8.783 for the 2 × 1, 1 × 1, and 0.5 × 0.5 grids, respectively. These results are consistent with those found by Lundgren et al. [68] who evaluated the sensitivity of grid resolution and the inverse relationship with standard deviations. On the other hand, indices like GDVI showed an increase in standard deviation as the grid was refined, suggesting a greater sensitivity to the interpolation scale. Still, all indices evaluated contributed to reducing uncertainty in OM estimates, except for GDVI at the densest grid, indicating the importance of selecting appropriate indices based on the spatial scale of the analysis.

For electrical conductivity, the frequency distributions of the standard deviations associated with ordinary kriging and cokriging are presented in Figure 8. Kriging showed a mean standard deviation of 1.39. In the first grid, cokriging with EVI provided the best result, reducing the standard deviation to 0.666, followed by GDVI (0.686). In grid 2, GDVI showed the lowest deviation (0.396), followed by NDSI and NDVI (both 0.674). In grid 3, EVI again exhibited the lowest deviation among the indices, with a standard deviation of 0.67, followed by SAVI (0.773) and NDVI and NDSI (0.785 and 0.785, respectively).

The frequency distribution of the standard deviations obtained from the mapping shows improvements in the estimates when biophysical indices are used as secondary variables. The best results for estimating EC and OM are shown in Figure 9A,B. OM values were achieved with the 0.5 × 0.5 grid, using the NDVI and NDSI indices to reduce uncertainty in organic matter estimation, with an approximate reduction of 43.52% in the mean standard deviation compared with kriging performed using only the primary variable. However, the frequency distribution of the standard deviations observed in the estimates indicates a significant improvement for the evaluated indices.

Regarding electrical conductivity, the EVI and GDVI indices stood out for having the lowest mean standard deviations, indicating a higher precision in the estimates. These results corroborate previous studies, such as those by Taghadosi and Hasanlou [64], who demonstrated the effectiveness of both indices as indirect indicators of salinity in vegetated areas. Similarly, Silva et al. [8] highlighted GDVI as a promising index for predicting soil electrical conductivity. These findings reinforce the potential of using spectral indices in salinity modeling, especially when applied to environments with vegetative cover.

3.5. Spatial Distribution Maps of Soil Attributes

To complement the uncertainty analysis, spatial distribution maps were generated to visualize the improvement in soil attribute predictions achieved by incorporating spectral indices as secondary variables in the cokriging models (Figure 10).

The spatial distribution maps of OM and EC illustrate the performance of kriging and cokriging approaches. While ordinary kriging (Figure 10A,B) produced more generalized spatial patterns, the integration of spectral indices in the cokriging models (Figure 10C–H) enhanced the spatial detail, particularly for EC, where zones of higher salinity were more sharply delineated. These improvements reflect the ability of multispectral indices to capture the spatial variability associated with vegetation vigor and salinity effects, contributing to more accurate and spatially consistent predictions of soil attributes in wastewater-irrigated systems.

4. Conclusions

This study highlighted the potential of using UAV-derived multispectral images integrated with geostatistical techniques for mapping the spatial variability in soil salinity and total organic carbon in areas irrigated with wastewater in the Brazilian semiarid region. The use of cokriging models incorporating spectral indices as secondary variables enhanced the accuracy of soil attribute predictions and reduced associated uncertainties.

Regarding the indices analyzed, NDVI, EVI, and NDSI stood out as the most promising, particularly NDVI and NDSI, which contributed significantly to reduce uncertainties in spatial interpolation, especially in grids with a higher resolution. These findings emphasize the importance of integrating remote sensing and geostatistics as a strategic tool for monitoring soil quality and supporting more sustainable agricultural practices, thereby enhancing water use efficiency in regions vulnerable to water scarcity. Moreover, the sensitivity analysis across different sampling grids provided valuable insights into the balance between sampling density and model performance, which can guide future operational decisions.

Additionally, the incorporation of NDVI and NDSI into the cokriging models led to reductions of up to 43% in prediction uncertainty for total organic carbon under higher sampling densities. The consistency of the correlations across different grid resolutions further supports the robustness of the proposed approach for practical field applications. By improving soil attribute mapping and supporting more efficient wastewater irrigation management, the proposed methodology contributes to more sustainable agricultural practices and enhances the resilience of forage cactus production systems in vulnerable semiarid environments.

Author Contributions

Conceptualization, E.G.F.A.d.S., T.A.B.A., and A.A.d.A.M.; methodology, E.G.F.A.d.S., T.A.B.A., A.A.d.A.M., L.d.B.d.S., and M.V.d.S.; software, E.G.F.A.d.S., T.A.B.A., M.C.G.d.L., R.E.d.M., and J.A.d.S.d.S.; validation, E.G.F.A.d.S., T.A.B.A., M.C.G.d.L., and R.E.d.M.; formal analysis, E.G.F.A.d.S., T.A.B.A., A.A.d.A.M., M.C.G.d.L., and R.E.d.M.; investigation, E.G.F.A.d.S., T.A.B.A., and A.A.d.A.M.; resources, A.A.d.A.M.; data curation, E.G.F.A.d.S., T.A.B.A., and A.A.d.A.M.; validation, E.G.F.A.d.S., T.A.B.A., M.C.G.d.L., and R.E.d.M.; writing—original draft preparation, E.G.F.A.d.S., T.A.B.A., A.A.d.A.M., L.d.B.d.S., M.V.d.S., M.C.G.d.L., R.E.d.M., and J.A.d.S.d.S.; writing—review and editing, E.G.F.A.d.S., T.A.B.A., L.d.B.d.S., M.V.d.S., M.C.G.d.L., and R.E.d.M.; visualization, A.A.d.A.M., and T.A.B.A.; supervision, A.A.d.A.M., T.A.B.A., and M.V.d.S.; project administration, A.A.d.A.M.; funding acquisition, A.A.d.A.M. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The data presented in this study are available upon request from the corresponding author.

Acknowledgments

We thank the Postgraduate Program in Agricultural Engineering (PGEA) of the Federal Rural University of Pernambuco (UFRPE) for supporting the development of this research.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1 Location of the study area in Pesqueira, Pernambuco state, Brazil, and the experimental site where UAV-based multispectral images were acquired.

Figure 2 Sample grid 1 (2 × 1 m) (A), grid 2 (1 × 1 m) (B) and grid 3 (0.5 × 0.5 m) (C).

Figure 3 Principal Component Analysis (PCA) between soil attributes and spectral indices in grid 1 (2 × 1 m) (A), grid 2 (1 × 1 m) (B), and grid 3 (0.5 × 0.5 m) (C).

Figure 4 Theoretical and experimental semivariograms for EC (A), OM (B), EVI (C), GDVI (D), MGRVI (E), NDSI (F), NDVI (G), SAVI (H), VARI (I).

Figure 5 Cross-semivariograms between Organic Matter (OM) and EVI (A), GDVI (B), NDSI (C), NDVI (D), SAVI (E).

Figure 6 Cross-semivariograms between electrical conductivity (EC) in layer A and EVI (A), GDVI (B), NDSI (C), NDVI (D), SAVI (E).

Figure 7 Frequency distribution for standard deviation of kriging and cokriging of organic matter (OM) with biophysical indices in grid 1 (A), grid 2 (B), and grid 3 (C).

Figure 8 Frequency distribution of kriging and cokriging standard deviation of electrical conductivity (EC) with biophysical indices in grid 1 (A), grid 2 (B) and grid 3 (C).

Figure 9 Frequency distribution for standard deviation of cokriging of organic matter (A) and electrical conductivity (B) with biophysical indices for the best results of each grid.

Figure 10 Spatial distribution maps of organic matter (OM) and soil electrical conductivity (EC) obtained by ordinary kriging (A,B) and cokriging models incorporating spectral indices as secondary variables (C–H).

Table 1

Physical and chemical characteristics of the soil (analyzed using a composite soil sample) in the experimental area, Pesqueira (Pernambuco state).

Depth	Sand	Silt	Clay	Bd	Pd	pH	EC	OC
m	%			g cm−3			dS m−1	g·kg−1
0–0.10	74.8	8.2	17	1.34	2.56	6.3	1.2	13.27
0.10–0.20	68.4	13.3	18.3	1.42	2.42	5.7	1.4	3.93

where Bd = bulk density; Pd = particle density; pH = hydrogen potential (H₂O); EC = electrical conductivity; OC = organic carbon.

Table 2

Spectral indices used for vegetation and salinity assessment based on multispectral images acquired by UAV.

Index	Name	Equation	Reference
NDVI	Normalized Difference Vegetation Index	${\displaystyle \frac{\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{R}}{\mathrm{N}\mathrm{I}\mathrm{R}+\mathrm{R}}}$	Rouse et al. [44]
NDSI	Normalized Difference Salinity Index	${\displaystyle \frac{\mathrm{R}-\mathrm{N}\mathrm{I}\mathrm{R}}{\mathrm{R}+\mathrm{N}\mathrm{I}\mathrm{R}}}$	Khan et al. [45]
MGRVI	Modified Green–Red Vegetation Index	${\displaystyle \frac{{\mathrm{G}}^2-{\mathrm{R}}^2}{{\mathrm{G}}^2+{\mathrm{R}}^2}}$	Bendig et al. [46]
VARI	Visible Atmospherically Resistant Index	${\displaystyle \frac{\mathrm{G}-\mathrm{R}}{\mathrm{G}+\mathrm{R}-\mathrm{B}}}$	Gitelson et al. [47]
SAVI	Soil-Adjusted Vegetation Index	${\displaystyle \frac{\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{R}}{\mathrm{N}\mathrm{I}\mathrm{R}+\mathrm{R}+\mathrm{L}}}(1+\mathrm{L})$	Huete [48]
EVI	Enhanced Vegetation Index	${\displaystyle \frac{2.5\ast (\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{R})}{\mathrm{N}\mathrm{I}\mathrm{R}+6.0\ast \mathrm{R}-7.5\ast \mathrm{B}}}$	Huete et al. [49]
GDVI	Generalized Difference Vegetation Index	${\displaystyle \frac{{\mathrm{N}\mathrm{I}\mathrm{R}}^2-{\mathrm{R}}^2}{{\mathrm{N}\mathrm{I}\mathrm{R}}^2+{\mathrm{R}}^2}}$	Wu et al. [17]

where G = green; B = blue; R = red; NIR = near-infrared; L = soil brightness correction factor (L = 0.5).

Table 3

Descriptive statistics of soil attributes and biophysical indices for different sampling grids.

	Mean	Minimum	Maximum	SD	CV	1° Q	3° Q	LL	UL
A. Grid 1 (2 × 1 m)
EC	1.83	0.118	5.68	1.563	0.854	0.448	3.245	−3.747	7.44
OM	33.831	11.563	61.431	9.337	0.276	26.53	40.108	6.163	60.475
EVI	0.426	0.211	0.712	0.094	0.22	0.359	0.48	0.179	0.66
GDVI	0.398	0.222	0.568	0.072	0.181	0.349	0.441	0.213	0.578
MGRVI	−0.224	−0.313	−0.123	0.04	−0.18	−0.25	−0.201	−0.324	−0.127
NDSI	−0.207	−0.281	−0.119	0.036	−0.175	−0.231	−0.186	−0.299	−0.118
NDVI	0.207	0.119	0.281	0.036	0.175	0.186	0.231	0.118	0.299
SAVI	0.313	0.169	0.468	0.062	0.198	0.271	0.348	0.154	0.465
VARI	−0.17	−0.239	−0.091	0.03	−0.176	−0.189	−0.15	−0.248	−0.091
B. Grid 2 (1 × 1 m)
EVI	0.433	0.211	0.663	0.084	0.194	0.38	0.48	0.231	0.628
GDVI	0.404	0.222	0.595	0.073	0.18	0.358	0.447	0.223	0.582
MGRVI	−0.225	−0.341	−0.103	0.043	−0.193	−0.253	−0.199	−0.334	−0.118
NDSI	−0.21	−0.309	−0.113	0.038	−0.18	−0.236	−0.185	−0.311	−0.11
NDVI	0.21	0.113	0.309	0.038	0.18	0.185	0.236	0.11	0.311
SAVI	0.316	0.169	0.478	0.061	0.192	0.277	0.353	0.162	0.468
VARI	−0.171	−0.254	−0.09	0.031	−0.183	−0.191	−0.149	−0.254	−0.087
C. Grid 3 (0.5 × 0.5 m)
EVI	0.44	0.2	0.721	0.092	0.208	0.38	0.496	0.205	0.671
GDVI	0.419	0.212	0.632	0.079	0.189	0.365	0.468	0.21	0.623
MGRVI	−0.235	−0.382	−0.084	0.054	−0.232	−0.273	−0.2	−0.382	−0.091
NDSI	−0.218	−0.332	−0.107	0.044	−0.203	−0.247	−0.187	−0.336	−0.098
NDVI	0.218	0.107	0.332	0.044	0.203	0.187	0.247	0.098	0.336
SAVI	0.329	0.161	0.516	0.067	0.202	0.283	0.369	0.155	0.498
VARI	−0.177	−0.28	−0.068	0.039	−0.221	−0.203	−0.151	−0.281	−0.074

where SD = Standard Deviation; CV = Coefficient of Variation; Li = Lower Limit; LS = Upper Limit; OM = Organic Matter; EC = Electrical Conductivity; EVI = Enhanced Vegetation Index; GDVI = Green Difference Vegetation Index; MGRVI = Modified Green–Red Vegetation Index; NDSI = Normalized Difference Salinity Index; NDVI = Normalized Difference Vegetation Index; SAVI = Soil-Adjusted Vegetation Index; and VARI = Visible Atmospherically Resistant Index.

Table 4

Spearman’s correlation between soil attributes and spectral indices in grids 1 (2 m × 1 m), 2 (1 m × 1 m), and 3 (0.5 m × 0.5 m), in the area irrigated with reused water in the Mimoso Alluvial Valley.

A. Grid 1
	OM	EC	EVI	GDVI	MGRVI	NDSI	NDVI	SAVI	VARI
OM	1
EC	−0.08	1
EVI	0.015	0.17	1
GDVI	−0.08	0.172	0.724	1
MGRVI	−0.02	0.04	0.656	0.31	1
NDSI	0.084	−0.16	−0.74	−1	−0.32	1
NDVI	−0.11	0.209	0.714	0.957	0.281	−0.96	1
SAVI	−0.08	0.162	0.737	0.999	0.32	−1	0.956	1
VARI	−0.04	0.016	0.608	0.318	0.989	−0.33	0.283	0.325	1
B. Grid 2
	OM	CEA	EVI	GDVI	MGRVI	NDSI	NDVI	SAVI	VARI
OM	1
EC	−0.08	1
EVI	0.015	0.17	1
GDVI	−0.08	0.172	0.596	1
MGRVI	−0.02	0.04	0.644	0.234	1
NDSI	0.084	−0.16	−0.6	−1	−0.24	1
NDVI	−0.11	0.209	0.593	0.953	0.214	−0.95	1
SAVI	−0.08	0.162	0.604	0.999	0.24	−1	0.953	1
VARI	−0.04	0.016	0.618	0.253	0.991	−0.26	0.227	0.257	1
C. Grid 3
	OM	EC	EVI	GDVI	MGRVI	NDSI	NDVI	SAVI	VARI
OM	1
EC	−0.08	1
EVI	0.015	0.17	1
GDVI	−0.08	0.172	0.306	1
MGRVI	−0.02	0.04	0.384	0.204	1
NDSI	0.084	−0.16	−0.31	−1	−0.21	1
NDVI	−0.11	0.209	0.295	0.956	0.18	−0.95	1
SAVI	−0.08	0.162	0.307	0.998	0.212	−1	0.955	1
VARI	−0.04	0.016	0.361	0.228	0.988	−0.24	0.204	0.236	1

where OM = Organic Matter; EC = Soil Electrical Conductivity; EVI = Enhanced Vegetation Index; GDVI = Green Difference Vegetation Index; MGRVI = Modified Green–Red Vegetation Index; NDSI = Normalized Difference Salinity Index; NDVI = Normalized Difference Vegetation Index; SAVI = Soil-Adjusted Vegetation Index; VARI = Visible Atmospheric Resistance Index.

Table 5

Parameters of the semivariograms adjusted for the soil variables and biophysical indices evaluated in the experimental area in grid 1 (2 × 1 m).

Variable	Model	C0	C0 + C1	R	DSD	R2
OM	Exponential	45.0	403.6	3.99	11.15%	0.911
EC	Exponential	1.238	2.58	6.81	47.98%	0.771
MGRVI	Exponential	0.001331	0.002872	56.25	46,34%	0.882
NDSI	Gaussian	0.001099	0.002288	8.3658	48.03%	0.987
NDVI	Gaussian	0.001099	0.002288	8.3658	48.03%	0.987
VARI	Exponential	0.000114	0.000918	3.95	12.41%	0.671
EVI	Gaussian	0.02583	0.05706	11.9165	45.27%	0.993
GDVI	Gaussian	0.0035	0.00793	8.2965	44.13%	0.937
SAVI	Gaussian	0.00221	0.00612	10.35	36.11%	0.930

where OM = Organic Matter; EC = Soil Electrical Conductivity; EVI = Enhanced Vegetation Index; GDVI = Green Difference Vegetation Index; MGRVI = Modified Green–Red Vegetation Index; NDSI = Normalized Difference Salinity Index; NDVI = Normalized Difference Vegetation Index; SAVI = Soil-Adjusted Vegetation Index; VARI = Visible Atmospheric Resistance Index; C₀ = Nugget Effect; C₀ + C₁ = Sill; R = Range; DSD = Dependence Spatial; R² = Coefficient of Determination.

Table 6

Parameters of the semivariograms adjusted for the spectral indices in grid 2 (1 × 1 m) (A) and 3 (0.5 × 0.5) (B).

	Model	C0	C0 + C1	R	DSD	R2
A. Grid 2 (1 × 1)
MGRVI	Gaussian	0.00363	0.00727	12.2976	49.93%	0.855
NDSI	Gaussian	0.001059	0.002128	7.3093	49.76%	0.85
NDVI	Gaussian	0.001059	0.002128	7.3093	49.76%	0.85
VARI	Exponential	0.001634	0.004288	23.04	38.10%	0.886
EVI	Gaussian	0.0663	0.1416	16.4718	46.82%	0.984
GDVI	Gaussian	0.00425	0.0128	8.581	33.20%	0.986
SAVI	Exponential	0.00254	0.01858	15.15	13.67%	0.954
B. Grid 3 (0.5 × 0.5)
MGRVI	Exponential	0.00057	0.0041	3.33	13.90%	0.556
NDSI	Gaussian	0.001499	0.002606	8.5766	57.52%	0.91
NDVI	Gaussian	0.001499	0.002606	8.5766	57.52%	0.91
VARI	Exponential	0.000258	0.002226	2.64	11.59%	0.515
EVI	Gaussian	0.01173	0.02356	10.6175	49.78%	0.756
GDVI	Gaussian	0.00399	0.00799	7.9848	49.93%	0.928
SAVI	Gaussian	0.003143	0.006296	8.0194	49.92%	0.927

where MGRVI = Modified Green–Red Vegetation Index; NDSI = Normalized Difference Salinity Index; NDVI = Normalized Difference Vegetation Index; VARI = Visible Atmospheric Resistance Index; EVI = Enhanced Vegetation Index; GDVI = Green Difference Vegetation Index; SAVI = Soil-Adjusted Vegetation Index; C₀ = Nugget Effect; C₀ + C₁ = Sill; R = Range; DSD = Dependence Spatial; R² = Coefficient of Determination.

Table 7

Fitting parameters of the cross-semivariograms between soil organic matter and spectral indices.

X-Semivar	Model	C0	C0 + C1	R	DSD	R2
OM × EVI	Gaussian	0.075763	0.287383	14.1488	26.36%	0.677
OM × NDSI	Gaussian	−0.0001	−0.0558	14.94	0.18%	0.711
OM × NDVI	Gaussian	0.0001	0.0558	14.94	0.18%	0.711
OM × SAVI	Gaussian	0.0001	0.0574	16.0908	0.17%	0.619
OM × GDVI	Gaussian	0.0093	0.1116	12.6613	8.33%	0.783
OM × MGRVI	EPP	-	-	-	-	-
OM × VARI	EPP	-	-	-	-	-

where MGRVI = Modified Green–Red Vegetation Index; NDSI = Normalized Difference Salinity Index; NDVI = Normalized Difference Vegetation Index; VARI = Visible Atmospheric Resistance Index; EVI = Enhanced Vegetation Index; GDVI = Green Difference Vegetation Index; SAVI = Soil-Adjusted Vegetation Index; C₀ = Nugget Effect; C₀ + C₁ = Sill; R = Range; DSD = Dependence Spatial; R² = Coefficient of Determination.

Table 8

Fitting parameters of the cross-semivariograms between soil electrical conductivity and spectral indices.

X-Semivar	Model	C0	C0 + C1	R	DSD	R2
EC × EVI	Gaussian	0.024	0.0804	12.9731	29.85%	0.775
EC × NDSI	Exponential	−0.01082	−0.03144	21.09	34.41%	0.583
EC × NDVI	Exponential	0.01082	0.03144	21.09	34.41%	0.583
EC × SAVI	Gaussian	0.013502	0.032223	11.5785	41.91%	0.468
EC × GDVI	Gaussian	0.00903	0.03886	7.3959	23.24%	0.655
EC × MGRVI	EPP	-	-	-	-	-
EC × VARI	EPP	-	-	-	-	-

where MGRVI = Modified Green–Red Vegetation Index; NDSI = Normalized Difference Salinity Index; NDVI = Normalized Difference Vegetation Index; VARI = Visible Atmospheric Resistance Index; EVI = Enhanced Vegetation Index; GDVI = Green Difference Vegetation Index; SAVI = Soil-Adjusted Vegetation Index; C₀ = Nugget Effect; C₀ + C₁ = Sill; R = Range; DSD = Dependence Spatial; R² = Coefficient of Determination.

References

1. Lima, G.S.A.; Ferreira, M.E.; Passos, J.S.; Costa, J.V.S.; Sales, J.C.; Santos, M.L. Estimativa de produtividade de arroz a partir de imagens multiespectrais registradas com aeronaves remotamente pilotadas. XX Simpósio Brasileiro de Sensoriamento Remoto; XX SBSR: Florianópolis, SC, Brasil, 2023.

2. Andrade, O.B.d.; Montenegro, A.A.d.A.; Silva Neto, M.A.d.; Sousa, L.d.B.d.; Almeida, T.A.B.; de Lima, J.L.M.P.; Carvalho, A.A.d.; Silva, M.V.d.; Medeiros, V.W.C.d.; Soares, R.G.F. . UAV-Based Classification of Intercropped Forage Cactus: A Comparison of RGB and Multispectral Sample Spaces Using Machine Learning in an Irrigated Area. AgriEngineering; 2024; 6, pp. 509-525. [DOI: https://dx.doi.org/10.3390/agriengineering6010031]

3. Khaliq, A.; Comba, L.; Biglia, A.; Ricauda Aimonino, D.; Chiaberge, M.; Gay, P. Comparison of Satellite and UAV-Based Multispectral Imagery for Vineyard Variability Assessment. Remote Sens.; 2019; 11, 436. [DOI: https://dx.doi.org/10.3390/rs11040436]

4. Mazzia, V.; Comba, L.; Khaliq, A.; Chiaberge, M.; Gay, P. UAV and Machine Learning Based Refinement of a Satellite-Driven Vegetation Index for Precision Agriculture. Sensors; 2020; 20, 2530. [DOI: https://dx.doi.org/10.3390/s20092530]

5. Messina, G.; Peña, J.M.; Vizzari, M.; Modica, G. A Comparison of UAV and Satellites Multispectral Imagery in Monitoring Onion Crop. An Application in the ‘Cipolla Rossa di Tropea’ (Italy). Remote Sens.; 2020; 12, 3424. [DOI: https://dx.doi.org/10.3390/rs12203424]

6. Zou, X.; Li, Z.; Li, Y.; Wang, Y.; Tang, W.; Deng, C.; Cui, J.; Fu, R. A Novel Method to Mitigate the Multipath Error for BDS-2 Dam Deformation Monitoring. Remote Sens.; 2021; 13, 1787. [DOI: https://dx.doi.org/10.3390/rs13091787]

7. Silva, G.I.N.d.; Jardim, A.M.d.R.F.; Santos, W.M.d.; Bezerra, A.C.; Alba, E.; Silva, M.V.d.; Silva, J.L.B.d.; Souza, L.S.B.d.; Marinho, G.T.B.; Montenegro, A.A.d.A. . Estimation of Biophysical Parameters of Forage Cactus Under Different Agricultural Systems Through Vegetation Indices and Machine Learning Using RGB Images Acquired with Unmanned Aerial Vehicles. Agriculture; 2024; 14, 2166. [DOI: https://dx.doi.org/10.3390/agriculture14122166]

8. Silva, D.C.d.; Lopes, P.M.O.; Silva, M.V.d.; Moura, G.B.A.; Nascimento, C.R.; Brito, J.I.B.; França e Silva, Ê.F.d.; Rolim, M.M.; Lima, R.P.d. Principal component analysis and biophysical parameters in the assessment of soil salinity in the irrigated perimeter of Bahia, Brazil. J. S. Am. Earth Sci.; 2021; 112, 103580. [DOI: https://dx.doi.org/10.1016/j.jsames.2021.103580]

9. Lourenço, V.R.; Montenegro, A.A.A.; Carvalho, A.A.; Sousa, L.B.; Almeida, T.A.B.; Almeida, T.F.S.; Vilar, B.P. Spatial variability of biophysical multispectral indexes under heterogeneity and anisotropy for precision monitoring. Rev. Bras. Eng. Agrícola Ambient.; 2023; 27, pp. 848-857. [DOI: https://dx.doi.org/10.1590/1807-1929/agriambi.v27n11p848-857]

10. Huang, S.; Huang, S.; Tang, L.; Hupy, J.P.; Wang, Y.; Shao, G. A commentary review on the use of normalized difference vegetation index (NDVI) in the era of popular remote sensing. J. For. Res.; 2021; 32, pp. 1-6. [DOI: https://dx.doi.org/10.1007/s11676-020-01155-1]

11. Saponaro, M.; Agapiou, A.; Hadjimitsis, D.G.; Tarantino, E. Influence of Spatial Resolution for Vegetation Indices’ Extraction Using Visible Bands from Unmanned Aerial Vehicles’ Orthomosaics Datasets. Remote Sens.; 2021; 13, 3238. [DOI: https://dx.doi.org/10.3390/rs13163238]

12. Gerardo, R.; de Lima, I.P. Applying RGB-Based Vegetation Indices Obtained from UAS Imagery for Monitoring the Rice Crop at the Field Scale: A Case Study in Portugal. Agriculture; 2023; 13, 1916. [DOI: https://dx.doi.org/10.3390/agriculture13101916]

13. Atanasov, A.I.; Mihova, G.M.; Atanasov, A.Z.; Vlădut, V. Long-Term Assessment of NDVI Dynamics in Winter Wheat (Triticum aestivum) Using a Small Unmanned Aerial Vehicle. Agriculture; 2025; 15, 394. [DOI: https://dx.doi.org/10.3390/agriculture15040394]

14. Chen, Y.; Huang, L.; Chen, C.; Xie, Y. Vegetation growth monitoring based on ground-based visible light images from different views. Front. Environ. Sci.; 2025; 12, 1439045. [DOI: https://dx.doi.org/10.3389/fenvs.2024.1439045]

15. Bandak, S.; Movahedi-Naeini, S.A.; Mehri, S.; Lotfata, A. A longitudinal analysis of soil salinity changes using remotely sensed imageries. Sci. Rep.; 2024; 14, 10383. [DOI: https://dx.doi.org/10.1038/s41598-024-60033-6]

16. Abbas, J.A.A. Soil salinity assessment by using spectral salinity indices in Al-Sweira project middle of the Iraqi alluvial plain. Int. J. Environ. Sci. Technol.; 2023; 20, pp. 10847-10860. [DOI: https://dx.doi.org/10.1007/s13762-022-04733-4]

17. Wu, W. The generalized difference vegetation index (GDVI) for dryland characterization. Remote Sens.; 2014; 6, pp. 1211-1233. [DOI: https://dx.doi.org/10.3390/rs6021211]

18. Santos, W.M.d.; Martins, L.D.C.d.S.; Bezerra, A.C.; Souza, L.S.B.d.; Jardim, A.M.d.R.F.; Silva, M.V.d.; Souza, C.A.A.d.; Silva, T.G.F.d. Use of Unmanned Aerial Vehicles for Monitoring Pastures and Forages in Agricultural Sciences: A Systematic Review. Drones; 2024; 8, 585. [DOI: https://dx.doi.org/10.3390/drones8100585]

19. Ferraz, E.X.L.; Bezerra, A.C.; Lira, R.M.d.; Cruz Filho, E.M.d.; Santos, W.M.d.; Oliveira, H.F.E.d.; Silva, J.A.O.S.; Silva, M.V.d.; Silva, J.R.I.d.; Silva, J.L.B.d. . What Is the Predictive Capacity of Sesamum indicum L. Bioparameters Using Machine Learning with Red–Green–Blue (RGB) Images?. AgriEngineering; 2025; 7, 64. [DOI: https://dx.doi.org/10.3390/agriengineering7030064]

20. Melo, C.G.B.d.; Rolim, M.M.; Cavalcanti, R.Q.; Silva, M.V.d.; Candeias, A.L.B.; Lopes, P.M.O.; Ortiz, P.F.S.; Lima, R.P.d. RGB Composition Obtained by a UAV in the Monitoring of Sugarcane Row Gaps Using the Biophysical Index. AgriEngineering; 2025; 7, 17. [DOI: https://dx.doi.org/10.3390/agriengineering7010017]

21. Du, L.; Yang, H.; Song, X.; Wei, N.; Yu, C.; Wang, W.; Zhao, Y. Estimating leaf area index of maize using UAV-based digital imagery and machine learning methods. Sci. Rep.; 2022; 12, 15937. [DOI: https://dx.doi.org/10.1038/s41598-022-20299-0]

22. Feng, H.; Tao, H.; Li, Z.; Yang, G.; Zhao, C. Comparison of UAV RGB Imagery and Hyperspectral Remote-Sensing Data for Monitoring Winter Wheat Growth. Remote Sens.; 2022; 14, 3811. [DOI: https://dx.doi.org/10.3390/rs14174158]

23. Feng, L.; Khalil, U.; Aslam, B.; Ghaffar, B.; Tariq, A.; Jamil, A.; Farhan, M.; Aslam, M.; Soufan, W. Evaluation of Soil Texture Classification from Orthodox Interpolation and Machine Learning Techniques. Environ. Res.; 2024; 246, 118075. [DOI: https://dx.doi.org/10.1016/j.envres.2023.118075]

24. Marinho, G.T.B.; Pandorfi, H.; da Silva, M.V.; Montenegro, A.A.d.A.; de Sousa, L.d.B.; Desenzi, R.; da Silva, J.L.B.; de Oliveira-Júnior, J.F.; Mesquita, M.; de Almeida, G.L.P. . Bioclimatic Zoning for Sheep Farming through Geostatistical Modeling in the State of Pernambuco, Brazil. Animals; 2023; 13, 1124. [DOI: https://dx.doi.org/10.3390/ani13061124]

25. Almeida, T.A.B.; Montenegro, A.A.A.; de Lima, J.L.M.P.; Farias, C.W.L.A.; Carvalho, A.A.; de Paiva, A.L.R. Spatiotemporal Variability of Groundwater Quality for Irrigation: A Case Study in Mimoso Alluvial Valley, Semiarid Region of Brazil. Water; 2025; 17, 410. [DOI: https://dx.doi.org/10.3390/w17030410]

26. Pessoa, L.G.M.; Freire, M.B.G.d.S.; Green, C.H.M.; Miranda, M.F.A.d.A.; Filho, J.C.G.; Pessoa, W.R.L.S. Assessment of soil salinity status under different land-use conditions in the semiarid region of Northeastern Brazil. Ecol. Indic.; 2022; 141, 109139. [DOI: https://dx.doi.org/10.1016/j.ecolind.2022.109139]

27. Marengo, J.A.; Alves, L.M.; Alvala, R.C.S.; Cunha, A.P.; Brito, S.; Moraes, O.L.L. Climatic characteristics of the 2010–2016 drought in the semiarid northeast Brazil region. An. Acad. Bras. Ciênc.; 2018; 90, pp. 1973-1985. [DOI: https://dx.doi.org/10.1590/0001-3765201720170206]

28. Medeiros, S.S.; Salcedo, I.H.; Santos, D.B.; Batista, R.O.; Júnior, J.A.S.; Lima, R.C.C.; Marin, A.M.P. Esgotamento Sanitário: Panorama Para o Semiárido Brasileiro; Instituto Nacional do Semiárido—INSA: Campina Grande, Paraíba, 2014.

29. Menezes, J.A.; Madureira, A.P.; Santos, R.B.D.; Duval, I.B.; Regoto, P.; Margonari, C.; Barata, M.M.L.; Confalonieri, U. Analyzing Spatial Patterns of Health Vulnerability to Drought in the Brazilian Semiarid Region. Int. J. Environ. Res. Public Health; 2021; 18, 6262. [DOI: https://dx.doi.org/10.3390/ijerph18126262] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34207877]

30. Mayer, M.C.; Medeiros, S.S.; Batista, M.M.; Barbosa, R.A.; Lambais, G.R.; Santos, S.L.; Van Haandel, A. Tratamento de esgoto na zona rural visando ao reúso agrícola no semiárido brasileiro. Rev. DAE; 2021; 69, pp. 104-114. [DOI: https://dx.doi.org/10.36659/dae.2021.023]

31. Azevedo Júnior, M.S.; Ferreira Neto, M.; Medeiros, J.F.; Sá, F.V.S.; Lima, Y.B.; Lemos, M.; Queiroz, J.L.F.; Batista, R.O. Growth and Biomass production of prickly pear in the second cycle irrigated with treated domestc sewage. Biosci. J.; 2020; 36, pp. 51-60. [DOI: https://dx.doi.org/10.14393/BJ-v36n1a2020-42175]

32. Mishra, S.; Kumar, R.; Kumar, M. Use of treated sewage or wastewater as an irrigation water for agricultural purposes-Environmental, health, and economic impacts. Total. Environ. Res. Themes; 2023; 6, 100051. [DOI: https://dx.doi.org/10.1016/j.totert.2023.100051]

33. Sdiri, W.; AlSalem, H.S.; Al-Goul, S.T.; Binkadem, M.S.; Ben Mansour, H. Assessing the Effects of Treated Wastewater Irrigation on Soil Physico-Chemical Properties. Sustainability; 2023; 15, 5793. [DOI: https://dx.doi.org/10.3390/su15075793]

34. Santos, M.R.; Araújo, J.S.; Donato, S.L.R.; de Souza, J.A.A.; de Lira, E.C.; Aspiazú, I. Forage Cactus Pear Cultivars Irrigated with Wastewater in a Semi-Arid Region. Water; 2024; 16, 2632. [DOI: https://dx.doi.org/10.3390/w16182632]

35. Hag Husein, H.; Lucke, B.; Bäumler, R.; Sahwan, W. A Contribution to Soil Fertility Assessment for Arid and Semi-Arid Lands. Soil Syst.; 2021; 5, 42. [DOI: https://dx.doi.org/10.3390/soilsystems5030042]

36. Ferraz, R.L.S.; Costa, P.S.; Dantas Neto, J.; Anjos, F.A.; Barreto, N.P.; Soares, L.S.; Magalhães, I.D.; Barbosa, M.A. Opuntia ficus-indica (L.) mill. (Cactaceae) in climate change scenarios and its potential for wastewater bioremediationin semi-arid regions: A systematic review and meta-analysis. Am. J. Exp. Agric.; 2017; 18, pp. 1-11. [DOI: https://dx.doi.org/10.9734/JEAI/2017/36730]

37. Cheng, M.; Jiao, X.; Shi, L.; Penuelas, J.; Kumar, L.; Nie, C.; Wu, T.; Liu, K.; Wu, W.; Jin, X. High-resolution crop yield and water productivity dataset generated using random forest and remote sensing. Sci. Data; 2022; 9, 641. [DOI: https://dx.doi.org/10.1038/s41597-022-01761-0] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/36271097]

38. Pizarro, S.; Pricope, N.G.; Figueroa, D.; Carbajal, C.; Quispe, M.; Vera, J.; Alejandro, L.; Achallma, L.; Gonzalez, I.; Salazar, W. . Implementing Cloud Computing for the Digital Mapping of Agricultural Soil Properties from High Resolution UAV Multispectral Imagery. Remote Sens.; 2023; 15, 3203. [DOI: https://dx.doi.org/10.3390/rs15123203]

39. Almeida, T.A.B.; Montenegro, A.A.d.A.; Mackay, R.; Montenegro, S.M.G.L.; Coelho, V.H.R.; Carvalho, A.A.d.; Silva, T.G.F.d. Hydrogeological trends in an alluvial valley in the Brazilian semiarid: Impacts of observed climate variables change and exploitation on groundwater availability and salinity. J. Hydrol. Reg. Stud.; 2024; 53, 101784. [DOI: https://dx.doi.org/10.1016/j.ejrh.2024.101784]

40. Empresa Brasileira De Pesquisa Agropecuária-Embrapa. Sistema Brasileiro de Classificação de Solos; 3rd ed. Empresa Brasileira de Pesquisa Agropecuária – Embrapa: Brasília, Brazil, 2013; 353p.

41. Allen, R.G.; Pereira, L.S.; Raes, D.; Smith, M. Crop Evapotranspiration: Guidelines for Computing Crop Water Requirements; Irrigation and Drainage Paper 56; United Nations Food and Agriculture Organization: Rome, Italy, 1998; 300p.

42. Teixeira, P.C.; Donagemma, G.K.; Fontana, A.; Teixeira, W.G. Manual de Métodos de Análise de Solo; Embrapa Solos: Rio de Janeiro, Brazil, 2017; 573p.

43. Yeomans, J.C.; Bremner, J.M. A rapid and precise method for routine determination of organic carbon in soil. Commun. Soil Sci. Plant Anal.; 1988; 19, pp. 1467-1476. [DOI: https://dx.doi.org/10.1080/00103628809368027]

44. Rouse, J.W.; Haas, R.H.; Schell, J.A.; Deering, D.W. Monitoring vegetation systems in the Great Plains with ERTS. Third Earth Resources Technology Satellite-1 Symposium; Freden, S.C.; Mercanti, E.P.; Becker, M. NASA: Washington, DC, USA, 1974; Volume I: Technical Presentations, NASA SP-351 pp. 309-317.

45. Khan, N.M.; Rastoskuev, V.V.; Sato, Y.; Shiozawa, S. Assessment of hydrosaline land degradation by using a simple approach of remote sensing indicators. Agric. Water Manag.; 2005; 77, pp. 96-109. [DOI: https://dx.doi.org/10.1016/j.agwat.2004.09.038]

46. Bendig, J.; Yu, K.; Aasen, H.; Bolten, A.; Bennertz, S.; Broscheit, J.; Gnyp, M.L.; Bareth, G. Combining UAV-based plant height from crop surface models, visible, and near infrared vegetation indices for biomass monitoring in barley. Int. J. Appl. Earth Obs. Geoinf.; 2015; 39, pp. 79-87. [DOI: https://dx.doi.org/10.1016/j.jag.2015.02.012]

47. Gitelson, A.A.; Kaufman, Y.J.; Stark, R.; Rundquist, D. Novel algorithms for remote estimation of vegetation fraction. Remote Sens. Environ.; 2002; 80, pp. 76-87. [DOI: https://dx.doi.org/10.1016/S0034-4257(01)00289-9]

48. Huete, A.R. A soil-adjusted vegetation index (SAVI). Remote Sens. Environ.; 1988; 25, pp. 295-309. [DOI: https://dx.doi.org/10.1016/0034-4257(88)90106-X]

49. Huete, A.R.; Liu, H.Q.; Batchily, K.; Leeuwen, W.V. A Comparison of Vegetation Indices over a Global Set of TM Images for EOS-MODIS. Remote Sens. Environ.; 1997; 59, pp. 440-451. [DOI: https://dx.doi.org/10.1016/S0034-4257(96)00112-5]

50. QGIS Development Team. QGIS Geographic Information System; Software Version 3.34.10-Prizren QGIS Association: Bern, Switzerland, 2021; Available online: https://www.qgis.org (accessed on 22 February 2025).

51. Warrick, A.W.; Nielsen, D.R. Spatial variability of soil physical properties in the field. Applications of Soil Physics; Hillel, D. Academic Press: New York, NY, USA, 1980; pp. 319-344.

52. Andrade, T.S.; Montenegro, S.M.G.L.; Montenegro, A.A.A.; Rodrigues, D.F. Variabilidade espaço-temporal da condutividade elétrica da água subterrânea na região semiárida de Pernambuco. Rev. Bras. Eng. Agrícola Ambient.; 2012; 16, pp. 496-504. [DOI: https://dx.doi.org/10.1590/S1415-43662012000500005]

53. Hoaglin, D.C.; Mosteller, F.; Tukey, J.W. Análise Exploratória de Dados: Técnicas Robustas; Edições Salamandra: Lisboa, Portugal, 1992; 446p.

54. Lê, S.; Josse, J.; Husson, F. FactoMineR: An R Package for Multivariate Analysis. J. Stat. Softw.; 2008; 25, pp. 1-18. [DOI: https://dx.doi.org/10.18637/jss.v025.i01]

55. Journel, A.G. Fundamentals of Geostatistics in Five Lessons; American Geophysical Union: Washington, DC, USA, 1989; ISBN 0-87590-708-3

56. Englund, E.; Sparks, A. Geo-EAS 1.2.1 user’s guide. EPA Report 60018-91/008; EPA-EMSL: Las Vegas, NV, USA, 1988.

57. Robertson, G.P. GS+: Geoestatistics for the Environmental Sciences—GS + User’s Guide; Gamma Design Software: Plainwell, MI, USA, 2004; 152p.

58. Cambardella, C.A.; Moorman, T.B.; Novak, J.M.; Parkin, T.B.; Karlen, D.L.; Turco, R.F.; Kornopka, A.E. Field scale variability of soils properties in central Iowa Soils. Soil Sci. Soc. Am. J.; 1994; 58, pp. 1501-1511. [DOI: https://dx.doi.org/10.2136/sssaj1994.03615995005800050033x]

59. Webster, R.; Oliver, M. Geostatistics for Environmental Scientists; 2nd ed. Wiley: Chichester, UK, 2007.

60. Goovaerts, P. Geostatistical modelling of uncertainty in soil science. Geoderma; 2001; 103, pp. 3-26. [DOI: https://dx.doi.org/10.1016/S0016-7061(01)00067-2]

61. Surfer; Version 13.6.618 Computer Software Golden Software: Colorado, CO, USA, 2016.

62. Montenegro, A.A.A.; Montenegro, S.M.G.L. Variabilidade espacial de classes de textura, salinidade e condutividade hidráulica de solos em planície aluvial. Rev. Bras. Eng. Agrícola Ambient.; 2006; 10, pp. 30-37. [DOI: https://dx.doi.org/10.1590/S1415-43662006000100005]

63. Guan, Y.; Grote, K.; Schott, J.; Leverett, K. Prediction of Soil Water Content and Electrical Conductivity Using Random Forest Methods with UAV Multispectral and Ground-Coupled Geophysical Data. Remote Sens.; 2022; 14, 1023. [DOI: https://dx.doi.org/10.3390/rs14041023]

64. Taghadosi, M.M.; Hasanlou, M. Trend analysis of soil salinity in different land cover types using Landsat time series data (case study Bakhtegan Salt Lake). Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci.; 2017; 42, pp. 251-257. [DOI: https://dx.doi.org/10.5194/isprs-archives-XLII-4-W4-251-2017]

65. Belenok, V.; Hebryn-Baidy, L.; Bielousova, N.; Gladilin, V.; Kryachok, S.; Tereshchenko, A.; Alpert, S.; Bodnar, S. Machine learning based combinatorial analysis for land use and land cover assessment in Kyiv City (Ukraine). J. Appl. Remote Sens.; 2023; 17, 014506. [DOI: https://dx.doi.org/10.1117/1.JRS.17.014506]

66. Fu, T.; Gao, H.; Liu, J. Comparison of Different Interpolation Methods for Prediction of Soil Salinity in Arid Irrigation Region in Northern China. Agronomy; 2021; 11, 1535. [DOI: https://dx.doi.org/10.3390/agronomy11081535]

67. Guo, L.; Fu, P.; Shi, T.; Chen, Y.; Zhang, H.; Meng, R.; Wang, S. Mapping field-scale soil organic carbon with unmanned aircraft system-acquired time series multispectral images. Soil Tillage Res.; 2020; 196, 104477. [DOI: https://dx.doi.org/10.1016/j.still.2019.104477]

68. Lundgren, W.J.C.; Silva, J.A.A.; Ferreira, R.L.C. A precisão da estimativa do erro da krigagem pela validação cruzada. Floresta Ambiente; 2017; 24, e00124114. [DOI: https://dx.doi.org/10.1590/2179-8087.124114]