Terrain Controls on NDVI Spatial Variability under Post-Harvest Conditions: A UAV-Based Geomorphometric and Machine Learning Approach in Mediterranean Croplands

1  Introduction

Soil erosion and land degradation represent major environmental challenges worldwide, threatening food security, hydrological stability and ecosystem functioning [1,2]. Loss of fertile topsoil reduces agricultural productivity, accelerates runoff generation, increases sediment loads and disrupts biogeochemical cycles [3,4]. These processes tend to be self-reinforcing: once soil structure is weakened, infiltration declines and erosion intensifies, leading to a progressive decline in soil quality and long-term landscape resilience [5,6].

Remote sensing has emerged as a fundamental tool for monitoring and quantifying soil degradation processes [7,8]. Satellite imagery enables the assessment of vegetation dynamics, bare soil exposure, surface roughness and hydrological pathways at regional scales, while advances in UAV remote sensing have introduced unprecedented levels of spatial detail [9]. Drones equipped with RGB, multispectral and thermal sensors allow detecting microtopographic features, rills, crusting patterns and vegetation stress with centimetric resolution, capturing erosion signatures that remain invisible to traditional satellite sensors [10]. Beyond their role as high-resolution mapping platforms, UAV-based remote sensing systems have become powerful tools for applied geomorphological analysis, allowing the explicit linking of surface processes, terrain structure and ecosystem response [11]. The combination of centimetric Digital Terrain Models and spectral indicators enables the characterization of geomorphological influence on hydrological connectivity, sediment redistribution and soil degradation at spatial scales directly relevant for field-based management [12,13]. In this context, we hypothesize that moderate R2 values may reflect the intrinsic complexity of NDVI response under post-harvest and management-dominated conditions, where human disturbance and soil properties exert a stronger influence than vegetation physiology [14–16].

Multi-scale approaches are particularly valuable, as geomorphological processes operate simultaneously across spatial scales, from micro-relief elements controlling infiltration and crusting to broader slope and convergence patterns governing runoff concentration and gully development [17,18]. Integrating UAV-derived geomorphometry with spectral indices therefore represents a key methodological advance for applied geomorphology, supporting both processes understanding and the spatial targeting of soil conservation measures. In Mediterranean environments, erosion rates are among the highest in Europe, driven by intense convective rainfall, steep slopes, shallow soils and recurrent drought [19,20]. These landscapes are highly vulnerable due to their climatic variability and the cumulative impact of soil loss over millennia [21,22]. Agricultural systems are particularly exposed: conventional tillage, bare-soil periods, soil compaction, and unsustainable management practices amplify runoff and sediment export, making Mediterranean croplands hotspots of degradation and a priority for monitoring and mitigation.

There is therefore an increasing need for data-driven approaches capable of disentangling the spatial drivers of soil erosion and identifying vulnerable areas before irreversible degradation occurs [23,24]. Artificial Intelligence (AI), and especially machine learning models, offer robust frameworks for integrating topographic, hydrological and spectral indicators, allowing the identification of persistent erosion drivers and early-stage stress signals (e.g., Gholami et al. [25], 2018; Arabameri et al. [26], 2020). These tools have the potential to support sustainable land management by highlighting where conservation practices would be most efficient and by quantifying the influence of terrain features on vegetation condition. However, the application of machine learning techniques does not necessarily lead to high predictive performance in all study areas, particularly in heterogeneous agricultural environments. In many cases, the strong intra-parcel variability, the complex interaction between topography, soil properties and management practices, and the direct influence of human activities such as tillage, traffic by heavy machinery or land levelling introduces non-linear and site-specific effects that are difficult to capture using spectral indices alone [27]. European Mediterranean agroecosystems are especially challenging in this regard due to the coexistence of diverse soil types, shallow and disturbed profiles, and highly variable microtopography [28,29]. Under such conditions, machine learning approaches often prove more effective as diagnostic tools for identifying influential relationships between controlling factors, rather than as purely predictive models aimed at maximizing performance metrics.

The objective of this study is to evaluate the capacity of UAV-derived geomorphometric variables to explain the spatial variability of NDVI in a Mediterranean cereal field affected by soil degradation, under conditions of high environmental and management-related complexity. The test area is located adjacent to a large gully system and was surveyed immediately after cereal harvest, when the soil surface was bare, compacted and extremely dry following the summer drought in Casabermeja (Málaga, southern Spain). Using Random Forest regression, the study aims to assess whether UAV-derived geomorphometric metrics across multiple spatial scales can be used as a diagnostic framework to analyze NDVI spatial variability as a surface spectral response under post-harvest conditions and to highlight the relative importance of terrain-related controls shaping surface spectral variability, supporting the identification of areas potentially affected by soil degradation processes. In this context, we hypothesize moderate R2 values and may reflect the intrinsic complexity of NDVI response under post-harvest and management-dominated conditions, where human disturbance and soil properties exert a stronger influence than vegetation physiology. Moreover, under post-harvest conditions, when vegetation cover is minimal, NDVI values should be largely controlled by soil surface characteristics, including moisture content, residue cover, surface roughness and microtopographic variability. In such contexts, spectral responses would primarily reflect soil-related processes rather than plant physiological status, providing an opportunity to analyse the spatial imprint of erosion, compaction and water redistribution on the soil surface.

2  Materials and Methods

2.1 Study Area

The study site is located near Casabermeja in Málaga (southern Spain; Fig. 1), within a transitional zone between the metamorphic reliefs of the Montes de Málaga and the calcareous Subbetic ranges [30]. The field lies on flysch materials composed of sandstones, marls and expansive clays, a geological setting well known for its susceptibility to erosion, landslides and rapid soil deformation during intense rainfall events. Local soils, mainly vertisols and cambisols [31], support a mixed agroecosystem of cereals, olives and almonds, but their shrink–swell behaviour creates significant instability when disturbed by tillage, road construction or other land-use pressures [32]. The Guadalmedina River, which drains the surrounding mountains, cuts deep valleys in these weak lithologies under a Mediterranean climate with mild temperatures and concentrated autumn–winter storms. These geological and climatic factors combine with human activities—paths, access cuts, nearby industrial areas and intensive cultivation—to accelerate degradation. The gully analysed in this study develops across a cereal field and exemplifies the region’s vulnerability, as expansive clays, steep microtopography and disturbed slopes favour rapid incision and soil losses [33]. This geomorphological and climatic context makes the study area particularly suitable for evaluating the potential of UAV-based geomorphometry to characterize soil degradation processes at field scale.

Figure 1: Study area. (a): Localisation in the Iberian Peninsula; (b): localisation in the Andalusian region; (c): orthomosaic with RGB camera; (d): NDVI index and study area for random forest training simulation; (e): digital terrain model (DTM); and; (f): 3D model and study area. Coordinate reference system used for research is ETRS89/UTM zone 30 N.

2.2 UAV Data Acquisition and Processing

UAV surveys were conducted over the study area using a DJI Mavic 3 Multispectral platform on 16 July 2025, coinciding with post-harvest conditions characterized by bare, compacted and extremely dry soils. The DJI Mavic 3 Multispectral platform integrates an RGB camera and a multispectral camera, enabling the acquisition of imagery for both photogrammetric terrain reconstruction and spectral index computation. Two independent flights were carried out: the RGB imagery was used for photogrammetric terrain reconstruction and DTM generation, while the multispectral imagery was used for reflectance mapping and NDVI computation. The RGB flight consisted of 360 calibrated images acquired at a ground sampling distance (GSD) of 2.54 cm, covering approximately 0.32 km2. Photogrammetric processing produced dense point clouds exceeding 50 million points (overlapping from 70%–80%), from which Digital Surface Models (DSM) and a Digital Terrain Model (DTM) were generated at native resolution. The multispectral flight comprised 1440 images from the green, red, red-edge and near-infrared bands, with a GSD of 4.27 cm over 0.29 km2. Radiometrically corrected reflectance maps and NDVI were generated using standard structure-from-motion and image calibration workflows. Twenty-four Ground Control Points (GCPs) were surveyed using an Emlid Reach RS3 GNSS receiver to improve the spatial accuracy and georeferencing of the UAV-derived products. The use of post-harvest conditions minimized vegetation interference in terrain reconstruction and maximized the sensitivity of NDVI to soil and microtopographic controls [34].

All raster products derived from RGB and multispectral processing were co-registered to a common spatial reference system (ETRS89/UTM Zone 30 N; EPSG: 25830) using Drone2Map (ESRI, USA). The Digital Terrain Model and geomorphometric layers were spatially aligned to the NDVI grid to ensure pixel-to-pixel correspondence. Resampling procedures preserved the original spatial resolution of each dataset. A cereal-field mask was applied to restrict the analysis to homogeneous land cover conditions and to exclude olive groves, roads and riparian areas that could introduce spectral or geomorphometric noise. The cereal-field mask was generated through manual digitization based on the UAV orthomosaic, allowing the exclusion of non-cereal elements such as olive groves, roads and riparian areas.

2.3 Geomorphometric Variable Generation

Geomorphometric predictors were derived from the UAV-based DTM to characterize terrain structure and hydrological behaviour at multiple spatial scales using ArcGis Pro (ESRI, USA). First-order terrain derivatives included slope and aspect, with aspect transformed into northness and eastness to represent directional solar exposure. Curvature metrics (general, planform and profile curvature) were computed to capture surface convexity, flow divergence and downslope acceleration. Hydrological and erosion-related indices were derived from flow-direction and flow-accumulation models, including the Topographic Wetness Index (TWI) and the Stream Power Index (SPI), which quantify potential moisture accumulation and erosive power, respectively. The Topographic Wetness Index (TWI) was calculated following the formulation of [35], while the Stream Power Index (SPI) was derived according to Moore et al. [36] (1991), representing potential soil moisture accumulation and erosive power, respectively. To explicitly account for scale-dependent geomorphological processes, a set of multi-scale local relief metrics was generated using neighbourhood radii of 1, 2, 4, 6 and 10 m. Multi-scale local relief metrics (lrelief2m, lrelief4m, lrelief6m and lrelief10m) were computed as the difference between the elevation of each pixel and the mean elevation within a circular neighbourhood of a given radius (2, 4, 6 and 10 m, respectively). These variables capture relative elevation deviations at different spatial scales, representing micro- to macro-topographic variability associated with soil redistribution, water retention and localized erosion processes. The variable “relief_local” corresponds to a generalized local relief metric calculated using a standard neighbourhood window, providing an overall indicator of terrain roughness at the plot scale. A total of 3000 sampling points were randomly distributed within the cereal-field mask (Fig. 2) to ensure representative coverage of slope gradients, curvature regimes and distances to the gully system. These points were generated using a random sampling approach within the masked area, avoiding spatial bias associated with regular grids. The selected sample size (n = 3000) represents a balance between statistical robustness and computational efficiency for Random Forest modelling at high spatial resolution. At each point, NDVI values and all geomorphometric predictors were extracted, generating a dataset comprising the response variable (NDVI) and 14 explanatory terrain variables.

Figure 2: Methodological workflow integrating UAV data acquisition, geomorphometric analysis and machine learning modelling.

2.4 Machine Learning Modelling

The relationship between geomorphometric variables and NDVI spatial variability was modelled using a Random Forest regression approach implemented in ArcGIS Pro (ESRI, USA). Model training and prediction were performed using the Forest-based Classification and Regression tool within the Spatial Statistics toolbox. Random Forest was selected due to its robustness to multicollinearity, its ability to capture non-linear relationships, and its suitability for high-dimensional environmental datasets. A total of 3000 samples were used, with a 90/10 training–validation split. The model was configured with 100 trees, a minimum leaf size of 10, and four randomly selected predictor variables at each split. Model performance was evaluated using the coefficient of determination (R2), root mean squared error (RMSE) and mean absolute error (MAE), as automatically reported by the ArcGIS Pro modelling framework. In addition to predictive accuracy, model interpretation focused on variable-importance metrics exported by the tool, which quantify the relative contribution of each geomorphometric predictor to NDVI variability. It should be noted that several geomorphometric predictors are inherently correlated, as they are derived from the same Digital Terrain Model. Therefore, variable-importance results should be interpreted with caution and considered as indicative rather than causal. This allowed the identification of terrain variables associated with NDVI spatial variability while maintaining a transparent and reproducible modelling workflow (Fig. 2).

3  Results

NDVI refers to the Normalized Difference Vegetation Index computed from NIR and red reflectance. Multi-scale local relief metrics (lrelief2m, lrelief4m, lrelief6m and lrelief10m) represent elevation deviations from the local mean calculated at increasing neighbourhood radii (2, 4, 6, and 10 m, respectively), while “relief_local” corresponds to general local relief relative to the surrounding terrain. TWI (Topographic Wetness Index) expresses the potential for soil moisture accumulation, and SPI (Stream Power Index) represents the erosive power of concentrated flow. Flowacc (flow accumulation) indicates the number of upstream contributing cells. Curvature (general curvature), pl_curvature (plan curvature) and p_curvature (profile curvature) describe surface shape and flow dynamics. Northness and eastness correspond to cosine- and sine-transformed aspect, representing directional solar exposure. Slope_deg represents slope in degrees.

3.1 Spatial Patterns of NDVI and Geomorphometric Variables

The statistical distribution of NDVI and geomorphometric predictors is summarized in Table 1. NDVI values range from −0.13 to 0.53, with a mean of 0.12 and a low variance (0.002), reflecting a relatively homogeneous soil status with localized patches of distinct spectral variability. Multi-scale local relief metrics (lrelief2mi, lrelief4me, lrelief6fi, lrelief10m and Relief local) show mean values close to zero and very small variances, which is consistent with their definition as deviations from local elevation means and indicates that most points lie close to locally neutral topographic positions, with only a subset associated with pronounced micro-relief features.

Hydrological and erosion-related variables present much larger dynamic ranges. Flow accumulation spans more than five orders of magnitude (0 to 570,021), with a mean of 562, while SPI ranges from 0 to almost 19,500, indicating the presence of strongly convergent flow paths and high-erosion potential zones within the field. Curvature metrics (general, plan and profile) are symmetrically distributed around zero with substantial variance, capturing the alternation between convex and concave landforms that control water convergence and divergence at the plot scale. Northness and eastness cover their full theoretical range [−1, 1], with slightly negative means, revealing a predominance of aspects oriented away from north and east and therefore a bias towards warmer, more exposed slopes. Slope is generally low (mean 2°) but reaches values above 40° in localized sectors close to the taluses and gully, confirming the presence of very steep micro-slopes associated with rill features that are likely to to influence soil moisture redistribution and surface conditions.

NDVI refers to the Normalized Difference Vegetation Index computed from NIR and Red reflectance. The multi-scale local relief variables (lrelief10m, lrelief6fi, lrelief4me and lrelief2mi) represent elevation deviations at macro-, fine-, meso- and micro-topographic scales, respectively, while “relief local” corresponds to general local relief relative to neighborhood mean elevation. TWI (Topographic Wetness Index; values of −9999 represent NoData cells and were removed before statistical analyses) expresses the potential for soil moisture accumulation, and SPI (Stream Power Index) quantifies the erosive power of concentrated flow. Flowacc denotes flow accumulation, representing the number of upstream contributing cells. Curvature indicates overall surface convexity or concavity, whereas planform and profile curvatures (pl_curvatu and p_curvatur) describe lateral flow divergence and downslope flow acceleration, respectively. Northness and eastness correspond to cosine- and sine-transformed aspect, capturing directional exposure to sunlight.

3.2 Random Forest Model Performance

The Random Forest regression achieved moderate predictive skill for NDVI within the cereal field, in line with previous studies on UAV-based NDVI spatial modelling under post-harvest conditions in complex terrain. The high spatial resolution of the NDVI dataset (4.27 cm) partly explains the moderate predictive performance, as fine-scale spectral variability introduces noise that is difficult to capture using terrain-based predictors alone. The model produced an R2 of 0.22 for the training data and 0.31 for the validation subset, showing a consistent but limited ability to capture fine-scale NDVI variability. Prediction errors remained stable, with a Root Mean Squared Error (RMSE) of 0.04 and Mean Absolute Error (MAE) of 0.03, indicating that the model reproduced broad NDVI gradients but not localized spectral noise or fine-scale variability in surface reflectance associated with soil conditions, residue cover and microtopography A summary of model diagnostics is presented in Table 2, showing similar error magnitudes between training and validation sets. Although the explained variance was low (0.00), this behaviour is expected in high-resolution NDVI modelling where spectral responses are influenced by a combination of degraded soils by heavy machinery, plant absence, illumination and microtopographic factors that are not fully captured by terrain-based predictors alone.

3.3 Variable Importance

Fig. 3 shows the spatial distribution of the main geomorphometric predictors derived from the UAV-based DTM. Fig. 4 presents the Random Forest variable-importance ranking used to identify the terrain variables most strongly associated with observed NDVI spatial variability. Slope-related metrics emerged as the most influential group of predictors within the model, representing nearly 29% of total importance, highlighting the role of gravitational water redistribution and runoff acceleration associated with lower NDVI values in steep Mediterranean agroecosystems. Eastness ranked second (≈12%), reflecting strong morning–afternoon radiation asymmetries that shape evapotranspiration and soil temperature patterns.

Figure 3: Spatial distribution of the main geomorphometric predictors derived from the UAV-based digital terrain model (DTM), including multi-scale local relief metrics (lrelief10m, lrelief6fine, lrelief4me, lrelief2mi and relief_local), hydrological indices (stream power index, SPI; topographic wetness index, TWI), flow accumulation (flowacc), and curvature metrics (general and profile curvature).

Figure 4: Relative importance of geomorphometric predictors in the random forest regression model used to explain NDVI spatial variability (model output). Importance values represent the relative contribution of each variable to the model. Slope_deg denotes slope steepness (degrees); eastness and northness are cosine- and sine-transformed aspect components representing solar exposure; SPI is the Stream Power Index; flowacc is flow accumulation; TWI is the topographic wetness index; curvature, pl_curvature and p_curvature correspond to general, planform and profile curvature, respectively; lrelief2mi, lrelief4me, lrelief6fi and lrelief10m represent multi-scale local relief metrics computed at increasing neighbourhood radii; relief_loc denotes general local relief.

Erosion-related metrics, including the Stream Power Index (SPI) and Flow Accumulation, contributed significantly (≈13% combined), indicating that areas with concentrated flow and potential soil removal displayed lower modelled NDVI (Random Forest output). Multi-scale local relief variables collectively accounted for more than 20% of total importance, demonstrating that subtle variations in micro- and meso-topography strongly influence moisture retention, root-zone structure and sediment deposition. Curvature metrics showed moderate but consistent effects, characterizing convex and concave terrain elements associated with drainage efficiency. The complete ranking of predictor importance appears in Table 3 and is visually summarised in Fig. 3.

The results highlight three main mechanisms linking geomorphology and NDVI spatial variability in post-harvest conditions. First, slope-driven water redistribution emerged as the primary limiting factor, with steeper sectors showing reduced observed NDVI due to lower infiltration and higher runoff and soil-loss potential. Second, solar radiation asymmetry, captured through eastness (and to a lesser degree northness), influenced spatial NDVI patterns by driving differential evapotranspiration and microclimatic contrasts across the field. Third, local terrain structure—reflected by relief layers and curvature—modulated the distribution of shallow soil pockets, nutrient accumulation, organic-matter retention and micro-drainage pathways. Hydrological and erosive indices (SPI and FlowAcc) reinforced the connection between reduced NDVI values and topographic convergence zones, particularly near the gully margin where erosive forces degrade soil structure. The combination of hydrological, radiative and multi-scale relief variables provides a coherent explanation of the observed NDVI distribution, despite the modest predictive performance of the model.

Based on the combined analysis of NDVI and geomorphometric variables, areas characterized by low NDVI values, high slope gradients and elevated flow-convergence indices (e.g., SPI and flow accumulation) can be interpreted as zones with higher potential soil degradation. These areas are mainly located in proximity to the gully system and along sectors with pronounced microtopographic variability, where erosive processes and soil redistribution are more active. This integrated interpretation supports the use of UAV-derived geomorphometry and spectral indicators as a diagnostic framework for identifying degradation-prone areas at field scale.

4  Discussion

In this study, NDVI patterns should therefore be interpreted as a surface spectral response reflecting soil conditions rather than vegetation dynamics. This interpretation is particularly relevant under post-harvest Mediterranean conditions, where bare soils expose the influence of geomorphological and hydrological processes on surface reflectance. The results highlight three main mechanisms linking geomorphology and surface spectral response within the studied Mediterranean cereal field. First, slope-driven water redistribution emerged as the primary limiting factor, with steeper sectors showing reduced NDVI values due to lower infiltration capacity, enhanced runoff generation and higher soil-loss potential. These conditions directly affect soil structure and moisture availability, limiting NDVI spatial variability under post-harvest conditions and extremely dry conditions [37,38]. Second, solar radiation asymmetry, captured through eastness and, to a lesser degree, northness, exerted a clear influence on NDVI spatial patterns. Differences in morning and afternoon radiation exposure drive contrasts in evapotranspiration rates and soil temperature, which are particularly relevant in Mediterranean environments characterized by prolonged summer droughts [39,40]. These radiative controls contribute to persistent spatial patterns of vegetation stress that remain detectable even when vegetation cover is minimal. Third, local terrain structure, reflected by multi-scale relief layers and curvature metrics, modulated the distribution of shallow soil pockets, sediment accumulation zones and micro-drainage pathways [41,42].

Hydrological and erosive indices, such as the Stream Power Index and flow accumulation, reinforced the link between surface spectral variability and zones of topographic convergence, particularly near the gully margins where erosive forces degrade soil structure [43]. The combined influence of hydrological, radiative and multi-scale relief variables provides a coherent explanation for the observed NDVI distribution, despite the modest predictive performance of the model [44].

The multi-scale local relief metrics provide insight into the hierarchical organization of geomorphological processes [18,45]. Small neighbourhood radii (1–2 m) capture microtopographic variability related to surface roughness, crusting and incipient rill formation. Intermediate scales (4–6 m) reflect localized runoff redistribution and the development of small drainage features, while larger scales (10 m) represent broader slope organization and flow convergence patterns. The higher importance of larger-scale metrics suggests that NDVI spatial variability is more strongly associated with terrain features controlling water redistribution at the field scale rather than with very fine-scale surface irregularities.

From a modelling perspective, the Random Forest regression achieved moderate predictive skill (R2 = 0.31), which is consistent with previous UAV-based NDVI modelling studies conducted under complex terrain and weak vegetation signal conditions. The limited explained variance reflects the fact that NDVI is influenced by multiple factors beyond topography, including soil compaction by heavy machinery, surface brightness effects, illumination variability and microtopographic noise that are not fully captured by terrain-based predictors alone. Nevertheless, the structure of variable importance revealed robust and interpretable geomorphological factors associated with NDVI spatial variability, highlighting the diagnostic value of machine learning approaches beyond pure prediction. To support the interpretation of the results and to synthesize the main terrain–NDVI relationships identified in this study, a conceptual framework is proposed linking UAV-derived geomorphometric association, eco-hydrological processes and NDVI response under post-harvest Mediterranean conditions (Fig. 5). We emphasize that the moderate predictive performance (R2 = 0.31) reflects the inherently complex and multi-factor nature of NDVI under post-harvest conditions, where the spectral signal is influenced not only by geomorphometry but also by soil properties, residue cover, surface crusting, illumination conditions and management practices. These factors introduce spatial variability that cannot be fully captured by terrain-based predictors alone. Therefore, the model should be interpreted as a diagnostic tool highlighting dominant spatial patterns rather than as a predictive model aimed at maximizing accuracy. The slightly higher validation R2 compared to the training R2 reflects the spatial heterogeneity of the dataset and the effects of random sampling in a high-resolution, noise-dominated environment. Under these conditions, small differences in sample composition between training and validation subsets can lead to variability in performance metrics, which should be interpreted with caution.

Figure 5: Conceptual framework illustrating the relationships between UAV-derived geomorphometric influence, eco-hydrological processes and modelled NDVI (random forest output) response under Mediterranean post-harvest conditions.

The modelling framework presented here can be expanded through several methodological directions. The comparison of additional machine learning models beyond Random Forest, such as Gradient Boosted Trees (e.g., XGBoost, LightGBM), Support Vector Regression, k-Nearest Neighbors regression, Neural Networks, hybrid CNN–RF models or Gaussian Process Regression for uncertainty quantification, may yield higher predictive performance or complementary insights [46]. Benchmarking these approaches would help identify optimal modelling strategies for UAV-based agronomic and geomorphological applications.

The primary objective of this study was to analyse the relationship between geomorphometric variables and NDVI spatial variability, rather than to produce a direct soil degradation map. However, the spatial distribution of geomorphometric predictors shown in Fig. 3 supports the interpretation that areas with low NDVI values, particularly those associated with steep slopes and flow-convergence zones, may correspond to zones with higher degradation potential. This interpretation should be understood as a diagnostic assessment based on NDVI patterns and geomorphometric controls, not as a direct mapping of soil degradation. We emphasize that this representation is intended as a diagnostic interpretation based on NDVI patterns and geomorphometric influence, rather than a direct mapping of soil degradation. Areas exhibiting low NDVI values, particularly in steep sectors and near flow-convergence zones, can be interpreted as locations with higher potential soil degradation. In addition to natural geomorphological controls, human disturbance plays a significant role in shaping surface spectral variability in Mediterranean agricultural systems [47,48]. Tillage operations and repeated machinery traffic generate regular geometric patterns, such as ploughing lines and tractor tracks, which alter soil structure, compaction and surface roughness [49,50]. These anthropogenic features can introduce spatial variability in NDVI that is not directly related to terrain attributes, potentially interfering with geomorphometric indices and contributing to the moderate predictive performance observed. This effect is particularly relevant in high-resolution UAV datasets, where fine-scale management-induced patterns become highly visible and may partially mask or interact with natural topographic signals. Moreover, future analyses should incorporate variables beyond topography, including soil electrical conductivity, soil texture maps, high-resolution soil moisture estimates, canopy height models derived from RGB photogrammetry and thermal imagery for stress detection. Multi-date modelling approaches would allow the evaluation of NDVI–geomorphometry relationships across different stages of the crop cycle, including pre-harvest, post-rainfall events, mid-season drought periods or after extreme rainfall events. Such temporal analyses would enable the assessment of the persistence of stress patterns, the role of slope–aspect combinations in heat accumulation and the evolution of erosion impacts before and after prolonged dry periods, strengthening the understanding of climate-driven crop responses. Future work should also explore the integration of UAV-derived data with coarser-resolution satellite imagery (e.g., Sentinel-2) and alternative vegetation indices to evaluate the scalability of the observed terrain–NDVI relationships.

5  Conclusions

This study shows that the combination of UAV-derived geomorphometry and machine learning provides an effective framework for analysing spatial patterns of NDVI variability and soil degradation in Mediterranean cereal systems affected by erosion, compaction and drought. Even under post-harvest conditions characterised by weak spectral signal dominated by soil surface conditions, terrain-related factors exerted a consistent control on NDVI spatial variability. Slope and solar exposure were identified as the most influential geomorphometric factors, while hydrological indices and multi-scale local relief metrics highlighted the influence of runoff convergence, microtopography and soil structure heterogeneity. These results confirm the importance of explicitly accounting for scale-dependent geomorphological processes when interpreting terrain–vegetation relationships in heterogeneous agricultural environments. From a methodological perspective, the moderate predictive performance of the Random Forest model reflects the intrinsic complexity of NDVI under dry and managed conditions, rather than a limitation of the approach. The analysis demonstrates the value of machine learning as a diagnostic tool for identifying associated controlling factors, rather than as a purely predictive model. Despite limitations related to the use of a single acquisition date and the absence of explicit soil and management variables, the proposed framework provides a robust baseline for future multi-temporal and multi-variable UAV-based studies of soil degradation in Mediterranean agroecosystems.

Acknowledgement: During the preparation of this manuscript, generative AI tools (ChatGPT and Gemini) were used exclusively for language refinement and for assisting in the drafting of schematic figures (e.g., workflow diagrams). All outputs were critically reviewed, modified and validated by the authors, who take full responsibility for the final content.

Funding Statement: This research was funded by the Spanish Ministerio de Ciencia, Innovación y Universidades, under the project titled “Desarrollo de productos basados en los nuevos sensores satelitales hiperespectrales europeos e ia para la caracterización de estresores en tierras de cultivo (HIPROESTRES)” (grant number PID2023-152656OB-I00), within the Programa Estatal de Investigación Científica, Técnica y de Innovación (2021–2023).

Author Contributions: Jesús Rodrigo-Comino: Conceptualization, Data curation, Formal analysis, Project administration, Funding acquisition, Visualization, Writing—original draft. María Teresa González-Moreno: Data curation, Formal analysis, Visualization, Writing—review & editing. Lucía Moreno-Cuenca: Writing—review & editing. Laura Cambronero-Ruiz: Writing—review & editing. Clemente Irigaray: Writing—review & editing. Francisco Serrano Bernardo: Writing—review & editing. Víctor Hugo Durán Zuazo: Writing—review & editing. Jesús Fernández-Gálvez: Writing—review & editing. Andrés Caballero-Calvo: Writing—review & editing. Víctor Rodríguez-Galiano: Conceptualization, Project administration, Resources, Writing—review & editing. All authors reviewed and approved the final version of the manuscript.

Availability of Data and Materials: Data and processing workflows are available from the corresponding author upon reasonable request.

Ethics Approval: Not applicable.

Conflicts of Interest: The authors declare no conflicts of interest.

References

1. An Y, Zhao W, Li C, Sofia Santos Ferreira C. Temporal changes on soil conservation services in large basins across the world. Catena. 2022;209(D5):105793. doi:10.1016/j.catena.2021.105793. [Google Scholar] [CrossRef]

2. Chakrabortty R, Ali T, Pal T, Pande CB, Elaksher AF, Abioui M. Climate change and land use dynamics: modeling soil erosion scenarios to achieve sustainable development goals. Earth Syst Environ. 2026;10(1):749–74. doi:10.1007/s41748-025-00631-0. [Google Scholar] [CrossRef]

3. Borrelli P, Ballabio C, Yang JE, Robinson DA, Panagos P. GloSEM: high-resolution global estimates of present and future soil displacement in croplands by water erosion. Sci Data. 2022;9(1):406. doi:10.1038/s41597-022-01489-x. [Google Scholar] [PubMed] [CrossRef]

4. Quinton JN, Fiener P. Soil erosion on arable land: an unresolved global environmental threat. Prog Phys Geogr Earth Environ. 2024;48(1):136–61. doi:10.1177/03091333231216595. [Google Scholar] [CrossRef]

5. Poesen J. Soil erosion in the anthropocene: research needs. Earth Surf Processes Landf. 2018;43(1):64–84. doi:10.1002/esp.4250. [Google Scholar] [CrossRef]

6. García-Ruiz JM, Nadal-Romero E, Lana-Renault N, Beguería S. Erosion in Mediterranean landscapes: changes and future challenges. Geomorphology. 2013;198(4):20–36. doi:10.1016/j.geomorph.2013.05.023. [Google Scholar] [CrossRef]

7. Arabameri A, Cerda A, Rodrigo-Comino J, Pradhan B, Sohrabi M, Blaschke T, et al. Proposing a novel predictive technique for gully erosion susceptibility mapping in arid and semi-arid regions (Iran). Remote Sens. 2019;11(21):2577. doi:10.3390/rs11212577. [Google Scholar] [CrossRef]

8. Kaiser A, Neugirg F, Rock G, Müller C, Haas F, Ries J, et al. Small-scale surface reconstruction and volume calculation of soil erosion in complex Moroccan gully morphology using structure from motion. Remote Sens. 2014;6(8):7050–80. doi:10.3390/rs6087050. [Google Scholar] [CrossRef]

9. Riddle L, Hill T, Gijsbertsen B. Geomorphology from ‘on high’: the use of drones/UAV technology in teaching soil erosion. JoGEA. 2021;2021(4):58–74. doi:10.46622/jogea.v4i1.3806. [Google Scholar] [CrossRef]

10. Peter KD, d’Oleire-Oltmanns S, Ries JB, Marzolff I, Ait Hssaine A. Soil erosion in gully catchments affected by land-levelling measures in the Souss Basin, Morocco, analysed by rainfall simulation and UAV remote sensing data. Catena. 2014;113(4):24–40. doi:10.1016/j.catena.2013.09.004. [Google Scholar] [CrossRef]

11. Alexiou S, Efthimiou N, Karamesouti M, Papanikolaou I, Psomiadis E, Charizopoulos N. Measuring annual sedimentation through high accuracy UAV-photogrammetry data and comparison with RUSLE and PESERA erosion models. Remote Sens. 2023;15(5):1339. doi:10.3390/rs15051339. [Google Scholar] [CrossRef]

12. Agarwal S, Sundriyal Y, Srivastava P. Dam in Himalaya induces geomorphic disconnectivity during extreme hydrological event: evaluating a case of 2013 Kedarnath Disaster. J Earth Syst Sci. 2022;131(4):263. doi:10.1007/s12040-022-01991-1. [Google Scholar] [CrossRef]

13. Keesstra S, Nunes JP, Saco P, Parsons T, Poeppl R, Masselink R, et al. The way forward: can connectivity be useful to design better measuring and modelling schemes for water and sediment dynamics? Sci Total Environ. 2018;644(10):1557–72. doi:10.1016/j.scitotenv.2018.06.342. [Google Scholar] [PubMed] [CrossRef]

14. de Castro AI, Shi Y, Maja JM, Peña JM. UAVs for vegetation monitoring: overview and recent scientific contributions. Remote Sens. 2021;13(11):2139. doi:10.3390/rs13112139. [Google Scholar] [CrossRef]

15. Douglas TJ, Coops NC, Drever MC. UAV-acquired imagery with photogrammetry provides accurate measures of mudflat elevation gradients and microtopography for investigating microphytobenthos patterning. Sci Remote Sens. 2023;7(10):100089. doi:10.1016/j.srs.2023.100089. [Google Scholar] [CrossRef]

16. Wang J, Huang H, Ariyasena HHS, Zhao J, Zhang X, Gao X, et al. A UAV-based method for root zone soil moisture modeling of different farmland scale with grain and economic crops. Agric Water Manag. 2025;321(9):109932. doi:10.1016/j.agwat.2025.109932. [Google Scholar] [CrossRef]

17. Cerdà A, Flanagan DC, le Bissonnais Y, Boardman J. Soil erosion and agriculture. Soil Tillage Res. 2009;106(1):107–8. doi:10.1016/j.still.2009.10.006. [Google Scholar] [CrossRef]

18. Tarolli P. High-resolution topography for understanding earth surface processes: opportunities and challenges. Geomorphology. 2014;216:295–312. doi:10.1016/j.geomorph.2014.03.008. [Google Scholar] [CrossRef]

19. Cerdan O, Govers G, Le Bissonnais Y, Van Oost K, Poesen J, Saby N, et al. Rates and spatial variations of soil erosion in Europe: a study based on erosion plot data. Geomorphology. 2010;122(1–2):167–77. doi:10.1016/j.geomorph.2010.06.011. [Google Scholar] [CrossRef]

20. García-Ruiz JM, Beguería S, Nadal-Romero E, González-Hidalgo JC, Lana-Renault N, Sanjuán Y. A meta-analysis of soil erosion rates across the world. Geomorphology. 2015;239(1–2):160–73. doi:10.1016/j.geomorph.2015.03.008. [Google Scholar] [CrossRef]

21. Amate JI, de Molina MG, Vanwalleghem T, Fernández DS, Gómez JA. Erosion in the Mediterranean: the case of olive groves in the south of Spain (1752–2000). Environ Hist. 2013;18(2):360–82. doi:10.1093/envhis/emt001. [Google Scholar] [CrossRef]

22. Bagarello V, Ferro V. Scale effects on plot runoff and soil erosion in a Mediterranean environment. Vadose Zone J. 2017;16(12):vzj2017.03.0059. doi:10.2136/vzj2017.03.0059. [Google Scholar] [CrossRef]

23. Busico G, Grilli E, Carvalho SCP, Mastrocicco M, Castaldi S. Assessing soil erosion susceptibility for past and future scenarios in semiarid Mediterranean agroecosystems. Sustainability. 2023;15(17):12992. doi:10.3390/su151712992. [Google Scholar] [CrossRef]

24. Diodato N, Fiorillo F, Rinaldi M, Bellocchi G. Environmental drivers of dynamic soil erosion change in a Mediterranean fluvial landscape. PLoS One. 2022;17(1):e0262132. doi:10.1371/journal.pone.0262132. [Google Scholar] [PubMed] [CrossRef]

25. Gholami V, Booij MJ, Nikzad Tehrani E, Hadian MA. Spatial soil erosion estimation using an artificial neural network (ANN) and field plot data. Catena. 2018;163(1):210–8. doi:10.1016/j.catena.2017.12.027. [Google Scholar] [CrossRef]

26. Arabameri A, Cerda A, Pradhan B, Tiefenbacher JP, Lombardo L, Bui DT. A methodological comparison of head-cut based gully erosion susceptibility models: combined use of statistical and artificial intelligence. Geomorphology. 2020;359(7):107136. doi:10.1016/j.geomorph.2020.107136. [Google Scholar] [CrossRef]

27. Rodrigo-Comino J, Cambronero-Ruiz L, Moreno-Cuenca L, González-Vivar J, González-Moreno MT, Rodríguez-Galiano V. Integrating UAV-LiDAR and field experiments to survey soil erosion drivers in Citrus Orchards using an exploratory machine learning approach. Water. 2025;17(24):3541. doi:10.3390/w17243541. [Google Scholar] [CrossRef]

28. Alewell C, Gupta S, Poulenard J, Niquille N, Kaiser A, Shokri N, et al. A first quantitative assessment of soil health at European scale considering soil genesis. J Plant Nutr Soil Sci. 2026;189(1):6–16. doi:10.1002/jpln.70034. [Google Scholar] [CrossRef]

29. Robinson DA, Fendrich A, Thomas A, Reinsch S, Leifeld J, Moore T, et al. Soil porosity prediction across Europe with a focus on soil particle density determination. Int Soil Water Conserv Res. 2026;14(2):100614. doi:10.1016/j.iswcr.2026.100614. [Google Scholar] [CrossRef]

30. Rodrigo-Comino J, Cambronero-Ruiz L, Gonzalez-Moreno MT, Keesstra SD, Vivar JG, Cerda A, et al. Sustainable development and integrated soil management: tools and keys for surveying the landscape evolution and degradation processes in Mediterranean mountainous regions. In: Environmental sustainability and global change. Amsterdam, The Netherlands: Elsevier; 2025. p. 19–31. doi:10.1016/b978-0-443-31596-1.00002-7. [Google Scholar] [CrossRef]

31. Rodrigo-Comino J. Los Suelos de la provincia de Málaga: revisión y actualización de las fuentes edafológicas según La clasificación de FAO-WRB. Málaga, Spain: Servicio de Publicaciones y Divulgación Científica, Universidad de Málaga; 2014. (In Spanish). [Google Scholar]

32. Senciales-González JM, Rodrigo-Comino J. Geomorfología de los montes de Málaga: pasado, presente y ¿futuro? BAETICA Estudios De Hist Mod Y Contemp. 2015(33):81–109. (In Spanish). doi:10.24310/baetica.2011.v0i33.111. [Google Scholar] [CrossRef]

33. González-Moreno MT, Rodrigo-Comino J. Geostatistical vegetation filtering for rapid UAV-RGB mapping of sudden geomorphological events in the Mediterranean areas. Drones. 2025;9(6):441. doi:10.3390/drones9060441. [Google Scholar] [CrossRef]

34. Rodrigo-Comino J, Al-Shammary AAG, Durán-Zuazo VH, Serrano-Bernardo F, Caballero-Calvo A, Rodríguez-Galiano V. The limits of RGB-based vegetation indexes under canopy degradation: insights from UAV monitoring of harvested cereal fields. Drones Auton Veh. 2026;3(1):10021. doi:10.70322/dav.2025.10021. [Google Scholar] [CrossRef]

35. Beven KJ, Kirkby MJ. A physically based, variable contributing area model of basin hydrology/Un modèle à base physique de zone d’appel variable de l’hydrologie du bassin versant. Hydrol Sci Bull. 1979;24(1):43–69. doi:10.1080/02626667909491834. [Google Scholar] [CrossRef]

36. Moore ID, Grayson RB, Ladson AR. Digital terrain modelling: a review of hydrological, geomorphological, and biological applications. Hydrol Process. 1991;5(1):3–30. doi:10.1002/hyp.3360050103. [Google Scholar] [CrossRef]

37. Auerswald K. Influence of initial moisture and time since tillage on surface structure breakdown and erosion of a loessial soil. Catena Suppl. 1993;24:93–101. [Google Scholar]

38. Bronick CJ, Lal R. Soil structure and management: a review. Geoderma. 2005;124(1–2):3–22. doi:10.1016/j.geoderma.2004.03.005. [Google Scholar] [CrossRef]

39. Kyratzis AC, Skarlatos DP, Menexes GC, Vamvakousis VF, Katsiotis A. Assessment of vegetation indices derived by UAV imagery for durum wheat phenotyping under a water limited and heat stressed Mediterranean environment. Front Plant Sci. 2017;8:1114. doi:10.3389/fpls.2017.01114. [Google Scholar] [PubMed] [CrossRef]

40. Vanderlinden K, Martínez G, Ramos M, Mateos L. Relevance of NDVI, soil apparent electrical conductivity and topography for variable rate irrigation zoning in an olive grove. Precis Agric. 2024;25(6):3086–108. doi:10.1007/s11119-024-10191-4. [Google Scholar] [CrossRef]

41. Chen HW, Chen CY, Yang PZ. Using drone-based multispectral imaging for investigating gravelly debris flows and geomorphic characteristics. Environ Earth Sci. 2024;83(8):247. doi:10.1007/s12665-024-11544-y. [Google Scholar] [CrossRef]

42. Wilson JP, Gallant JC. Terrain analysis: principles and applications. Hoboken, NJ, USA: John Wiley & Sons, Inc.; 2000. [Google Scholar]

43. Langhammer J. UAV monitoring of stream restorations. Hydrology. 2019;6(2):29. doi:10.3390/hydrology6020029. [Google Scholar] [CrossRef]

44. Demir S, Başyiğit L. Evaluating bare soil properties and vegetation indices for digital farming applications from UAV-based multispectral images. DataSCI. 2023;6(4):5–10. doi:10.3390/rs70404026. [Google Scholar] [CrossRef]

45. Sørensen R, Zinko U, Seibert J. On the calculation of the topographic wetness index: evaluation of different methods based on field observations. Hydrol Earth Syst Sci. 2006;10(1):101–12. doi:10.5194/hess-10-101-2006. [Google Scholar] [CrossRef]

46. Gewali UB, Monteiro ST, Saber E. Machine learning based hyperspectral image analysis: a survey. arXiv:1802.08701. 2018. [Google Scholar]

47. Alameda D, Villar R. Moderate soil compaction: implications on growth and architecture in seedlings of 17 woody plant species. Soil Tillage Res. 2009;103(2):325–31. doi:10.1016/j.still.2008.10.029. [Google Scholar] [CrossRef]

48. López-Vicente M, Álvarez S. Stability and patterns of topsoil water content in rainfed vineyards, olive groves, and cereal fields under different soil and tillage conditions. Agric Water Manag. 2018;201(3):167–76. doi:10.1016/j.agwat.2018.02.004. [Google Scholar] [CrossRef]

49. Basic F, Kisic I, Mesic M, Nestroy O, Butorac A. Tillage and crop management effects on soil erosion in central Croatia. Soil Tillage Res. 2004;78(2):197–206. doi:10.1016/j.still.2004.02.007. [Google Scholar] [CrossRef]

50. Boulal H, Gómez-Macpherson H, Gómez JA, Mateos L. Effect of soil management and traffic on soil erosion in irrigated annual crops. Soil Tillage Res. 2011;115(5):62–70. doi:10.1016/j.still.2011.07.003. [Google Scholar] [CrossRef]