Remote Sensing for Within-Plot Soil Variability Assessment Using NDVI Dispersion Metrics

Abstract

While within-plot soil variability strongly influences crop development and agronomic management efficiency, practical and transferable methods for characterizing this variability remain limited. This study proposes a remote sensing framework to assess soil-driven crop heterogeneity using Normalized Difference Vegetation Index (NVDI)-derived parameters and phenological analysis. Crop dynamics were monitored using Sentinel-2 imagery to identify key phenological stages (start of season, maximum development, and senescence), while high-resolution PlanetScope imagery was used to analyze within-plot spatial variability. NDVI dispersion metrics were used as proxies for spatial variability in crop development across soil units and under homogeneous management conditions to isolate soil effects. Results showed that NDVI standard deviation effectively captured spatial heterogeneity in crop development, revealing greater variability in Glacis Slope soils compared with the more homogeneous Platform unit. Crop senescence emerged as the most informative stage for detecting soil-driven variability, while emergence patterns were influenced by sowing dates and early establishment conditions. Analyses under homogeneous management confirmed that soil properties controlling soil moisture dynamics govern spatial variability in crop response. The proposed methodology provides a robust, low-cost, and transferable approach for identifying within-plot variability and supports precision agriculture by enabling site-specific management decisions based on satellite-derived indicators.

Introduction

Within-plot soil variability is a well-recognized factor influencing crop development, yield stability, and agronomic management efficiency (Bouma and Finke, 1993; Figueroa Jáuregui et al., 2018). Spatial differences in soil properties such as texture, structure, organic matter content, salinity, and water availability often lead to heterogeneous crop responses within the same field, even under uniform management practices (Auerswald et al., 1997). This heterogeneity represents a non-negligible challenge for agricultural production systems, as it can reduce resource-use efficiency and limit the effectiveness of uniform management strategies. Consequently, the identification and characterization of within-plot soil variability have become central objectives in precision agriculture, particularly to support variable rate applications (VRA) (Schimmelpfennig and Lowenberg-DeBoer, 2021).

Remote sensing has consolidated its role as an operational tool to address this challenge by providing spatially continuous and temporally repeated information on crop development (Khorram et al., 2012). Satellite-based observations enable the detection of spatial patterns in crop performance that would otherwise require extensive and costly field surveys. In this context, the availability of open-access satellite imagery has significantly increased the feasibility of low-cost, scalable approaches for monitoring agricultural systems. Such approaches are especially relevant in irrigated cereal systems, where relatively small variations in soil properties can translate into marked differences in crop development within individual plots.

Alternative remote sensing approaches, including synthetic aperture radar (SAR), LiDAR data, and proximal or handheld sensors, can provide complementary information on soil and vegetation properties and have demonstrated their potential in agricultural applications. However, these techniques often involve higher data-processing complexity, specialized expertise, or additional operational costs, which may constrain their routine implementation in farm-scale management. In contrast, optical satellite imagery offers an operational balance between information content, spatial coverage, and ease of implementation, making it particularly suitable for scalable assessments.

Vegetation indices derived from optical satellite imagery have been widely applied to characterize crop condition and spatial variability within agricultural fields (Domínguez Carrera, 2020). These indices are obtained by combining and comparing the spectral bands detected by the sensor, which range from the typical four bands to several tens or even hundreds in multispectral or hyperspectral imagery (Roman-Gonzalez and Vargas-Cuentas, 2013). This combination describes the spectral response of vegetation and is commonly used to monitor crop growth, assess biomass development, and identify zones of reduced vigor. Although vegetation indices primarily reflect plant characteristics, crop growth and spatial patterns are strongly influenced by soil conditions; therefore, spatial variability in vegetation indices can indirectly reflect underlying soil heterogeneity through crop–soil interactions. This indirect relationship makes vegetation indices a practical proxy for assessing within-plot variability in many agronomic contexts.

Among the available vegetation indices, the Normalized Difference Vegetation Index (NDVI), originally introduced by (Rouse et al., 1974), remains one of the most widely used metrics in agricultural remote sensing. Its widespread adoption is largely explained by its simple formulation, intuitive interpretation, and broad availability from both open-access satellite missions, such as Sentinel-2, and high-resolution commercial constellations, such as PlanetScope. NDVI has been consistently applied across a wide range of spatial scales and cropping systems, providing a robust and comparable baseline for monitoring crop dynamics (Amezketa Lizarraga et al., 2011; Girimonte and García Fronti, 2020; Stoy et al., 2022). While it is well established that NDVI exhibits limitations, particularly saturation under dense canopies at peak biomass (Xue and Su, 2017), its consistency, temporal continuity, and operational simplicity make it well suited for methodological frameworks focused on practical applicability. Importantly, this study does not rely on NDVI to discriminate subtle differences in high-biomass conditions but rather exploits its spatial dispersion and temporal evolution to characterize within-plot heterogeneity.

Within-plot variability fluctuates throughout the crop cycle as soil–crop interactions evolve. Several studies distinguish three key phenological stages to interpret crop response and its environmental controls: the start of the season (SOS), corresponding to crop emergence and early establishment; the period of maximum development (MAX), associated with peak canopy cover; and the end of the season (EOS), marking crop senescence (Ren et al., 2017). These phenological stages represent distinct phases of crop–soil interaction. Early growth is strongly influenced by soil physical properties affecting emergence and establishment, whereas later stages reflect differences in water availability and soil water-holding capacity that become more evident as crops approach senescence (Stadler et al., 2015). Standardizing analysis around these phenological stages could facilitate consistent comparisons among crops, years, and soil units.

Despite the extensive use of NDVI in agricultural studies, many existing approaches for characterizing within-field variability rely on complex processing workflows, multiple data sources, or intensive field measurements, which can limit their reproducibility and operational transferability (Gopp et al., 2017; Whetton et al., 2017). In particular, relatively few studies focus explicitly on simple, robust metrics of within-plot heterogeneity and on the identification of crop stages that maximize the detectability of soil-driven variability (Verhulst et al., 2009). There is therefore a need for methodologies that prioritize simplicity, comparability across seasons, and applicability using freely available satellite data, while still providing meaningful insights into soil–crop interactions at the plot scale.

Freely available satellite imagery, such as Sentinel-2 (EOS Data Analytics, 2026), provides continuous monitoring of crop dynamics at the plot scale and has greatly expanded opportunities for low-cost agricultural assessments. When finer spatial detail is required, high-resolution imagery from commercial constellations, such as PlanetScope (Planet Labs PBC, 2025), offers enhanced capability for detecting within-plot variability and is increasingly accessible for research and applied purposes. Together, these complementary data sources enable multi-scale analysis of crop development and soil–crop interactions.

Under the hypothesis that spatial variability of NDVI may reflect underlying soil heterogeneity, dispersion metrics such as the standard deviation can reveal uneven crop development patterns. The temporal evolution of NDVI across phenological stages further influences the detection and interpretation of soil-driven variability.

Building on this hypothesis, this study aims to develop a practical and transferable framework for characterizing within-plot soil variability using a two-stage remote sensing approach based on satellite-derived NDVI. By integrating NDVI spatial heterogeneity, expressed through dispersion metrics, with its temporal dynamics at key phenological stages (SOS, MAX, and EOS), the framework aims to identify the crop development stages most informative for detecting soil variability. Applied to irrigated cereal systems, it also evaluates the potential of NDVI-based indicators to support precision agriculture management.

Materials and Methods

The overall methodological framework, summarized in Figure 1, integrates plot-scale variability analysis, controlled within-plot assessment, and multivariate evaluation to identify the soil and environmental factors controlling crop response.

FIGURE 1

Description of Study Area

The study area comprises the 2,324 ha Sardeta ravine basin, located in the municipalities of Montesusín and Robres (Huesca province, Aragón, Spain), within the Los Monegros region, approximately 28 km south of Huesca, with an elevation ranging between 300 and 350 m above sea level. This location within Spain is shown in Figure 2. The area is characterized by an arid climate typical of the Ebro Valley and Lower Aragón, with hot summers, cold winters, and scarce rainfall throughout most of the year. In Montesusín, land consolidation and irrigation modernization have been carried out (mostly from surface irrigation to sprinkler irrigation), while in Robres this has not yet been done, resulting in the predominance of small, terraced plots with surface irrigation.

FIGURE 2

Crop distribution was characterized using official agricultural statistics from the SINGEAR database (Gobierno de Aragón, 2026), with mean values for 2019–2023. In both municipalities, winter cereals are dominated by barley, which represents 33.23% and 47.17% of the cultivated area in Robres and Montesusín, respectively, compared with 18.01% and 11.47% for wheat. In irrigated areas of Montesusín, alfalfa (25.61% of the total area) and maize (17.58%) are also widely cultivated. Maize is frequently grown as a second crop following barley, representing 48.34% of the barley area, whereas double cropping after wheat is considerably less common (16.02% of the wheat area) and in some years accounts for less than 3% of the cultivated surface. This difference is mainly associated with the longer growth cycle of wheat, which delays the establishment of the second crop and reduces its agronomic viability. These patterns indicate that barley and maize are the predominant cropping systems in the study area and were therefore selected as representative crops for the analysis.

The soils of the basin developed on limestone sediments and are classified as Xerorthent (poorly developed) and Calcixerept (moderately developed). A soil map of the basin is shown in Figure 2, identifying seven distinct homogeneous edaphic units (Usón Murillo et al., 2023). Table 1 describes these units. Due to the complexity and low representativeness of units 2 and 3 (Shallow Soils and Saline Association), these units were excluded from the variability analysis. Additionally, units 4 and 6 (Valley Bottom and Floodplain) were grouped into a single unit belonging to the same soil subgroup.

Agronomic practices, including sowing dates, irrigation management, and fertilization, were uniform within each plot; therefore, the observed spatial variability in crop development can be primarily attributed to differences in soil properties.

The applied methodological framework was designed to characterize soil-driven crop variability under irrigated semi-arid conditions through a progressive analysis from plot-scale assessment to controlled within-plot evaluation.

NDVI-Based Assessment of Plot-Scale Variability Among Soil Units

Among available remote sensing approaches, the Normalized Difference Vegetation Index (NDVI) was selected due to its simplicity, robustness, and wide availability from optical satellite imagery. NDVI provides a consistent and operational indicator of crop development and spatial variability while requiring minimal data processing, making it suitable for cost-effective and transferable agricultural applications. It has been widely used as an indirect indicator of soil-driven crop variability, showing strong relationships with crop productivity, soil physicochemical properties, and spatial yield patterns (Amezketa Lizarraga et al., 2011; Gopp et al., 2017; Gopp and Savenkov, 2019; Lobell et al., 2003). Although NDVI may exhibit saturation under dense canopy conditions, particularly during peak crop development, this limitation does not compromise its suitability for the objectives of this study. The analysis focuses on spatial variability and temporal dynamics rather than absolute biomass estimation and includes phenological stages where saturation effects are less critical. Moreover, the use of dispersion metrics (standard deviation) enhances the detection of within-plot heterogeneity even under high canopy cover conditions.

Given the objective of developing a practical and transferable methodology for within-plot variability assessment, NDVI derived from optical satellite imagery was selected as the most efficient and cost-effective indicator. Its wide availability, minimal processing requirements, and low implementation cost make it particularly suitable for research and operational applications under limited resource conditions.

NDVI time series were derived from Sentinel-2 satellite Level-2A surface reflectance imagery using Sentinel Hub EO Browser (Sentinel Hub, 2026). Sentinel-2 provides multispectral observations with a spatial resolution of 10 × 10 m per pixel, enabling the characterization of crop development at the plot scale. For each study plot, polygon boundaries previously defined in QGIS (version 3.18) were uploaded to EO Browser to extract NDVI values for all available satellite acquisitions during the 2019–2023 period. NDVI was computed from the red and near-infrared spectral bands and exported as pixel-level data in comma-separated values (CSV) format.

The analysis was conducted using agricultural plots for which detailed soil characterization was available, allowing the direct linkage between crop development and soil properties. A total of 12 plots were initially selected based on soil profile information; in four plots complete soil profile pits were excavated and sampled by horizons, while the remaining plots were characterized through auger sampling. For the satellite-based analysis, only plots located within the most representative soil units of the study basin were considered. Two edaphic units characterized by low spatial representativeness were excluded, resulting in a final dataset of 11 plots distributed across Platforms (three plots), Glacis–Slope (four plots), Redoxmorphic Floodplain (one plot), and Valley Bottom–Floodplain (three plots).

To ensure the reliability of the NDVI time series, cloud-contaminated observations were removed using the cloud detection algorithm implemented in EO Browser, which is based on the Braaten–Cohen–Yang method (Sentinel Hub, 2026). Only satellite acquisitions with 0% cloud coverage within each plot were retained to minimize potential contamination of NDVI values. Although this strict selection may reduce the temporal density of the image series and potentially lead to some degree of under-sampling in the NDVI time series, the remaining observations were sufficient to identify the main phenological stages analyzed in this study.

Pixel-level NDVI data were imported into QGIS (version 3.18), where zonal statistics were computed for each plot and acquisition date. Mean, median, standard deviation, and the 10th and 90th percentiles were calculated from all pixels within each plot, transforming spatially distributed information into parcel-level metrics describing crop development and within-plot variability. These statistics were associated with the acquisition date of each satellite image, enabling the construction of parcel-level NDVI time series for analyzing crop temporal dynamics. NDVI dispersion metrics, primarily represented by the standard deviation, were used to characterize within-plot spatial heterogeneity. These indicators enhance the detection of soil-driven variability and remain sensitive to spatial differences in crop performance even when NDVI values approach saturation.

The resulting NDVI time series were analyzed to characterize crop dynamics and identify crop cycles throughout the study period. The temporal evolution of NDVI values enabled the differentiation between perennial and annual crops based on their distinct growth patterns. Perennial crops such as alfalfa exhibited multiple growth cycles within a year, resulting in recurrent NDVI peaks that hindered the consistent identification of phenological stages (Figure 3A). In contrast, annual cereal crops showed a single, well-defined growth cycle, typically following the sequence barley–second maize crop (Figure 3B). Consequently, the analysis focused on winter cereals (mainly barley) and summer crops (mainly maize), while perennial crops were excluded from further analysis.

FIGURE 3

Crop phenological stages were determined from NDVI time series using a threshold-based approach combined with a temporal consistency criterion to identify the start of the season (SOS), maximum development (MAX), and end of the season (EOS). Although previous studies define phenological transitions based on relative thresholds of seasonal NDVI amplitude (Ren et al., 2017; Yin et al., 2025), these approaches do not fully account for signal variability associated with intensive irrigation management, heterogeneous crop emergence, or transient disturbances commonly observed in semi-arid irrigated systems. Therefore, threshold values and detection criteria were adapted to local agronomic conditions to ensure robust identification of crop development stages.

In the study area, crops are intensively managed under irrigation, typically reaching high NDVI values and exhibiting rapid senescence following the cessation of irrigation prior to harvest. Moreover, NDVI temporal series may present short-term fluctuations caused by soil background effects, weed growth during early stages, spatial heterogeneity within plots, or temporary irrigation interruptions. These factors may produce transient NDVI increases or decreases that do not represent actual phenological transitions.

Early vegetation signals detected before crop establishment may represent weed growth or residual vegetation rather than the primary crop. This can confound the identification of phenological transitions when using vegetation index time series (Sarvia et al., 2026). However, crop establishment in the study area generally occurs rapidly due to warm conditions and intensive agronomic management, particularly under double-cropping systems, and the widespread use of pre-emergence herbicides typically limits early weed infestation.

To avoid misclassification of phenological transitions, a temporal consistency criterion based on consecutive observations was defined:

The SOS stage was determined as the fourth consecutive observation following three successive NDVI increases with values exceeding 0.4 (approximately 40% of seasonal amplitude). This threshold was selected to ensure sustained crop establishment while minimizing the influence of early signal variability. The criterion was particularly relevant for maize, where typical row spacing (∼0.7 m) results in substantial soil exposure and irregular NDVI values during initial growth stages.

The MAX stage corresponded to the maximum NDVI value recorded during the crop cycle, representing peak canopy development. Although NDVI may saturate under dense canopy conditions, this stage provides a consistent reference for crop maturity under irrigated management.

The EOS stage was defined as the fourth consecutive observation following three successive NDVI decreases after the maximum value, when NDVI values fall below 0.9 (approximately 90% of seasonal amplitude), representing the onset of crop senescence. This criterion captures the initial decline in vegetation activity rather than advanced senescence, allowing the detection of spatial differences in crop drying related to soil water availability and water-holding capacity. The use of consecutive observations further ensured that detected transitions reflected sustained phenological changes rather than temporary NDVI fluctuations.

Although this approach enables comparison among soil units at the plot scale, differences in agronomic management between farmers may partially obscure soil-driven crop responses. This limitation is addressed through a controlled within-plot analysis presented in the following section.

Statistical analyses were performed using IBM SPSS Statistics (version 26). Analysis of variance was used to evaluate differences in NDVI-derived metrics among soil units, and mean comparisons were conducted using the Tukey B test when significant differences were detected (p < 0.05).

Comparison of Soil Units Under Homogeneous Management

To minimize the influence of agronomic management and isolate soil-driven crop responses, we conducted a controlled within-plot analysis. This approach allows the isolation of soil effects under homogeneous management conditions.

A large agricultural plot (21 ha) located in the municipality of Montesusín (41.887777°N, 0.386029°W) was selected because it includes two representative soil units (Platform and Glacis–Slope) with contrasting edaphic characteristics. The field is managed uniformly with respect to sowing date, fertilization, irrigation, phytosanitary treatments, and harvesting operations. Under these homogeneous management conditions, differences in crop development can be primarily attributed to soil variability.

To characterize crop performance within each soil unit while avoiding non-representative areas such as field edges, headlands, or traffic zones, six polygons of approximately 1,000 m² were randomly delineated within each unit, resulting in twelve sampling areas. This random selection ensured that the sampled zones represented dominant soil conditions while minimizing localized disturbances. Figure 4 illustrates the spatial distribution of the sampling polygons within the study plot, with striped areas corresponding to the Glacis–Slope unit and dotted areas to the Platform unit.

FIGURE 4

Crop cycles within the selected plot were first identified using Sentinel-2 NDVI time series derived from EO Browser, following the procedure described in Section NDVI-Based Assessment of Plot-Scale Variability Among Soil Units. Phenological stages (SOS, MAX, and EOS) were determined from these time series. Subsequently, high-resolution PlanetScope imagery (3 m spatial resolution) corresponding to the identified dates was acquired, allowing these stages to be analyzed at a finer spatial resolution. This two-step approach enabled precise temporal identification of crop development while improving the spatial characterization of within-plot variability.

NDVI values were extracted for all crop cycles recorded between 2019 and 2023, including winter crops (barley) and summer crops (maize). The use of a consistent phenological framework enabled comparison of crop responses between soil units under identical management conditions across multiple growing seasons.

Statistical analyses were performed using IBM SPSS Statistics (version 26). Data distribution was evaluated prior to analysis, and a two-factor univariate analysis of variance was conducted to assess the effects of soil unit and crop type on NDVI median values and NDVI standard deviation at the three phenological stages.

NDVI-Driven Multivariate Analysis of Soil-Crop Relationships

To further interpret the soil-driven variability observed under homogeneous management conditions, an analysis was conducted to identify the soil properties most strongly associated with crop performance, using NDVI-derived parameters as indicators of crop response. This approach supports the interpretation of spatial variability in terms of soil–crop interactions.

The analysis incorporated variables grouped into three categories: soil properties (soil unit, clay and sand content, water table depth, depth of redoximorphic features, surface horizon organic matter, and salinity indicators); plot characteristics (plot area, mean elevation, slope, and irrigation system); and crop-related variables (crop type, growing season, and NDVI-derived indicators).

Spearman’s rank correlation analysis was used to evaluate the relationships among the variables. Although most variables followed a normal distribution, Spearman’s coefficient was selected to assess monotonic relationships, which are expected in soil–crop systems where responses may not follow strictly linear patterns. This approach provides a robust assessment of directional associations while accommodating potential non-linear trends. Correlation coefficients (ρ) were considered significant at the 95% confidence level.

In addition, principal component analysis (PCA) was performed to explore multivariate relationships and identify the main factors explaining variability within the soil–crop system. This approach facilitates the interpretation of complex interactions by reducing data dimensionality and highlighting the variables that contribute most strongly to soil–crop system variability.

Results and Discussion

Plot-Scale Variability Among Soil Units

To evaluate crop response across soil units at the plot scale, NDVI-derived metrics were analyzed at key phenological stages (SOS, MAX, and EOS). In addition, the seasonal MEAN NDVI was considered to provide an integrated indicator of crop performance throughout the growing cycle. NDVI median values were used as an indicator of overall crop development, whereas the NDVI standard deviation was used to characterize within-plot spatial heterogeneity. A two-factor analysis of variance was conducted to assess the effects of soil unit and crop type on NDVI metrics across growth stages.

NDVI values did not differ significantly among soil units or between crop types across the evaluated growth stages. These results indicate that none of the soil units imposed significant limitations on crop development and that growth patterns were comparable between summer and winter crops. Significant differences were observed only among phenological stages, with the lowest values at emergence, a peak at maximum development, and a decline during senescence (Figure 5).

FIGURE 5

The absence of significant differences in NDVI median values (0.6577 for summer crops and 0.6502 for winter crops) is consistent with NDVI saturation under dense canopy conditions (Tucker, 1979). This occurs because red reflectance saturates earlier than near-infrared, reducing sensitivity at high biomass levels (Salvador-Castillo et al., 2021), a behavior linked to the strong absorption of chlorophyll in the red spectral region (Gitelson et al., 2003). In this context, the limited discrimination at maximum canopy development reinforces the importance of focusing on SOS and EOS stages, as well as on dispersion metrics, to better capture soil-driven variability.

In contrast, NDVI standard deviation differed significantly among soil units (Figure 6). The Glacis Slope unit exhibited the highest variability, whereas the Platform unit showed the lowest values, with Floodplain and Redoximorphic Floodplain displaying intermediate behavior (letter groupings in Figure 6). A tendency toward higher variability in winter crops compared with summer crops was also observed (p = 0.071), suggesting greater spatial heterogeneity under less uniform water supply conditions. This pattern may reflect more homogeneous irrigation management in maize, whereas barley shows greater dependence on rainfall, leading to increased spatial variability in crop response; alternatively, it could also be influenced by earlier or stronger NDVI saturation under the denser maize canopy, potentially reducing the sensitivity of the index to spatial differences.

FIGURE 6

The crop developmental stage also influenced spatial variability. NDVI standard deviation increased from SOS (0.043) to EOS (0.061), with the highest variability observed during senescence. This suggests that crop senescence is the most suitable stage for detecting soil-driven differences in crop response. In contrast, variability at SOS was more strongly affected by differences in sowing dates and early establishment conditions.

Variability in crop emergence timing among plots was assessed by comparing the SOS dates across all evaluated crops. Substantial differences were observed, with emergence varying by up to 3 months for winter crops and by approximately 2 months for summer crops. These differences reflect variations in sowing dates and management practices among plots and may introduce uncertainty when interpreting soil-driven crop responses. To minimize these effects, subsequent analyses were conducted under homogeneous management conditions.

Within-Plot Variability Under Homogeneous Management

To isolate the influence of soil properties from management effects, a plot with homogeneous agronomic management was analyzed in detail. Uniform irrigation, fertilization, and sowing practices across the field enabled evaluation of soil-driven variability under consistent conditions.

The cropping sequence within the selected plot over the last 5 years was analyzed to provide temporal context for the study. Figure 7 shows the growth cycles of representative winter and summer crops, including barley (2021 and 2023) and maize (2019–2023), together with the timing of the SOS, MAX, and EOS stages. Main-crop maize and double-crop maize following barley are distinguished, with the latter exhibiting a delayed growth cycle due to later planting. This temporal framework confirms the consistency of phenological development across seasons and supports the use of these stages to evaluate soil-driven variability.

FIGURE 7

Crop response is discussed across the key phenological stages.

SOS (Emergence)

NDVI median differed significantly between soil units and crops (Figure 8A), with higher values in the Platform unit and in winter crops. These differences may reflect higher winter soil moisture, which increases reflectance (Remer et al., 2001), and the wider spacing of maize, which reduces early vegetation cover and NDVI values (Espinosa-Espinosa et al., 2017). The coarser texture of the Platform unit may also favor faster establishment.

FIGURE 8

Summer crops exhibited a higher NDVI standard deviation (Figure 8B), indicating greater heterogeneity during periods of low canopy development, whereas the Platform unit showed more homogeneous early growth.

MAX (Maturity)

NDVI median values differed significantly between soil units and crops (Figure 9A), with higher values observed in the Platform unit and summer crops. This pattern indicates more uniform canopy development on the Platform and greater canopy density under summer conditions.

FIGURE 9

NDVI standard deviation differed significantly between soil units (Figure 9B), with greater variability on the Glacis Slope, reflecting its higher soil heterogeneity, whereas the Platform exhibited more uniform crop development. Lower spatial variability was associated with higher NDVI values, reinforcing patterns observed at the plot scale.

EOS (Senescence)

NDVI median values differed significantly only between crops (Figure 10A), with higher values observed in summer crops, reflecting slower canopy drying compared with winter crops. Higher moisture conditions during late-season senescence may also contribute to maintaining NDVI values and reducing soil-related differences (Medida et al., 2023; Ren et al., 2017).

FIGURE 10

NDVI standard deviation showed significant effects of soil unit, crop type, and their interaction (Figure 10B), indicating that senescence is the most suitable stage for detecting crop variability. Spatial variability was greater on the Glacis Slope than on the Platform and was more pronounced in winter crops. These patterns reflect differences in soil water-holding capacity that become more evident after irrigation ceases and evapotranspiration demand increases, thereby enhancing variability in the more heterogeneous unit.

Overall, under homogeneous management, NDVI median proved to be a reliable indicator of crop growth, while NDVI standard deviation captured spatial variability in crop development. Eliminating management effects revealed clearer differences between crops and soil units, with soil influence becoming most evident during senescence, confirming this stage as optimal for detecting soil-driven variability.

Multivariate Analysis of Soil-Crop Relationships

To better understand the mechanisms underlying the observed variability, relationships among soil properties and crop response were examined using correlation analysis and principal component analysis (PCA).

Correlation analysis using Spearman’s rank coefficients was conducted to explore relationships among soil properties, plot characteristics, and crop response indicators derived from NDVI metrics. Only statistically significant relationships are presented to facilitate interpretation of the main drivers of crop variability, as summarized in Figure 11. Relationships are described below according to three variable groups.

FIGURE 11

Soil-related variables showed strong interconnections reflecting their influence on soil hydrological behavior and crop development. Soil unit was associated with several edaphic and groundwater-related attributes, confirming its role as an integrative indicator of soil conditions across the watershed. Texture-related properties exhibited contrasting patterns consistent with their control over soil water dynamics: sand content was associated with conditions favoring permeability and drainage, whereas clay content reflected greater water retention and a tendency toward salt accumulation. Organic matter content emerged as a stabilizing factor, being associated with improved soil structure and more homogeneous crop development. Relationships involving salinity and redoximorphic features further indicate the influence of groundwater conditions and soil moisture regimes.

Topographic variables related to plot morphology were closely linked to soil physical properties and hydrological behavior. Slope was associated with increased crop variability, suggesting uneven soil moisture redistribution along the terrain. Elevation showed relationships consistent with differences in drainage conditions and soil development across the landscape. The irrigation system was associated with variability patterns, reflecting differences in water application uniformity and management intensity.

Crop-related indicators derived from NDVI confirmed the sensitivity of remote sensing metrics to soil-driven heterogeneity. NDVI standard deviation, used as an indicator of within-plot crop variability, showed strong monotonic associations with terrain and soil moisture–related variables. In contrast, NDVI mean values were inversely related to variability, indicating that greater crop vigor is associated with more homogeneous canopy development.

Collectively, the correlation patterns indicate that spatial variability in crop development is primarily controlled by soil-water dynamics. The interaction of texture, topography, groundwater influence, salinity, and organic matter governs soil moisture distribution, which in turn determines crop response. These relationships support the use of NDVI standard deviation as a robust indicator of soil-driven heterogeneity in agricultural landscapes.

Principal Component Analysis (PCA) was performed to visualize how variables group together and to provide a clearer interpretation of the relationships identified in the correlation analysis. This multivariate approach complements the correlation results by revealing underlying gradients that control crop variability.

The first two principal components explained 37.33% of the total variance, reflecting the complexity of the system and the diversity of variables considered, including soil properties, terrain attributes, and crop indicators collected across different crops and years.

Component 1 was primarily associated with textural fractions and elevation, with clay and sand showing opposite contributions. This axis represents a gradient of soil texture and permeability that influences water movement and drainage conditions.

Component 2 was strongly related to soil unit, slope, and organic matter content, representing a gradient of soil heterogeneity and moisture retention capacity across the landscape.

The loading plot (Figure 12) reveals meaningful groupings of variables that support the interpretation of soil hydrological processes. Variables located close to each other are positively correlated, whereas those positioned on opposite sides of the origin are negatively correlated. Sand content clustered with water table depth, reflecting high permeability and enhanced vertical water movement. The irrigation system was grouped with salinity, indicating that areas affected by soil salinity tend to maintain traditional surface irrigation practices. Elevation was clustered with redoximorphic feature depth and organic matter content, consistent with the coarser texture and improved drainage conditions of the Platform unit. The soil unit was grouped with NDVI standard deviation, reinforcing the role of this indicator as a proxy for soil-driven crop variability. To facilitate interpretation, the loadings of the variables on the first two principal components are presented in Table 2.

FIGURE 12

Overall, the PCA highlights the dominant influence of soil texture, terrain, and groundwater on soil moisture dynamics, revealing patterns consistent with the associations identified in the correlation analysis. These processes govern crop spatial variability, particularly in finer-textured areas where salinity accumulation and redoximorphic conditions are more likely to occur.

Taken together, these results lead to the following conclusions.

Conclusion

NDVI-derived metrics proved effective for assessing crop development and spatial variability at the plot scale. NDVI mean values reflected overall crop performance, whereas NDVI standard deviation captured within-plot variability associated with soil heterogeneity.

The use of key phenological stages (SOS, MAX, and EOS) simplified the analysis and facilitated interpretation of crop response. Variability in agronomic management partially masked soil-related differences at the watershed scale, but these became clearer under homogeneous management conditions.

NDVI standard deviation revealed differences in crop development between soil units, with more homogeneous growth in the Platform unit and greater heterogeneity in the Glacis Slope unit. Crop senescence (EOS) showed the greatest spatial variability, making it the most suitable stage for detecting soil-driven differences, while emergence timing (SOS) also showed potential for supporting precision management of the current crop.

Taken together, these findings demonstrate that NDVI variability metrics provide a robust approach for identifying soil-driven crop heterogeneity and supporting site-specific management.

The main limitations relate to differences in agronomic management and to the inherent constraints of NDVI under dense canopy conditions. Future research should evaluate additional vegetation indices and sensing approaches, including the integration of synthetic aperture radar and optical observations, to provide complementary information on soil moisture and structural properties linked to spatial variability. This methodology should also be tested across diverse agroecosystems to further support precision agriculture and sustainable land management.

References

• 1

  Amezketa LizarragaE.Urdanoz MeladoV.BarinagarrementeriaI.Albizua HuarteL.BerkaneY.PorteroC.et al (2011). Validación De Índices Espectrales Para Detectar Salinidad Edáfica En Cebada Mediante Sensores Electromagnéticos Terrestres. Teledetección Bosques Y Cambio Climático, 277–280.

  • Google Scholar

• 2

  AuerswaldK.SippelR.KainzM.DemnzelM.ScheiizostA.SinowskiW.et al (1997). The Crop Response to Soil Variability in an Agroecosystem. Adv. GeoEcology30, 39–53.

  • Google Scholar

• 3

  BoumaJ.FinkeP. A. (1993). “Origin and Nature of Soil Resource Variability,” in Proceedings of Soil Specific Crop Management: A Workshop on Research and Development Issues. Editors RobertP. C.RustR. H.LarsonW. E. (John Wiley & Sons, Ltd), 1–13. 10.2134/1993.SOILSPECIFICCROP.C1

  • CrossRef
  • Google Scholar

• 4

  Domínguez CarreraE. M. (2020). Diseño E Implementación De Un Sistema Óptico Para Captación De Imágenes Multiespectrales a Ser Usado En La Detección De Áreas Infectadas Por Hongos En Hojas De Plantas. Quito: ESCUELA POLITÉCNICA NACIONAL.

  • Google Scholar

• 5

  EOS Data Analytics (2026). Sentinel-2 Images. Available online at: https://eos.com/find-satellite/sentinel-2/ (Accessed February 17, 2026).

  • Google Scholar

• 6

  Espinosa-EspinosaJ. L.Espinosa-EspinosaJ. L.Palacios-VélezE.Tijerina-ChávezL.Flores-MagdalenoH.Quevedo-NolascoA.et al (2017). Sistema De Monitoreo Satelital Para El Seguimiento Y Desarrollo De Cultivos Del Distrito De Riego 038. Tecnol. ciencias del agua8, 95–104. 10.24850/J-TYCA-2017-01-07

  • CrossRef
  • Google Scholar

• 7

  Figueroa JáureguiM. de L.Martínez MenezM. R.Ortiz SolorioC. A.Fernández ReynosoD. (2018). Influence of Factors of Soil Formation in the Properties of Soils in Mixteca, Oaxaca, México. Terra Latinoam.36, 287–299. 10.28940/terra.v36i3.259

  • CrossRef
  • Google Scholar

• 8

  GirimonteP.García FrontiJ. (2020). El Índice NDVI Y La Clasificación De Áreas Sembradas Aprendizaje Automático No Supervisado “k-means.”. Rev. Investig. Model. matemáticos Apl. la Gest. la Econ.I, 1–15.

  • Google Scholar

• 9

  GitelsonA. A.GritzY.MerzlyakM. N. (2003). Relationships Between Leaf Chlorophyll Content and Spectral Reflectance and Algorithms for Non-Destructive Chlorophyll Assessment in Higher Plant Leaves. J. Plant Physiol.160, 271–282. 10.1078/0176-1617-00887

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 10

  Gobierno de Aragón (2026). SINGEAR (Sistema Informático De Gestión Estadística Agroalimentaria De Aragón). Available online at: https://www.aragon.es/estadisticas-agrarias/singear (Accessed 17 February, 2026).

  • Google Scholar

• 11

  GoppN. V.SavenkovO. A. (2019). Relationships Between the NDVI, Yield of Spring Wheat, and Properties of the Plow Horizon of Eluviated Clay-Illuvial Chernozems and Dark Gray Soils. Eurasian Soil Sci.52, 339–347. 10.1134/S1064229319030050/TABLES/3

  • CrossRef
  • Google Scholar

• 12

  GoppN. V.NechaevaT. V.SavenkovO. A.SmirnovaN. V.SmirnovV. V. (2017). Indicative Capacity of NDVI in Predictive Mapping of the Properties of Plow Horizons of Soils on Slopes in the South of Western Siberia. Eurasian Soil Sci.50, 1332–1343. 10.1134/S1064229317110060/METRICS

  • CrossRef
  • Google Scholar

• 13

  KhorramS.NelsonS. A. C.KochF. H.van der WieleC. F. (2012). “Remote Sensing,” in Springerbriefs in Space Development. 1st ed. New York: Springer. 10.1007/978-1-4614-3103-9

  • CrossRef
  • Google Scholar

• 14

  LobellD. B.AsnerG. P.Ortiz-MonasterioJ. I.BenningT. L. (2003). Remote Sensing of Regional Crop Production in the Yaqui Valley, Mexico: Estimates and Uncertainties. Agric. Ecosyst. Environ.94, 205–220. 10.1016/S0167-8809(02)00021-X

  • CrossRef
  • Google Scholar

• 15

  MedidaS. K.Prasuna RaniP.Suneel KumarG. V.Geetha SireeshaP. V.KranthiK. C.VinushaV.et al (2023). Detection of Water Deficit Conditions in Different Soils by Comparative Analysis of Standard Precipitation Index and Normalized Difference Vegetation Index. Heliyon9, 1–17. 10.1016/J.HELIYON.2023.E15093

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 16

  Planet Labs PBC (2025). Satellite Imagery Analytics. Available online at: https://www.planet.com/products/satellite-imagery-of-earth/ (Accessed 21 March, 2025).

  • Google Scholar

• 17

  RemerL. A.WaldA. E.KaufmanY. J. (2001). Angular and Seasonal Variation of Spectral Surface Reflectance Ratios: Implications for the Remote Sensing of Aerosol over Land. IEEE Trans. Geoscience Remote Sens.39, 275–283. 10.1109/36.905235

  • CrossRef
  • Google Scholar

• 18

  RenJ.CampbellJ. B.ShaoY. (2017). Estimation of SOS and EOS for Midwestern US Corn and Soybean Crops. Remote Sens.9, 1–14. 10.3390/RS9070722

  • CrossRef
  • Google Scholar

• 19

  Roman-GonzalezA.Vargas-CuentasN. I. (2013). Análisis De Imágenes Hiperespectrales. Revista Ingenieria & Desarrollo35, 14–17. Available online at: https://hal.science/hal-00935014v1 (Accessed February 10, 2026).

  • Google Scholar

• 20

  RouseJ. W.Jr.HaasR. H.SchellJ. A.DeeringD. W. (1974). Monitoring Vegetation Systems in the Great Plains with ERTS.Texas: NASA. Goddard Space Flight Center 3d ERTS-1 Symp., Vol. 1, Sect. A.

  • Google Scholar

• 21

  Salvador-CastilloJ. M.Bolaños-GonzálezM. A.Palacios-VélezE.Palacios-SánchezL. A.López-PérezA.Muñoz-PérezJ. M. (2021). Estimation of Fractional Vegetation Cover and Canopy Nitrogen Content in Corn by Remote Sensing. Terra Latinoam.39, 1–11. 10.28940/TERRA.V39I0.899

  • CrossRef
  • Google Scholar

• 22

  SarviaF.De PetrisS.FarboA.ChiesaE.OrusaT.PariziaF.et al (2026). Mapping Post-Crop Biomass for Grazing in Alessandria Province, North West Italy: Moving Towards Sustainable Agriculture Through the Support of Remote Sensing. Agric. Syst.234, 104694. 10.1016/j.agsy.2026.104694

  • CrossRef
  • Google Scholar

• 23

  SchimmelpfennigD.Lowenberg-DeBoerJ. (2021). “Precision Agriculture Adoption, Farm Size and Soil Variability,” in Precision Agriculture’21. Editor StaffordJ. (Wageningen Academic Publishers), 769–776. 10.3920/978-90-8686-916-9_92

  • CrossRef
  • Google Scholar

• 24

  Sentinel Hub (2026). EO Browser. Available online at: https://apps.sentinel-hub.com/eo-browser/?zoom=8&lat=42.03258&lng=12.35413&themeId=DEFAULT-THEME&toTime=2026-02-17T11%3A48%3A45.826Z (Accessed February 17, 2026).

  • Google Scholar

• 25

  StadlerA.RudolphS.KupischM.LangensiepenM.van der KrukJ.EwertF. (2015). Quantifying the Effects of Soil Variability on Crop Growth Using Apparent Soil Electrical Conductivity Measurements. Eur. J. Agron.64, 8–20. 10.1016/J.EJA.2014.12.004

  • CrossRef
  • Google Scholar

• 26

  StoyP. C.KhanA. M.WipfA.SilvermanN.PowellS. L. (2022). The Spatial Variability of NDVI Within a Wheat Field: Information Content and Implications for Yield and Grain Protein Monitoring. PLoS One17, 1–18. 10.1371/JOURNAL.PONE.0265243

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 27

  TuckerC. J. (1979). Red and Photographic Infrared Linear Combinations for Monitoring Vegetation. Remote Sens. Environ.8, 127–150. 10.1016/0034-4257(79)90013-0

  • CrossRef
  • Google Scholar

• 28

  Usón MurilloA.DíezM.SampérizM. (2023). Uso De Teledetección Para La Delimitación De Unidades De Suelo Homogéneas. XXXIII Reunión Nac. Suelos Pamplona-Iruña, 2023 Libro resúmenes. 150-151.

  • Google Scholar

• 29

  VerhulstN.GovaertsB.SayreK. D.DeckersJ.FrançoisI. M.DendoovenL. (2009). Using NDVI and Soil Quality Analysis to Assess Influence of Agronomic Management on Within-Plot Spatial Variability and Factors Limiting Production. Plant Soil317, 41–59. 10.1007/S11104-008-9787-X/FIGURES/6

  • CrossRef
  • Google Scholar

• 30

  WhettonR.ZhaoY.ShaddadS.MouazenA. M. (2017). Nonlinear Parametric Modelling to Study How Soil Properties Affect Crop Yields and NDVI. Comput. Electron. Agric.138, 127–136. 10.1016/J.COMPAG.2017.04.016

  • CrossRef
  • Google Scholar

• 31

  XueJ.SuB. (2017). Significant Remote Sensing Vegetation Indices: A Review of Developments and Applications. J. Sens.2017, 1–17. 10.1155/2017/1353691

  • CrossRef
  • Google Scholar

• 32

  YinP.LiX.HeiskanenJ.PellikkaP. (2025). Comparative Analysis of Global Urban Land Surface Phenology Between the MODIS and VIIRS Products and Extraction Methods. J. Environ. Manage.375, 1–14. 10.1016/J.JENVMAN.2025.124326

  • Pubmed Abstract
  • CrossRef
  • Google Scholar