Combining electromagnetic induction and satellite-based NDVI data for improved determination of management zones for sustainable crop production

Accurate delineation of management zones is essential for optimizing resource use and improving yield in precision agriculture. Electromagnetic induction (EMI) provides a rapid, non-invasive method to map soil variability, while the Normalized Difference Vegetation Index (NDVI) obtained with remote sensing captures aboveground crop dynamics. Integrating these datasets may enhance management zone delineation but presents challenges in data harmonization and analysis. This study presents a workflow combining unsupervised classification (clustering) and statistical validation to delineate management zones using EMI and NDVI data in a single 70 ha field of the patchCROP experiment in Tempelberg, Germany. Three datasets were investigated: (1) EMI maps, (2) NDVI maps, and (3) a combined EMI–NDVI dataset. Historical yield data and soil samples were used to refine the clusters through statistical analysis. The results demonstrate that four EMI-based zones effectively captured subsurface soil heterogeneity, while three NDVI-based zones better represented yield variability. A combination of EMI and NDVI data resulted in three zones that provided a balanced representation of both subsurface and aboveground variability. The final EMI–NDVI-derived map demonstrates the potential of integrating multi-source datasets for field management. It provides actionable insights for precision agriculture, including optimized fertilization, irrigation, and targeted interventions, while also serving as a valuable resource for environmental modeling and soil surveying.

1 Introduction

Reliable and accurate agricultural management zones that capture within-field variability affecting crop development can play a pivotal role in sustainable agriculture. Management zones can be used in the context of precision agriculture to optimize farming practices, increase yields, and reduce the use of available resources (Gebbers and Adamchuk, 2010; Janrao et al., 2019). This is not only valuable for profit maximization (Adhikari et al., 2022), but is also vital to meet future climate change and food security challenges (Antle et al., 2017; Chartzoulakis and Bertaki, 2015; Bongiovanni and Lowenberg-Deboer, 2004), such as Goal 2 (Zero Hunger) and Goal 15 (Life on Land) of the United Nations Sustainable Development Goals (SDGs) (Hou et al., 2020; UN, 2021). Generally, management zones aim to consider the impact of various factors that can influence crop productivity and result in yield gaps, a key one being soil heterogeneity and health (Licker et al., 2010). Soil systems can be relatively static in time (Arshad et al., 1997) and are fundamental due to their multifunctional role and impact on ecosystem services (Hamidov et al., 2018). Within these systems, soil properties such as texture, organic matter content, cation exchange capacity, and bulk density greatly influence soil moisture dynamics, salinity, nutrient availability, and other variables affecting crop yield (Kibblewhite et al., 2008; Dobarco et al., 2021) and are thus a good target for management zone delineation. However, soil heterogeneity is not solely responsible for yield losses, and effective management zones should also incorporate other influencing factors to provide a comprehensive and holistic management solution.

Traditional methods for soil characterization to support management zone delineation (Brogi et al., 2021; Geologischer Dienst NRW, 2025) generally rely on laborious in situ sampling and laboratory analysis, which may fail in capturing soil variability with sufficient detail (Kuang et al., 2012). In recent years, advances in proximal soil sensing, defined as methods that utilize sensors positioned near or in direct contact with the soil (Adamchuk et al., 2017), have provided valid alternatives to direct soil sampling (Pradipta et al., 2022). In particular, non-invasive agro-geophysical methods such as electromagnetic induction (EMI) have proven suitable for management zone delineation due to their high mobility (Binley et al., 2015; Garré et al., 2021) and the fact that the measured apparent electrical conductivity (ECa) of the soil is related to key soil properties, such as soil salinity, soil water content, texture, compaction, and organic matter content (Corwin and Lesch, 2003; Abdu et al., 2008; Altdorff et al., 2017; Jadoon et al., 2015; Robinet et al., 2018; Zhu et al., 2010; von Hebel et al., 2018). Modern EMI devices are able to efficiently provide soil information for multiple depth ranges thanks to multi-coil instrumentation (Rudolph et al., 2015; von Hebel et al., 2014; Blanchy et al., 2024; Lueck and Ruehlmann, 2013; Corwin and Scudiero, 2019), especially when supported by a moderate amount of ground-truth data (Brogi et al., 2019). However, the use of EMI alone can show limitations in capturing local aspects that have an impact on yield but that are not strongly influenced by soil variability. For instance, pest and weed infestations can drastically reduce crop productivity, and these factors may not correlate directly with soil variability (Becker et al., 2022; López-Granados, 2011). Additionally, climate change impacts, such as altered precipitation patterns and temperature fluctuations, can affect crop health and yield in ways that EMI cannot detect (Pradipta et al., 2022). Finally, it is also important to stress that accurate EMI mapping generally requires optimal conditions like bare soil, favorable weather, and the absence of confounding factors (James et al., 2003).

An alternative to proximal soil sensing for the delineation of management zones is the use of remote sensing approaches, which enables efficient large-scale data acquisition without the need for direct physical access to the investigated area (Weiss et al., 2020). By using sensors mounted on satellites, airplanes, or drones, remote sensing monitors parameters related to crop health and development (Jin et al., 2019; Liaghat and Balasundram, 2010). For example, vegetation indices such as the Normalized Difference Vegetation Index (NDVI) are generally well-established, simple, and effective proxies for crop health (Carfagna and Gallego, 2005; Stamford et al., 2023; Wang et al., 2020; Xue and Su, 2017). High-resolution (< 5 m) data products from satellites are being increasingly used in precision agriculture (Mohammed et al., 2020; Trivedi et al., 2023). Also, remote sensing platforms like PlanetScope, Sentinel-2, and Landsat offer frequent revisit times, thus providing sufficient temporal resolution to track changes in plant health throughout the growing season (Hunt et al., 2019; Skakun et al., 2021). Despite these advantages, remote sensing data are affected by cloud cover or other suboptimal meteorological conditions (Wilhelm et al., 2000) and primarily capture aboveground information on plant health and biomass; they can thus struggle to provide direct information about the interplay between soil conditions and crop development.

Several studies have explored a combination of EMI and remote sensing methods for the delineation of management zones. For example, von Hebel et al. (2021) combined EMI and drone-based NDVI measurements and found that EMI-based management zones offered consistent insights into soil texture and water content, while the added value of NDVI strongly depended on the timing of the drone measurements and thus on the specific crop conditions. In a similar study, Esteves et al. (2022) showed that integration of EMI and NDVI from Sentinel-2 (10 m resolution) effectively provided zones with distinct soil and crop nutrient characteristics. However, they reported a negative relationship between ECa and NDVI due to local magnesium imbalances and vegetation stress. In addition to EMI and remote sensing, historical yield maps can help in identifying yield trends across years and different cultivated crops. For example, Ali et al. (2022) integrated 7 years of yield data with Landsat-based NDVI and soil sampling over a 9 ha field but ultimately could obtain only a limited subdivision of the field into two management zones with a relatively low resolution of 30 m. Overall, previous studies have made important contributions towards integrating EMI and NDVI data for improved management zone delineation (Corwin and Scudiero, 2019; Ciampalini et al., 2015). However, the results can be influenced by data resolution and acquisition timing as well as by local management and soil–plant interactions, with some studies suggesting that EMI alone can offer sufficient insights into soil patterns (Esteves et al., 2022; von Hebel et al., 2021). Nonetheless, the added value of NDVI holds unexplored potential due to the higher spatial and temporal resolution of recent satellite platforms (Breunig et al., 2020; Georgi et al., 2018).

Machine learning clustering algorithms have been widely used to delineate management zones from spatially distributed datasets such as EMI or NDVI (Saifuzzaman et al., 2019; Castrignanò et al., 2018; Chlingaryan et al., 2018; Zhang and Wang, 2023). For example, Wang et al. (2021) used supervised random forest classification for combining EMI data with environmental covariates to predict soil salinity. Similarly, Brogi et al. (2019) employed supervised learning to combine EMI with soil sampling and generate high-resolution soil maps for a 1 km² agricultural area. However, the results of supervised classification approaches may depend on the interpreter and often need expert knowledge as well as extensive ground-truth data for training (Liakos et al., 2018; Usama et al., 2019). K-means and ISODATA clustering are unsupervised methods used to delineate management zones (Bijeesh and Narasimhamurthy, 2020; Ylagan et al., 2022; Tagarakis et al., 2013), but these approaches can be sensitive to initial conditions and struggle to handle nonlinear relationships in datasets (Geng et al., 2020; Li et al., 2018). Thus, more advanced methods such as self-organizing maps (SOMs) have been successfully used to analyze complicated data structures provided by proximal and remote sensing data (Romero-Ruiz et al., 2024; Moshou et al., 2006; Taşdemir et al., 2012). A remaining key challenge of unsupervised methods is the definition of the optimal number of clusters. Widely used approaches such as the elbow and silhouette method (Saputra et al., 2020) often struggle when applied to nonlinearly distributed or spatially complex datasets (Schubert, 2023) and may thus require subjective judgment or expert knowledge (Liang et al., 2012). To address this challenge, the multi-cluster average standard deviation (MCASD) approach that relies on an evaluation of the intra-cluster variability has recently been introduced (O'Leary et al., 2023) and successfully applied to the integration of complex spatial datasets (O'Leary et al., 2024). However, many of these novel approaches have seen limited applications in agricultural contexts (Khan et al., 2021), and the added value of delineating management zones from datasets of different origin remains unaddressed (Koganti et al., 2024).

Within this context, this study expands on previous research by combining high-resolution multi-coil EMI and satellite-based NDVI data within a harmonized framework, applying consistent normalization, and validating the resulting zones with multiyear yield data and dense soil sampling. The potential of delineating management zones by integrating EMI and NDVI is explored for a single 70 ha agricultural field near Berlin, Germany. Management zones were derived using three data sources: (i) ECa maps from nine different depths of investigation (DOI) obtained with EMI between 2022 and 2024, (ii) seven NDVI images obtained from PlanetScope in 2019, and (iii) a combination of EMI and NDVI data. Management zones were delineated using SOMs, while the optimal number of clusters was obtained with the MCASD method. In a following step, the number of clusters was refined using post hoc analysis using a large dataset of soil samples and yield maps at 10 m resolution from 2011 to 2019. Finally, to what extent management zones derived from EMI, NDVI, or a combination of both represent soil characteristics and yield patterns was evaluated using visual inspection and statistical analysis.

Figure 1 Overview of the patchCROP study area in Tempelberg, Brandenburg (ESRI, 2020). The yellow border indicates the boundary of the investigated field, whereas the green boxes indicate the 30 patches of the patchCROP landscape experiment. The inset map shows the location of the study site within Germany; the red dot indicates the site location in Tempelberg.

2 Materials and methods

2.1 Study area

The study site is part of the patchCROP (patchCROP, 2020) landscape laboratory of the Leibniz Centre for Agricultural Landscape Research (ZALF) near Tempelberg, Brandenburg, Germany (52.4426° N, 14.1607° E; altitude 68 m). It is located in the transition zone between humid oceanic and dry continental climate. The long-term average temperature from 1980 to 2020 was 8.3 °C and the mean annual precipitation for the same period was 533 mm (DWD, 2021; Koch et al., 2023). The investigated field has an area of approximately 70 ha (Fig. 1). Until 2020, this field was managed as a single unit. In March 2020, the patchCROP experiment was established to study the impact of landscape diversification through the use of smaller field sizes, site-specific crop rotations, different field management practices, and the use of new technologies including proximal soil sensing, remote sensing, and robotic technologies (Grahmann et al., 2021). For this, 30 patches of 72 m × 72 m were established within the investigated field (Donat et al., 2022) (Fig. 1). In terms of geomorphology, the site is described as a young moraine landscape shaped by past glaciations and characterized by an undulating relief and heterogeneous soil characteristics (Koch et al., 2023; Öttl et al., 2021; Meyer et al., 2019). The topsoil is predominantly sandy, but a more clayey layer is present at different depths in the subsoil (Hernández-Ochoa et al., 2024).

Figure 2 Workflow diagram showing the integration of proximal (EMI) and remote sensing (NDVI) data for unsupervised clustering using MCASD and SOMs. Yield and soil datasets were used for post hoc validation and refinement of management zones.

2.2 Data collection and processing

The overall methodology of this study is summarized in Fig. 2. This flowchart highlights the role of EMI and NDVI datasets in the clustering process and the use of multiyear yield maps and soil samples for validation and refinement of the resulting management zones. More details are provided in the subsequent sections.

2.2.1 Electromagnetic induction (EMI) measurements

Frequency-domain EMI devices generate a fixed-frequency alternating current in a transmitter coil, which produces a primary magnetic field. This primary magnetic field induces eddy currents in the soil, thus generating a secondary magnetic field. The primary and secondary magnetic fields are sensed by a receiver coil. The quadrature component of the ratio between the primary and secondary magnetic fields is directly proportional to the apparent electrical conductivity (ECa) of the ground (Keller and Frischknecht, 1966; Ward and Hohmann, 1988; McNeill, 1980). The measured ECa is strongly affected by soil properties such as salinity, water content, clay content (and thus texture), compaction, and to a lesser degree organic matter content and cation exchange capacity (Corwin and Lesch, 2005; Robinet et al., 2018). The depth sensitivity of EMI measurements depends on coil spacing and coil orientation. Larger spacing results in increased depths of investigation (DOI), while the coil orientation influences the sensitivity to the shallow or deep subsurface (Lavoué et al., 2010; Simpson et al., 2009).

Table 1 Details of the two EMI devices with coil number, orientation, separation, DOI, and frequency.

In this study, two EMI devices were used simultaneously: a CMD Mini Explorer (GF Instruments, Brno, Czech Republic) with three receiver coils oriented in a vertical coplanar configuration (VCP) and a custom-made CMD Mini Explorer Special Edition equipped with six receiver coils oriented in a horizontal coplanar configuration (HCP). The VCP configuration is most sensitive to the shallow subsurface, with decreasing sensitivity as depth increases. In contrast, the HCP configuration is less sensitive to the shallow subsurface, with sensitivity peaking at a depth of approximately 0.4 times the coil separation (McNeill, 1980). As a rule of thumb, the DOI for the VCP setup is approximately 0.75 times the coil separation. For the HCP setup, the DOI is approximately 1.5 times the coil separation. For the setup used here, the resulting DOI ranges from 0–24 to 0–270 cm. Details of the EMI setup are summarized in Table 1.

Due to the ongoing patchCROP experiment with small patches using variable cropping systems, it was not possible to cover the entire field in a single EMI campaign. EMI data were thus collected in four campaigns conducted between August 2022 and October 2024. During each campaign, the EMI devices were placed in sleds and warmed up for approximately 30 min before use. The sleds were then pulled by an all-terrain vehicle (ATV) at a speed of approximately 6 to 8 km h⁻¹. Data collection occurred at a frequency of 0.2 s, resulting in an inline spatial resolution of 0.25 to 0.50 m. A track spacing of ∼ 2.5 m was used within the experimental patches and a track spacing between 5 and 45 m (typically well below 10 m) was used in the rest of the field. A Real Time eXtended (RTX) center point differential global positioning system (DGPS) (Trimble Inc., Sunnyvale, CA, United States) was used to record the position of the sleds with centimeter accuracy. For more information about the setup for EMI measurements, the reader is referred to von Hebel et al. (2018).

The measured ECa values were filtered using a Python-based method similar to the approach of von Hebel et al., (2014), which has been successfully applied in several studies (Brogi et al., 2019; von Hebel et al., 2021; Kaufmann et al., 2020; Schmäck et al., 2022). The first filter removes values that are deemed too high or too low based on user-defined thresholds (−50 and 50 mS m⁻¹ in this study). A second filter divides the data into a user-defined number of bins (20 in this study) and removes the data from bins with a low fraction of measurements (< 1 % in this study). In a third step, a spatial filter is used to identify and discard ECa values that deviate from adjacent positions by more than a given amount (1 mS m⁻¹ in this study) to avoid unrealistically high lateral ECa variations. After the application of these three filters, ∼ 5 % of the measured ECa values were removed, although this value varied between measurement campaigns.

Given that the EMI data were acquired in four campaigns with different environmental conditions (e.g., soil water content, soil temperature), each EMI acquisition campaign was separately normalized by using a standardized z-score normalization method as used by Rudolph et al. (2015):

$$\begin{array}{}\text{(1)}& {\text{ECa}}_{z,i}=({\text{ECa}}_i-{\mathit{\mu }}_i)/{\mathit{\sigma }}_i,\end{array}$$

where ECa_z,i is the normalized ECa value for the ith campaign, ECa_i is the measured ECa value for the ith campaign, μ_i is the mean ECa value of the ith campaign, and σ_i is standard deviation of ECa values for the ith campaign. Following normalization, manual cleaning was conducted in ArcMap v10.8.2 (ESRI, Redlands CA, USA) to remove points typically occurring at the start and end of each campaign or in short periods where the EMI system was left stationary. In the final step, the normalized data for each of the nine coil configurations were interpolated to a regular 3 m by 3 m grid using ordinary kriging with a Gaussian semivariogram and merged into a single multidimensional raster mosaic dataset.

2.2.2 Remotely sensed NDVI data

High-resolution PlanetScope Level 3B satellite images from the 2019 growing season (winter rye) were used to obtain NDVI maps. Between 1 January 2019 and 31 July 2019, 48 cloud-free images were available. Seven of these images were selected to represent crop development during the growing season. PlanetScope image products are pre-processed and have already undergone radiometric and atmospheric corrections. No additional pre-processing was required. The PlanetScope sensor captures spectral information in four bands: blue (B1), green (B2), red (B3), and near-infrared (NIR – B4) with a spatial resolution of 3 m. The Normalized Difference Vegetation Index (NDVI) was calculated using the reflectance in the red (R) and near-infrared bands (NIR).

$$\begin{array}{}\text{(2)}& \text{NDVI}=(\text{NIR}-R)/(\text{NIR}+R)\end{array}$$

The resulting NDVI values range from −1 to 1, where values close to 1 indicate healthy vegetation, and values close to zero or negative values generally represent non-vegetated surfaces; senescent, stressed, or unhealthy plants or dry vegetation; or features such as clouds and water that exhibit lower NIR reflectance (Sishodia et al., 2020).

2.2.3 Yield data

Georeferenced yield maps of nine growing seasons (2011–2019) were used. These yield maps were generated using a yield monitoring system (CLAAS Quantimeter, Hersewinkel, Germany) mounted on two different combine harvesters. From 2011 to 2013, data were collected using a CLAAS 580. From 2014 onwards, a CLAAS Lexion 770 TT was used. In the 2011–2019 period, the field was cultivated with either winter rye (2011, 2013, 2014, 2016, 2017, and 2019) or rapeseed (2012, 2015, and 2018). For additional details on data processing and yield map generation, readers are referred to Donat et al. (2022). The original yield data from Donat et al. (2022) were available as georeferenced yield data points with a spacing of ∼ 10 m. These points were interpolated to a regular grid with 10 m resolution using ordinary kriging.

2.2.4 Soil sampling and data on soil characteristics

Extensive soil sampling campaigns were conducted between 2020 and 2024, focusing on the experimental patches within the 70 ha field. At 160 locations, soil samples up to 100 cm depth were obtained using a Pürckhauer soil auger with an 18 mm inner diameter. The soil properties analyzed in this study included the depth of soil texture transition, defined as the depth (in cm) at which the sandy top layer ends (EOS (end of sandy layer) in the following), as well as the soil texture (percentage of sand, silt, and clay) of the top sandy layer and the layer below. Soil texture was determined by using the wet-sieving and sedimentation method (DIN ISO, 2002). The particle size distribution was defined according to the IUSS Working Group 150 WRB guidelines (IUSS Working Group, 2015). When multiple subsamples for a single layer were available at a given location, weighted averages of the sand, silt, and clay fraction for the whole layer were obtained using the thickness of each subsample.

2.3 Clustering for delineation of management zones

Three different data combinations were created and investigated: (a) EMI maps, (b) time series of NDVI maps, and (c) a combination of the EMI maps and NDVI maps. Before clustering, a standard pre-processing step of normalization was applied to each dataset to ensure that variables with different ranges and units contribute equally to the classification process. The choice of normalization method can be particularly important when combining datasets with different scales, such as EMI and NDVI, to prevent dominance of one dataset over the other and to maintain the integrity of the input features. In this study, a min–max scaling was applied, where all values were rescaled to a standard range between 0 and 1 (Patro and Sahu, 2015). For EMI, a single normalization was applied to the nine ECa_z maps. In this case, the min–max normalization used the minimum (ECa_z  min) and maximum value (ECa_z  max) from all nine maps:

$$\begin{array}{}\text{(3)}& {\text{ECa}}_z^{\prime }={\displaystyle \frac{{\text{ECa}}_z-{\text{ECa}}_{z\text{min}}}{{\text{ECa}}_{z\text{max}}-{\text{ECa}}_{z\text{min}}}},\end{array}$$

where ECa_z is the original value, and ${\text{ECa}}_z^{\prime }$ is the normalized value. For NDVI, each of the seven NDVI maps was normalized independently:

$$\begin{array}{}\text{(4)}& {\text{NDVI}}_i^{\prime }={\displaystyle \frac{{\text{NDVI}}_i-{\text{NDVI}}_{i,\text{min}}}{{\text{NDVI}}_{i,\text{max}}-{\text{NDVI}}_{i,\text{min}}}},\end{array}$$

where ${\text{NDVI}}_i^{\prime }$ is the normalized value for the ith map, NDVI_i is the original value of NDVI of the ith map, and NDVI_i,min and NDVI_i,max are the minimum and maximum values of the ith NDVI map. This difference in normalization was necessary to preserve the depth-dependent structure of EMI data, as ECa represents a bulk measurement where each reading is influenced by adjacent depths. In contrast, NDVI measurements are independent and acquired at different time points and thus reflect temporal variations in vegetation dynamics.

In this study, a self-organizing map (SOM), an unsupervised machine learning classification technique, was used for clustering (Kohonen, 2013). SOM is a centroid-based clustering technique, similar in some respects to K-means clustering (Celebi et al., 2013). While K-means clustering assigns each data point to a cluster based on the minimum distance to the cluster centroid in the data space, SOM utilizes an artificial neural network to organize and visualize high-dimensional data (Valentine and Kalnins, 2016). The key distinction lies in how SOM projects the data onto a two-dimensional grid while preserving the topological relationships of the input data. Each data vector in SOM is assigned to a numerical cluster, where the cluster center is representative of all the data points associated with it. These cluster centers, which have dimensions similar to the input data vectors, adjust iteratively during the training process to better represent the underlying data distribution. This approach allows SOM to effectively map complex data patterns while maintaining the spatial relationships between clusters.

The multi-cluster average standard deviation (MCASD) approach was used to determine the optimal number of clusters for SOM. This method evaluates the stability of the cluster centers in the data space over multiple clustering attempts as the number of clusters increases. This metric assumes that an appropriate number of clusters for a dataset is any at which the cluster centers do not vary significantly when the clustering algorithm is run multiple times. In this study, MCASD analysis was tested with a maximum number of 20 clusters with 100 SOM clustering runs for each number of clusters to calculate the MCASD stability metric. The number of clustering runs was determined during preliminary testing, where it was observed that most datasets stabilized in terms of variability between 70 and 80 iterations. To ensure consistency and reproducibility, we adopted 100 runs per cluster number. Upon completion of MCASD analysis, the highest number of clusters with a low MCASD metric was selected, as this represents the maximum resolution of the spatial variability that can be obtained through clustering (O'Leary et al., 2023). This clustering process was performed in MATLAB v2023a (MathWorks, Natick, Massachusetts, USA).

2.3.1 Statistical analysis

To assess the differences between clusters derived from the three datasets, a one-way analysis of variance (ANOVA) was conducted in SPSS (IBM, Chicago, IL, United States). ANOVA was used to determine whether there were significant differences between clusters in terms of soil properties or yield using a significance threshold of p < 0.05. Following the ANOVA, a Tukey's HSD (honestly significant difference) test was used as a post hoc analysis to determine which of the clusters were significantly different. In this step, the depth of the sandy layer, the texture of the overlying layer, the texture of the layer below, and the yield data were used. Thus, this step is complementary to the previous cluster selection step with MCASD, which did not consider soil and yield data. Clusters that did not exhibit significant differences were merged during a reclassification step, refining the clustering results to ensure that each final cluster was distinct and statistically meaningful in terms of both the input datasets and soil properties and yield. The latter was confirmed using two-tailed t tests between matching layers of adjacent soil classes in the reclassified map.

3 Results and discussion

3.1 ECa_z, NDVI, and yield maps

The ECa_z, NDVI, and yield maps highlight unique aspects of field heterogeneity and offer insights into subsurface soil properties, aboveground crop performance, and their combined effects on productivity. In the following, these input datasets for management zone delineation are briefly introduced.

Figure 3 Comparison of apparent electrical conductivity (ECa) maps before and after z-score normalization for the HCP 50 cm configuration with (a) the non-normalized ECa map, where the zoomed-in section highlights the influence of varying environmental conditions such as soil moisture and temperature, leading to inconsistent patterns, and (b) the z-score-normalized ECa map, which minimizes the influence of these external factors.

3.1.1 EMI maps

Nine ECa maps with 3 m resolution were obtained from the interpolation of the nine coil configurations recorded during the EMI measurements. The results for one coil configuration (HCP 050 cm) are exemplarily shown in (Fig. 3) before and after normalization. The study area was measured under varying conditions in terms of soil temperature, soil moisture, and effect of agricultural management. This resulted in differences of average ECa and spatial patterns (Fig. 3a). Although it is well known that temperature affects measured ECa (Pedrera-Parrilla et al., 2016; Vogel et al., 2019), it was not possible to perform a comprehensive temperature correction in this study due to the lack of sufficient soil temperature data. Moreover, it has been shown that temperature correction has limitations compared to normalization methods when the dataset is composed of various depths of investigation and is affected by multiple agricultural management practices (Brogi et al., 2019; Rudolph et al., 2015). Thus, z-score normalization was applied for each measurement campaign to reduce the differences between data measured on different days. Figure 3b shows the normalized EMI map for the same coil configuration as shown in Fig. 3a. The normalization successfully harmonized the data, minimizing the influence of varying soil moisture and temperature during acquisition, resulting in more consistent spatial patterns that better represent subsurface soil properties. However, some localized artifacts in the normalized maps still persist. For example, areas near the field boundaries or experimental patches exhibit subtle inconsistencies that may be influenced by edge effects or localized disturbances. Despite these minor limitations, the normalized ECa maps provide a robust foundation for further analysis and management zone delineation.

Figure 4 Normalized apparent electrical conductivity (ECa_z) maps derived from electromagnetic induction (EMI) measurements using multiple coil separations in vertical coplanar (VCP) and horizontal coplanar (HCP) configurations (see Table 1 for more details). These maps highlight the spatial variability of subsurface soil properties, with higher ECa_z values (red) indicating areas of higher moisture retention or finer soil textures and lower ECa_z values (blue) corresponding to sandy soils with lower conductivity.

Figure 4 shows the nine normalized ECa_z maps for the VCP and HCP configurations. These maps display heterogeneous patterns of ECa, primarily attributed to variations in soil characteristics in space and with depth. A prominent feature is the elongated channel extending from the northeast to the southwest of the field, which represents areas with lower ECa_z values. This feature is associated with sandy soils that generally hold less water and nutrients, indicating a coarse-textured zone with lower electrical conductivity. In contrast, the northwest and southeast regions of the field exhibit medium to high ECa_z values, which may reflect areas of higher moisture content and finer soil particles, such as loamy textures. Additionally, in the northeastern part of the field, a more heterogeneous area with short-scale variations can be observed where the ECa_z values vary considerably between the nine maps. For the shallow VCP configurations, this area shows low ECa_z values, which are indicative of sandy soils or dry conditions near the surface. For the deeper HCP configurations, this same area shows higher ECa_z values, suggesting an increase in soil moisture or finer soil texture at greater depths. This pattern highlights the layered soil heterogeneity in this region, with subsurface properties differing significantly from the surface. Overall, the EMI data reveal a high degree of spatial variability and provide valuable insights into subsurface soil variability, which is critical for precision agricultural management.

3.1.2 NDVI maps

All available PlanetScope satellite images for the growing season 2019 (winter rye) were visually evaluated to assess their usability. Before April 2019, no meaningful patterns in NDVI were observed due to the relatively short height (10 to 20 cm) and low biomass of winter rye and the lack of water- or nutrient-induced stress in this early growth stage. Moreover, images from July 2019 were excluded from the analysis as the crop had reached maturity, and no further growth or development was evident. By this time, the physiological activity of the plants had ceased, and harvesting was completed on 4 August 2019.

Figure 5 Seven NDVI maps derived from PlanetScope satellite imagery representing the temporal variability in vegetation development during the 2019 growing season. The images, dated from 5 April to 24 June 2019, capture critical crop growth stages, including flowering and maturity.

Table 2 Summary of remotely sensed NDVI imagery and corresponding dates after sowing.

After this initial analysis, seven NDVI images spanning the period between April and June, hence from flowering to maturity, were selected for further analysis (Fig. 5). The descriptive statistics of the NDVI data are given in Table 2 and show a high degree of temporal variation. Following crop development during the growing season, the mean NDVI peaked on 30 April 2019 (221 d after sowing). Afterwards, NDVI values gradually declined as the crop approached maturity, which is consistent with physiological changes during growth of winter rye (Hatfield and Prueger, 2010). Figure 5 illustrates the temporal development of the spatial variation of NDVI, highlighting the spatial heterogeneity of crop performance within the field (especially Fig. 5d–g), where areas of lower NDVI are associated with poorer crop performance and areas of higher NDVI indicate healthier crops. Generally, the key patterns in crop performance are in good agreement with the patterns observed in the EMI maps. Areas with persistently low NDVI values generally correspond to areas with low ECa_z, and areas with high NDVI values mostly correspond to areas with high ECa_z. However, differences between patterns in NDVI and EMI can also be found. This is expected given that the dynamic changes in crop vigor and vegetation health shown by NDVI are not solely related to subsurface soil conditions captured by EMI. For example, specific areas with low NDVI values were observed in regions of medium to high ECa_z, possibly reflecting localized crop stress due to non-soil-related factors such as disease, waterlogging, or nutrient imbalances.

Figure 6 Nine interpolated yield maps (2011–2019) for the patchCROP field showing spatial variability of crop yield at a 10 m resolution. The maps illustrate yield distributions for winter rye (2011, 2013, 2014, 2016, 2017, 2019) and rapeseed (2012, 2015, 2018). High-yield areas (green) and low-yield areas (red) reflect the inherent field heterogeneity. Variability is observed both within and across years, influenced by crop type, management practices, and environmental conditions. The yield range for each year is provided in decitonnes per hectare (dt ha⁻¹).

3.1.3 Yield maps

Figure 6 presents 9 years (2011–2019) of yield maps interpolated at a 10 m resolution to represent spatial variability across the field. The maps illustrate distinct patterns of high- and low-productivity areas. Yield variability is consistent across multiple years, although variations in measured yield can be observed between years. The years 2012 and 2013 show lower-quality yield data due to incomplete datasets (Donat et al., 2022) caused by equipment issues and environmental challenges during data collection. Despite these limitations, they were retained for spatial context as they still exhibited consistent patterns with other years, and the maps successfully captured the general spatial yield trends and heterogeneity of the field. These years were not weighted differently during validation, and the potential influence of this lack of weighting was mitigated by evaluating multiyear trends and conducting year-by-year comparisons (see Sect. 3.4). The high- and low-yield zones align with known intrinsic field characteristics, such as soil texture, moisture retention, and nutrient availability (Grahmann et al., 2024). These yield patterns will serve as validation for comparing the management zones derived from EMI and NDVI data, as both datasets aim to explain the variability in productivity.

Figure 7 Clustering results for the patchCROP experimental site. (a) MCASD analysis showing appropriate cluster numbers for EMI data. (b) Spatial distribution of original EMI clusters (ESRI, 2020). (c) Spatial distribution of refined EMI clusters after post hoc analysis. (d) MCASD analysis for NDVI data. (e) Spatial distribution of original NDVI clusters. (f) Spatial distribution of refined NDVI clusters after post hoc analysis. (g) MCASD analysis for the combined (EMI + NDVI) dataset. (h) Spatial distribution of the original clusters based on the EMI and NDVI data. (i) Spatial distribution of the refined clusters for the combined dataset after post hoc analysis.

3.2 Clustering of EMI and NDVI

The MCASD analysis for the three datasets provided a robust method to determine the optimal number of clusters (Fig. 7). The analysis suggested a maximum of five clusters for the EMI data (Fig. 7b). These clusters reflect differences in subsurface properties such as soil texture, moisture, and compaction. Cluster 1 corresponds to areas with the highest ECa_z values, which gradually decrease with each subsequent cluster. Cluster 5 represents the lowest ECa_z values. For NDVI (Fig. 7e), a maximum of four clusters was selected. While a five-cluster solution was initially identified as viable for NDVI, increasing the number of clusters beyond four did not significantly reduce variability. This made the four-cluster solution more practical and efficient for representing spatial variability in the NDVI data. Cluster 1 identifies areas with relatively high NDVI values, indicative of healthy, dense vegetation and higher crop performance. NDVI values progressively decrease with higher cluster numbers, with cluster 4 showing the lowest values, representing stressed or less productive areas. The combined EMI and NDVI dataset resulted in four clusters (Fig. 7h). Visual inspection suggests that both the EMI- and NDVI-based patterns are preserved in the combined dataset, likely due to the min–max scaling applied to standardize each dataset before MCASD analysis (see Appendix A). Clusters 1 and 2 represent areas with high values for both ECa_z and NDVI, while cluster 4 identifies zones with low values for both variables, integrating both aboveground and subsurface variability effectively.

3.3 Post hoc analysis

Starting from the optimal number of clusters identified with MCASD, a post hoc analysis based on the nine available yield maps and the point-scale soil samples was conducted. The aim was to verify that the clusters are statistically separated not only in terms of the input data (i.e., EMI, NDVI, or a combination of EMI and NDVI), but also in terms of yield and soil characteristics (i.e., texture of the first and second layers, depth to the second layer). For the EMI-based clusters, 18 soil sampling locations were within cluster 4 and only four of these had an EOS layer within 100 cm depth. The other 14 locations had an EOS layer below the sampling depth of 100 cm and thus no textural values for the lower layer. Thus, the EOS layer depth of cluster 4 was assumed to be below 100 cm, and the texture of the lower layer was excluded from further analysis to have a more consistent characterization of the prevailing soil characteristics.

Table 3 Average values of yield (dt ha⁻¹) and soil properties for the management zones (MZs) derived from EMI, NDVI, and a combination of EMI and NDVI.

NA: not available.

Post hoc analysis indicated that not all clusters were significantly different from each other in terms of yield or soil characteristics. Based on the results of the post hoc analysis, clusters were either left separated when yield or soil characteristics were statistically different (p < 0.05) or grouped together when no statistical separation was identified. For example, clusters 1, 2, and 3 of the EMI-based classification had at least one significant difference in texture, EOS layer, or yield. On the contrary, clusters 4 and 5 did not show statistically significant differences for any of the investigated properties. Thus, clusters 4 and 5 were merged together and the resulting EMI-based cluster map had four clusters with statistically significant separation of input data (i.e., EMI), yield, and soil characteristics. A more detailed breakdown of this post hoc analysis and the resulting merging decisions is provided in Appendix B.

After this post hoc analysis, the resulting refined maps (Fig. 7c, f, and i) have clusters that are statistically separated in terms of the input dataset (i.e., EMI and NDVI) but also in terms of the target variables, which are yield and soil characteristics. Therefore, they are referred to as management zones instead of clusters from this point onwards. These management zone maps appear to be a simplification of the original clustered maps (Fig. 7b, e, and h), but they provide a more holistic understanding of the field by integrating belowground (EMI) and aboveground (NDVI) information with yield and soil data.

3.4 Assessment of management zones derived from different datasets

For each management zone of the maps derived from EMI, NDVI, and a combination of EMI–NDVI, Table 3 shows the average yield between 2011 and 2019 and average soil characteristics, specifically the depth of the soil texture transition (EOS) and the textural fractions (percentages of sand, silt, and clay) of two layers up to 100 cm depth. The average yields of Table 3 vary considerably between different years and follow a general trend of decreasing yields with increasing cluster number. Thus, yields decrease with decreasing ECa_z and NDVI.

Figure 8 Box plots illustrating rye yield (dt ha⁻¹) for 2019 across management zones (MZs) derived from (a) EMI, (b) NDVI, and (c) a combination of EMI and NDVI datasets.

Figure 8 shows the variation in rye yield (dt ha⁻¹) for the management zones derived from different data sources for the year 2019, which is considered representative for most previous years while also allowing a direct comparison with the NDVI data for the 2019 growing season. For the EMI-based management zones (Fig. 8a), the yield distributions for zones 1–3 are relatively similar, with overlapping interquartile ranges and medians. This indicates that the EMI-based management zones are more reflective of subsurface soil properties than yield variability for this particular field. However, zone 4 showed significantly lower yields, corresponding to sandy soils with poor moisture retention (see Table 3). The NDVI-based management zones (Fig. 8b) demonstrate stronger differentiation in yield distribution and a more consistent decline in yield between zones, reflecting the ability of NDVI to capture aboveground vegetation vigor and crop health. In particular, zone 2 reflects an intermediate yield zone between zones 1 and 3, showcasing the ability of NDVI to differentiate changes in crop performance. The management zones derived from combining EMI and NDVI (Fig. 8c) offer narrower interquartile ranges, particularly in zone 2, compared to NDVI-based management zones. This indicates that the integration of EMI and NDVI provides a more consistent and stable representation of yield variability, combining subsurface soil properties with aboveground dynamics. Although NDVI alone offers slightly more pronounced yield differentiation, the combined dataset balances both subsurface and vegetation-related factors effectively, making it a robust approach for management zone delineation. Similar box plots for additional years are provided in Appendix C.

Figure 9 Final management zone maps derived from (a) EMI, (c) NDVI, and (e) a combination of EMI and NDVI datasets. Each zone represents areas with similar subsurface and/or aboveground characteristics. (b, d, f) Corresponding soil profiles for each management zone, detailing soil texture (sand, silt, clay %); dotted lines between zones indicate the depth of textural change (layer 1: above EOS; layer 2: below EOS), and the error bar represents the standard error.

The refined management zones can be associated with a typical soil profile based on the average soil characteristics (Fig. 9). The soil profiles show the textural properties of the first two soil layers and the depth of the interface between these layers (EOS) up to a depth of 100 cm. In some profiles, the EOS layer reaches 100 cm, and thus the textural properties of the second layer are not available. In the case of the EMI-based zones (Fig. 9a and b), zone 1 is characterized by generally higher ECa_z values and identifies areas with a substantial average clay content, especially in the second soil layer (18.6 %). Moreover, the sandier top layer is rather shallow and reaches a depth of around 54 cm. Moving from zone 1 to zone 4, ECa_z generally decreases. At the same time, the depth of the top layer (EOS) becomes deeper while the clay and silt content of the soil decreases and the sand content increases. In zone 4, the average clay content up to 100 cm is 4.1 %, while the sand content is 86.2 %. In the case of the NDVI-based management zones (Fig. 9c and d), the three zones appear to be more indicative of crop development, which results in typical soil profiles with differences that seem less pronounced compared to the case of EMI-based zonation. In this case, NDVI is generally higher in cluster 1 and lowest in cluster 3. The change in soil characteristics between zones follows a similar trend compared to that of EMI-based zones. The depth of the interface between soil layers 1 and 2 increases from 62.7 to 87.4 cm from zone 1 to 3, while the sand content of both layers also increases (from 68.6 % to 87.2 % and 58.1 % to 66.1 %, respectively). The management zones derived from the combined EMI–NDVI dataset (Fig. 9e and f) have typical soil profiles that are similar to those based on NDVI. Also, the sand, silt, and clay contents of the first soil layer appear to be rather similar. However, the range of the depth of the interface between soil layers 1 and 2 is higher for the EMI–NDVI clustered map (63.8 to 100 cm) compared to that of NDVI-based profiles (62.7 to 87.4 cm). At the same time, the difference in texture between the second soil layer of clusters 1 and 2 is stronger in the profiles based on a combination of EMI and NDVI data (see Table 3). These two factors show that the management zones from EMI and NDVI have a relatively high variation between soils of different management zones, which is an improvement compared to the case of the NDVI-based management zones.

Table 4 Results of the pairwise t tests for yield and soil properties between management zones derived from EMI, NDVI, and EMI–NDVI. Bold font indicates significant differences.

NA: not available.

In a final step, statistical validation of the management zones was conducted using pairwise t tests to evaluate the degree of significant differences in yield and soil properties across consecutive zones. The results are summarized in Table 4. A pairwise t test for neighboring zones derived from EMI indicated that the yield of 2012, 2013, and 2016 was not significantly different between zone 1 and zone 2 (p = 0.603, 0.060, 0.253), while the yield of 2012 was not significantly different between zones 2 and 3 (p = 0.209). All other pairwise comparisons indicated significant differences in mean yield. The textural composition of layer 1 was significantly different between all EMI-derived zones. On the contrary, the depth of the top layer was not significantly different between zones 2 and 3 (p = 0.167). In addition, the composition of soil layer 2 was not significantly different between zones 1 and 2 (p of 0.498 for sand, 0.636 for silt, and 0.805 for clay).

The pairwise t tests between neighboring zones based on NDVI indicated that differences in yield among all investigated years were statistically significant. On the contrary, both the depth of the top layer and the composition of soil layer 2 were not significantly different between zones 1 and 2 (p of 0.147 for depth, 0.558 for sand, 0.986 for silt, and 0.627 for clay). These results show that EMI-based zones subdivided the area into one additional class and provided a more comprehensive representation of soil properties up to 100 cm compared to the NDVI-based zones for the investigated field. At the same time, the NDVI-based zones offered a better representation of yield from 2011 to 2019.

The pairwise t test between neighboring zones based on the combined EMI–NDVI dataset showed that the three zones were significantly different for both yield and soil characteristics. This indicates that integrating EMI and NDVI datasets allows for the delineation of zones that are robust in representing both yield variability and soil heterogeneity. Moreover, a visual inspection of the management zone maps (Fig. 9) shows that both maps based solely on EMI or NDVI are affected by west–east-oriented patterns due to measurement direction for EMI and tractor lines in NDVI. These features are not present in the management zone map that integrates EMI and NDVI, suggesting that it also provides a representation of the field that is less affected by external factors. These results underscore the added value of integrating complementary datasets to capture the full spectrum of variability within the field, supporting more informed and effective precision agriculture practices.

3.5 Limitations and perspectives for future work

This study successfully demonstrated the integration of EMI and NDVI datasets for the delineation of management zones, but some limitations are still present and should be addressed in future research. The EMI data were collected during different campaigns under varying environmental conditions (e.g., soil temperature and moisture) and thus required z-score normalization to minimize variability. While effective in this study, this approach may not fully account for certain external factors such as the impact of different management practices in different parts of the field. Similarly, the NDVI dataset was limited to the 2019 growing season as (a) PlanetScope imagery became accessible for this region only in 2019, and (b) the subdivision of the field into differently cultivated patches from 2020 prevented the use of later satellite products. Nonetheless, the choice of PlanetScope imagery (3 m resolution) enabled capturing detailed within-field variability in NDVI, which was particularly important in this study area due to the spatial heterogeneity introduced by soil variation. If coarser-resolution imagery such as Sentinel-2 (10 m) were used instead, smaller-scale patterns in crop development or soil-related variation would have been less detectable due to spatial averaging. This could have reduced the effectiveness of the SOM clustering in identifying distinct management zones. However, for more homogeneous or large-scale fields, Sentinel-2 data could be a practical and freely accessible alternative (Kaya et al., 2025). Another limitation of this study is that the 2019 dataset was considered to be representative of the investigated area. However, a single season of NDVI data may not fully capture interannual variability driven by climatic conditions or crop management practices (Scudiero et al., 2018). Incorporating NDVI data from multiple years in future studies could enable a more comprehensive analysis of temporal dynamics and their impact on management zone delineation to capture yield and soil variability.

A further limitation of the study design was the distribution of soil sampling locations. Although the 160 sampling points provided valuable insights, leveraging EMI-based maps to guide targeted soil sampling could improve spatial representativeness. Additionally, while EMI in this study had a depth of investigation of up to 270 cm, soil sampling was limited to 100 cm depth, potentially missing soil heterogeneity that can affect crops.

Regarding data processing, min–max scaling was a suitable method in this study due to the relatively smooth and filtered input data for both EMI and NDVI. However, this scaling approach is known to be sensitive to outliers and data range extremes (Pedregosa et al., 2011). For datasets with greater variability or different pre-processing methods, alternative scaling approaches such as standardization or robust scaling could be more appropriate (de Amorim et al., 2023). Another factor was the proper application of data normalization prior to clustering, which was essential for obtaining meaningful results in this study (see Appendix A). Future studies should assess the impact of different scaling and normalization strategies on clustering outcomes, especially in settings with noisier or unfiltered sensor data.

In this study, clustering relied on a combination of multi-cluster average standard deviation (MCASD) to determine the optimal number of clusters and self-organized maps (SOMs). While cluster variability was addressed using the MCASD across 100 SOM runs to a large extent, future studies may benefit from incorporating additional stability metrics such as the adjusted rand index (ARI) or cluster overlap measures to better assess classification consistency. The availability of yield and soil data supported the refinement of SOM-based clusters, enabling the merging of groups that were not agronomically distinct. These datasets helped to ensure that the final management zones were both data-driven and interpretable. However, in scenarios where such ground-truth data are limited or unavailable, the initial clusters may still offer useful insights, albeit with greater uncertainty in their agronomic interpretation. Thus, the presented post hoc validation step added confidence in the results but is not strictly required.

The SOM algorithm and the statistical methods used in this study (ANOVA, Tukey's HSD, and t tests) do not explicitly account for spatial autocorrelation, which is inherently present in the interpolated geospatial datasets used here. This may influence statistical outcomes or lead to less spatially coherent clusters in some cases. For instance, kriging interpolation introduces a spatial structure that may challenge the assumption of independence underlying post hoc statistical tests. However, the use of multiyear yield trends and high-resolution soil data helped reduce uncertainty in post hoc validation. Future studies may benefit from incorporating spatially explicit methods, such as spatially constrained clustering, variogram-based diagnostics, or spatial ANOVA, to better account for spatial dependence during both classification and validation stages. In addition to these methodological considerations, future studies should focus on improving the temporal consistency of data collection and increasing the density and depth of soil sampling. The quantification of uncertainty in management zone delineation could also be investigated, for example through ensemble clustering or by incorporating uncertainty from spatial inputs such as EMI interpolation. Finally, long-term monitoring using datasets from multiple years could provide insights into the temporal stability of management zones and their relationship with yield.

The detailed management zone maps complemented with soil characterization obtained in this study should in a next step be integrated into agroecosystem models. This enables simulating and predicting the impact of different management strategies under future environmental and climatic conditions, thus helping to optimize irrigation, fertilization, and other field management practices, further supporting decision-making for sustainable and resource-efficient agriculture.

4 Conclusions

This study integrated proximal soil sensing (EMI) and remote sensing (NDVI) data to delineate high-resolution management zones in a 70 ha agricultural field. Self-organizing maps (SOMs), an advanced unsupervised machine learning technique, were combined with statistical validation methods to identify spatial areas with similar aboveground and belowground properties. Historical yield maps and detailed soil information up to a depth of 100 cm were used to refine and validate the clustering results, ensuring both their accuracy and practical applicability.

To address the variability introduced by environmental conditions during data collection, EMI measurements from multiple campaigns were standardized using z-score normalization, ensuring consistent input for further analysis of the investigated field. Similarly, NDVI data from the 2019 growing season were selected as they represented an uninterrupted crop cycle prior to the subdivision of the investigated field in multiple patches. Before clustering, data were appropriately normalized. The multi-cluster average standard deviation (MCASD) method was applied to determine the optimal number of clusters for different datasets. The optimal number of clusters was determined to be five using the EMI data, four for the NDVI data, and four for the combination of EMI and NDVI datasets. However, statistical validation through Tukey's post hoc analysis using independent yield maps and soil samples reduced the cluster numbers to four, three, and three, respectively. This ensured that the clusters were not only computationally distinct with respect to the input data, but also significantly different in terms of soil characteristics and yield data, thereby increasing their practical relevance in precision agriculture.

Results showed that EMI-based management zones provided a better representation of subsurface properties, particularly soil texture and the depth at which textural changes occur, which underlines the utility of EMI for guiding soil management practices. In comparison, NDVI-based management zones aligned more closely with topsoil characteristics and yield maps, effectively capturing aboveground variability. In general, the integration of EMI and NDVI datasets provided a more comprehensive representation of the spatial variability of both soil characteristics and yield, resulting in management zones that linked both subsurface soil conditions and aboveground vegetation performance. These combined zones effectively explained productivity patterns by bridging the gap between soil properties and crop health.

The product of this study is a high-resolution management zonation map which would provide significant added value in precision and sustainable agriculture. Moreover, it can help in setting up agroecosystem models for the simulation of crop performance and yield and in guiding future soil sampling campaigns. Finally, the workflow proposed in this study can provide a robust blueprint for unsupervised clustering of proximal soil sensing and remote sensing data in agriculture, and future studies should explore the scalability of this methodology in different climatic conditions or other crop systems, as well as investigating additional data sources to further enhance the representation of within-field heterogeneity in soil and crops.

Appendix A: Influence of data normalization

Figure A1 shows a visual comparison of management zone delineation using different normalization approaches. These are (Fig. A1a) EMI-based clustering of ECa_z maps, (Fig. A1b) combined EMI–NDVI clustering with dataset-wise normalization (i.e., normalized by using the minimum and maximum values for all the available data), and (Fig. A1c) combined EMI–NDVI clustering with dataset-wise normalization of EMI data and separate column-wise normalization of NDVI data. As apparent in Fig. A1b, the EMI measurements dominate the clustering results when an inappropriate normalization is used. On the contrary, the normalization strategy used here (Fig. A1c) provides a clustering result where both EMI and NDVI meaningfully contribute.

Figure A1 Comparison of management zone delineation using different normalization approaches: (a) EMI-based clustering without normalization, (b) combined EMI and NDVI clustering with dataset-wise normalization, and (c) combined EMI and NDVI clustering with individual normalization, where EMI data were normalized as a dataset, while NDVI data were normalized column-wise.

Appendix B: Additional results for post hoc analysis

For the EMI dataset (VCP + HCP, nine coils), the MCASD analysis suggested five clusters. The results of the post hoc analysis are shown in Table B1. Statistically significant differences between two clusters are indicated by an O, whereas an X indicates no significant differences. When two clusters have no statistically significant difference for any of the evaluated properties, they are merged. Therefore, clusters 4 and 5 were merged into a new cluster, cluster 4. For the NDVI dataset, the MCASD analysis suggested four clusters and the results of the post hoc analysis (Table B2) merged clusters 3 and 4 into a new cluster, cluster 3. For the combined dataset (EMI + NDVI), the MCASD analysis suggested four clusters and the results of the post hoc analysis (Table B3) merged clusters 1 and 2 into a new cluster, cluster 1.

Table B1 Post hoc analysis of soil characteristics and yield for the EMI-based clusters leading to cluster merging. Statistically significant (O) or nonsignificant differences (X) are provided between clusters for soil texture, EOS layer, and yield.

Table B2 Post hoc analysis of soil characteristics and yield for the NDVI-based clusters leading to cluster merging. Statistically significant (O) or nonsignificant differences (X) are provided between clusters for soil texture, EOS layer, and yield.

Table B3 Post hoc analysis of soil characteristics and yield for the clusters based on EMI and NDVI leading to cluster merging. Statistically significant (O) or nonsignificant differences (X) are provided between clusters for soil texture, EOS layer, and yield.

Appendix C: Differences in yield between derived management zones

Figure C1 presents box plots illustrating yield variability (dt ha⁻¹) for two additional years: winter rye in 2017 (Fig. C1a) and rapeseed in 2018 (Fig. C1b). Results are presented for management zones derived from three clustering approaches: EMI-based (left), NDVI-based (middle), and combined EMI + NDVI (right). These 2 years were selected as additional representative examples, as the overall yield variation across the full 9-year dataset followed the same trend. In the EMI-based management zones, yield distribution is relatively similar across the first three zones, with a noticeable drop in the fourth zone. In contrast, NDVI-based and EMI + NDVI zones show a progressive decline in yield across clusters, indicating a clearer trend of decreasing productivity.

Figure C1 Yield distribution across final management zones based on EMI, NDVI, and combined EMI–NDVI datasets.