Management Zones Delineation through Clustering Techniques Based on Soils Traits, NDVI Data, and Multiple Year Crop Yields

Abstract

Availability of georeferenced yield data involving different crops over years, and their use in future crop management, are a subject of growing debate. In a 9 hectare field in Northern Italy, seven years of yield data, including wheat (3 years), maize for biomass (2 years), sunflower, and sorghum, and comprising remote (Landsat) normalized difference vegetation index (NDVI) data during central crop stages, and soil analysis (grid sampling), were subjected to geostatistical analysis (semi-variogram fitting), spatial mapping (simple kriging), and Pearson’s correlation of interpolated data at the same resolution (30 m) as actual NDVI values. Management Zone Analyst software indicated two management zones as the optimum zone number in multiple (7 years) standardized yield data. Three soil traits (clay content, total limestone, total nitrogen) and five dates within the NDVI dataset (acquired in different years) were shown to be best correlated with multiple- and single-year yield data, respectively. These eight parameters were normalized and combined into a two-zone multiple soil and NDVI map to be compared with the two-zone multiple yield map. This resulted in 83% pixel agreement in the high and low zone (89 and 10 respective pixels in the soil and NDVI map; 73 and 26 respective pixels in the yield map) between the two maps. The good agreement, which is due to data buffering across different years and crop types, is a good premise for differential management of the soil- and NDVI-based two zones in future cropping seasons.

1. Introduction

Field management zones in precision agriculture (PA) are defined as “the sub-regions within the same piece of land showing similar yield influencing factors within which different crop management practices are carried out at the right time and place to optimize crop productivity and minimize adverse environmental impact” [1]. Well delineated homogenous areas should have the same yield trend across years and different cultivated crops [2]. These management zones are not only used for site-specific, real-time crop management, but also for soil sampling to characterize the soil nutrient levels for variable rate fertilization [3]. Under uniform crop management practices, fertilizer use efficiency is quite low [4,5]. Therefore, efficient fertilizer management practices as those under PA are expected to enhance the efficiency of nutrients supplied with fertilizers, as a premise for improved crop yield. Specifically, crop input rates, i.e., fertilizers, should be applied based on the crop requirements, estimated through the information of within-field spatial variability, such as nutrient status or productivity potential, within a specific region of the field [6,7]. It is recognized that there is often high spatial variability within the same field, due to nutrients or other factors or interaction among them, which influence the final crop yield. The high or low yield factors could be soil, environmental, topographical, anthropogenic, or biological entities [8]. Specifically, the interaction of various factors makes it more complex to understand the real story behind them. Therefore, efficient methods should be used to measure the within-field spatial variability in an appropriate way for delineating site-specific zones [9].

Geostatistical analysis has been extensively used to estimate the spatial variability of crop yields, soil, and remote sensing data using different kriging methods with small to larger size grids [10]. Generally, geostatistics is used to estimate the spatial autocorrelation within a dataset in terms of variogram parameters and spatial mapping. More to this, spatial variability within a dataset is estimated through best-fit modelling [11].

Various methods have been developed to measure the within-field spatial variability for delineating management zones. Most of them rely on single types of data for analysis of multiple year yield mapping [12,13,14], various soil properties [15,16,17], remotely sensed data [12,18,19,20,21], topographic factors [22], and soil apparent electrical conductivity (ECa) [23,24]. All these methods, directly or indirectly, are related to crop productivity [25]. As PA continuously evaluates the different methods for setting spatial management zones, and no generally accepted method exists yet, cluster analysis may be a basis for unresolved questions in delineating management zones [26]. This statistical cluster analysis uses various types of data sources to develop within-field areas at similar characteristics (soil and crop-related parameters or environmental impact) for the quantification of spatial variability trends and to optimize crop productivity through site-specific crop management practices [27]. Franzen and Kitchen (1999) [28] used different typologies of data such as topography, satellite imagery, soil electrical conductivity, crop yields, and soil properties, to develop potential zones for variable nitrogen rates. However, delineating management zones based on crop yields, soil, and satellite data have all been considered as logical methods in PA [2,29,30].

Based on this, in a field subjected to arable cropping with a georeferenced record of crop yield data for a number of years, the objectives of this study were: (1) to characterize the spatial variability of crop yields, soil properties, and remotely sensed NDVI data influencing crop productivity through geostatistical analysis; (2) to delineate feasible crop management zones with a fuzzy clustering method, using previous crop yields, soil properties, and NDVI data; (3) to interpret the obtained management zones for site-specific management in view of future crop management practices.

2. Materials and Methods

2.1. Study Site

The study was conducted on an 8.94 ha area located near Budrio, Bologna, in Northern Italy (49°43’41” N, 70°67’52” E, 15 m above sea level). The area falls in the Mediterranean North Environmental Zone [31]. The uppermost field part has a sandy texture, compared to the central and lower parts where the soil is richer in clay. The experimental area falls in a temperate–warm climate zone, where average temperature is around 15 °C on a yearly basis, whereas precipitation is mostly registered in the cold semester.

2.2. Data Collection

2.2.1. Crop Data

In this field’s history, georeferenced yield data were collected at harvest for seven years between 2005 and 2020. They were three years of winter cereals and four years of spring/summer cereals, also including a dicot. The period considered (7 years of data over 16 years) derives from the availability of georeferenced yield data, limited only to crops that are harvested with a machine equipped with a yield monitoring system. The specific crops, their seeding and harvesting dates, and subsequent crop cycles are reported in

Table 1. Seeding and harvesting dates, and crop duration of cultivated crops.

Crop & Year	Botanical Name	Seeding	Harvesting	Duration (Days)
BW 2005	Triticum aestivum L.	19 October (2004)	24 June	248
SU 2006	Helianthus annuus L.	10 April	12 September	155
BW 2007	Triticum aestivum L.	9 October (2006)	18 June	252
BW 2009	Triticum aestivum L.	21-October (2008)	30 June	252
SO 2010	Sorghum bicolor L.	21 April	28 August	129
MA 2018	Zea mays L.	18 June	25 September	99
MA 2020	Zea mays L.	28 March	8 August	133

BW, bread wheat; SU, sunflower; SO, sorghum; MA, maize.

. Good management practices were used for the cultivation of each specific crop, based on local conditions.

For five years, grain yield (GY) data were collected by a New Holland CR 9080 (CNH Industrial N.V., Basildon, UK), equipped with an assisted guiding system based on real-time kinematic GPS, yield mapping, consisting of a Pektron flow meter (Pektron Group Ltd., Derby, UK), and an Ag Leader moisture sensor (Ag Leader Technology, Ames, IA, USA). Among seven years, maize crop was cultivated for biomass production during cropping seasons of 2018 and 2020. The crop was harvested with a New Holland FR Forage Cruiser 780, fully integrated with auto-guidance package RTK correction signals. The sensor of feed roller displacement measured crop throughput, giving the yield data. The resistive-type moisture-sensing system provided real-time moisture data displayed and registered in the in-cab IntelliView™ system. Raw yield data were filtered using Yield Editor software (Version 2.0.7; USDA-ARS, Columbia, MI, USA), outliers were removed, and moisture was adjusted at 13% for wheat, 9% for sunflower, 14% for sorghum and 0% for maize. The final yield points of each crop were intersected with the field boundary layer. Then average of GY data points were 12152 for wheat per year, 6455 and 6575 for sunflower and sorghum, respectively, whereas an average of 27,808 data points for maize per year were retained for further analysis.

2.2.2. Soil Sampling

Spatial soil variability was determined through the grid soil sampling technique. At each grid, one representative soil sample was collected at 0.0–0.3 and 0.3–0.6 m soil depth, and the position of each sampling site was recorded with a Trimble GPS. In this way, 21 soil samples were collected to assess the soil physical-chemical properties over the entire field during the year 2019, as shown in Figure 1. A manual soil probe was used for collecting the soil samples. At the end, all samples were air-dried at 40 °C and sieved to 2 mm diameter.

2.2.3. NDVI Index

Remotely sensed NDVI data from a Landsat 5-Thematic Mapper (TM), Landsat 7-Enhanced Thematic Mapper Plus (ETM+) and a Landsat 8-Operational Land Imager (OLI) were downloaded using US Geological Survey (USGS)-Earth Explorer platform. Collection Level-1 (C1L1), Tier 1, Precision Terrain (LITP) was the adopted inventory structure that provided the highest quality data for the time-series records, with 16 days average revisiting time (temporal resolution). More, the USGS-Landsat data collection platform was used, having 30 m spatial resolution, for vegetation monitoring during the growing period of surveyed crops [32]. Landsat satellite scene selection was done under a clear sky and good quality of pixels at field level. Multispectral NDVI time-series data were downloaded over the entire following prospections: 28 February–15 June for BW 2005, 20 April–20 August for SO 2006, 21 February–20 May for BW 2007, 15 February–30 May for BW 2009, 11 May–29 July for SO 2010, 9 July–21 September for MA 2018, 28 April–26 July for MA 2020. Depending on available data, different Landsat missions were used during effective crop growing periods as shown in

Table 2. Prospection of remotely sensed NDVI imagery and corresponding dates after sowing (DAS).

Crop & Year	LS- Mission	Date	DAS
BW 2005	LS-7 ETM+	15 April	178
	LS-5 TM	23 April	186
	LS-7 ETM+	1 May	194
	LS-7 ETM+	2 June	226
SU 2006	LS-5 TM	26 April	16
		15 July	96
		31 July	112
		16 August	127
BW 2007	LS-7 ETM+	4 March	146
	LS-7 ETM+	5 April	178
	LS-5 TM	29 April	202
	LS-5 TM	15 May	218
BW 2009	LS-7 ETM+	21 February	123
	LS-5 TM	18 April	179
	LS-5 TM	20 May	211
SO 2010	LS-5 TM	8 June	48
		24 June	64
		10 July	80
MA 2018	LS-7 ETM+	24 July	36
	LS-7 ETM+	9 August	52
	LS-8 OLI	17 August	60
	LS-7 ETM+	10 September	84
MA 2020	LS-7 ETM+	26 May	59
	LS-8 OLI	19 June	83
	LS-8 OLI	5 July	99
	LS-8 OLI	21 July	115

BW, bread wheat; SU, sunflower; SO, sorghum; MA, maize; LS, Landsat; DAS, days after sowing.

. A total of 26 images during the study period were retained as being the most interesting ones (Table 2), and used in subsequent analysis. The NDVI index was calculated through the near-infrared and red portion of the electromagnetic spectrum, based on the formula shown in Equation (1):

$$\mathrm{NDVI}=\frac{\left(R_{NIR}-R_{Red}\right)}{\left(R_{NIR}+R_{Red}\right)}$$

(1)

where, $R_{NIR}$ and $R_{Red}$ represent the reflectance values of near-infrared and red (visible), respectively.

2.3. Analysis Methods

2.3.1. Descriptive Statistics

All surveyed datasets were subjected to descriptive statistics to characterize them based on their mean, minimum, maximum, standard deviation, skewness, and kurtosis. In this study, averaged soil depth data (0–0.6 m) was used, resulting from the separate analysis of 0–0.3 and 0.3–0.6 m depth samples [14].

2.3.2. Soil Analysis

The 42 soil samples resulting from the 21 positions per 2 depths (0–0.3 and 0.3–0.6 m) were subjected to the following determinations: sand, silt, clay, pH, total limestone (CaCO₃), organic carbon (C), total nitrogen (N), and electrical conductivity of the soil saturated extract (ECe) with a 1:2.5 (w/w) soil-to-water ratio. The soil texture was determined by the pipette method [33]. Soil pH was determined in a 1:2.5 (w/w) soil-to-water ratio. Total limestone was determined volumetrically [34]. Total organic C and total N were determined by a CHNS-O Elemental Analyser 1110, Thermo Scientific GmbH, Dreieich, Germany. The complete dataset is reported in the Supplementary Materials (Table S2).

2.3.3. Geostatistical Analysis

The geostatistical analysis was performed on crop yields, soil, and NDVI data to explain the spatial dependence (SpD) in their datasets in relation to the variogram parameters: nugget (C₀), the error at 0 distance; sill (C₀ + C), the maximum value at y-axis that increases with increasing lag distance (h); range (a), the maximum distance at which data points are still correlated. The degree of SpD was calculated by nugget-to-sill ratio and interpreted according to Cambardella et al. [35]: it was indicated as being strong with SpD < 25%; intermediate, in the range 25–75%; weak with SpD > 75%. To archive the best prediction accuracy, we employed the iterative cross-validation technique, seeking the highest coefficient of determination (R²) and minimal error parameters: mean absolute error (MAE), root mean square error (RMSE), and absolute standard error (ASE). Based on this, we used the best fitting semi-variogram model among the circular, spherical, exponential, Gaussian, and stable models [36].

2.3.4. Spatial Maps

Spatial maps of crop yields, soil, and NDVI data were delineated by simple kriging (SK) interpolation with a 10 m grid size. SK was chosen because it gives maximum R² and minimum error parameters [36].

For yield mapping over single and multiple crops, standardized interpolated data laid over 937 grid points were used, through converting the actual GY (t/ha) into a relative percentage, which gives an 100% field average over entire dataset [14]. The procedure is given below:

Over a single crop, the formula of Equation (2) was used to delineate the spatial maps over a single crop of the specified year:

$$S_i=\left(\frac{y_i}{\overline{y}}\right)\times 100$$

(2)

where, $S_i$ = standardized yield (percentage) at point (i), $y_i$ = interpolated yield at point i (t/ha), and $\overline{y}$ = average interpolated yield (t/ha) over the entire dataset across the field.

Over multiple crops, the spatial map was delineated by calculating the averaged values (relative percentage) laid over seven years of crop yield with this formula, shown in Equation (3).

$${\overline{S}}_i=\frac{{{\displaystyle \sum }}_{t=1}^{n}S{i}_I}{n}$$

(3)

where, $\overline{S}$_i = mean interpolated yield (percentage) over n years, i.e., as the seven-year yield, at point (i); S$i_I$ = interpolated yield (percentage) at point (i), summed over the seven years surveyed, $n$ = number of years surveyed (=7).

Over multiple crops, the temporal map was delineated based on the “coefficient of variation (CV %)” of each grid point over seven years of crop yield (Equation (4)) [37]. However, a temporal variability map was produced to determine the consistency of crop yields over seven years of cropping history.

$$CV{S}_i=\frac{{\left(\frac{\left(n{{\displaystyle \sum }}_{t=1}^{t=n}S{i}_t{}^2-{\left({{\displaystyle \sum }}_{t=1}^{t=n}S{i}_t\right)}^2\right)}{n\left(n-1\right)}\right)}^{0.5}}{\overline{S_i}}\times 100$$

(4)

where, $CV{S}_i$ = coefficient of variation of the standardized yield at point (i) over n years; Si_t = standardized yield (percentage), at point (i); $\overline{S_i}=$ mean standardized yield at point (i).

To define the similarity among spatial maps, they were developed with five classes, using the equal interval classification method. Map production and geostatistical data analysis were carried out with the ArcGIS software (Version 10.3, ESRI, Redlands, CA, USA) under the reference system WGS 84/UTM zone 32 N.

2.3.5. Pearson’s Correlation

For summarizing and data clustering purpose, we aggregated the interpolated crop and soil data with a 30 m wide grid to a similar extent of NDVI imagery (30 m resolution). The attribute point layer of the measured data was joined with the polygon grid by using the mean function to calculate the mean values fallen in each grid, resulting in 99 data points. Thereafter, we determined the Pearson’s correlations at the same layer extent between yield, soil, and NDVI data.

2.3.6. Delineation of Management Zones

After geostatistical and correlation matrices, the fuzzy c-means clustering algorithm method was used to classify the data and evaluate the feasible number of clusters by the Management Zone Analyst (MZA) software [30,38]. Using this specific software, the following clustering parameters were defined:

• Measure of similarity: Euclidean for univariate, Diagonal for independent and Mahalanobis for multivariate
• Fuzziness exponent = 1.30.
• Maximum no. of iterations = 300.
• Convergence criterion = 0.0001.
• Maximum no. of zones = 6.
• Minimum no. of zones = 2.

The optimum number of clusters is achieved at minimal values of the two indices FPI (fuzziness performance index) and NCE (normalized classification entropy), which together indicate data dispersion [38]. Over seven years of crop yield data, FPI and NCE values were minimized when classified into two clusters (

Table 3. Values of FPI and NCE based on standardized crop yields over seven years for the delineation of 2, 3, and 4 zones.

Number of Clusters	Crop Yields
	FPI	NCE
2	0.0369	0.0151
3	0.1007	0.0536
4	0.1342	0.0809

FPI (fuzziness performance index), NCE (normalized classification entropy).

). Hence, the two-zone option was selected based on the minimal FPI and NCE indices, because more than two zones were redundant.

In order to select the parameters of soil and NDVI indices that could be the best candidates for delineating the management zones, based on the multiple soil traits and NDVI indices, we selected the three soil traits (clay, CaCO₃, and N), which showed the best correlations with multiple crops’ yield data. For the selection of NDVI data, we chose data from one acquisition date per year, which presented high and significant correlations with yield data in each specific year (15 April 2005, 29 April 2007, 21 February 2009, 24 June 2010, and 24 July 2018). However, no NDVI data was chosen for SU 2006 and MA 2018 due to insignificant correlations.

Thereafter, we developed two-zone maps based on each of the eight soil traits and NDVI data, using the MZA clustering method. We went further and developed a final two-zone map based on best multiple soil traits and NDVI data: dealing with datasets having different units, we used normalized data (average = 0, SD = 1) to remove the weight of each parameter during cluster analysis. Similarly, a two-zone map over seven years of multiple crops’ yield was also developed, which was used as a reference map.

Since the main goal of this study was to optimize crop productivity, the two-zone maps based on single or multiple parameters were compared with the two-zone yield map (reference map). According to the procedure used, each pixel of the map (independent variable/s) was compared with the same cell of the two-zone yield map (Figure 2). Pixel agreement was calculated as being the number of pixels belonging to the same zone (i.e., number of pixels belonging to a high yield zone (e.g., map of SU_2006) and a high single soil property zone (e.g., map of carbon zones), plus the number of pixels belonging to low yield zone and a low soil property zone, divided by the total number of pixels, to provide an index of the overall accuracy, i.e., a degree of spatial agreement between the compared data. In contrast to this, disagreement was established if there was no zone consistency between pixels of independent variable/s and yield maps.

2.3.7. Weather Data

The weather data was downloaded by the nearest station of the Hydro-meteorological Service of the Emilia-Romagna Region, under the ARPAE datahub platform. The wet and dry periods during the growing seasons of the surveyed crops were calculated as the difference between precipitation (p) plus, only in MA 2018, irrigation, and crop evapotranspiration (ET_C), this latter determined according to Allen et al. [39]. The data are reported in the Supplementary Materials (Table S1).

2.3.8. Statistical Analysis

Crop yields, soil traits, and NDVI data were subjected to descriptive statistics. Pearson’s correlation (r) was used to evaluate the relationships between soil properties and NDVI data on one side, and single and multiple yields on the other side. One way analysis of variance (ANOVA) was run with the Statistica 10 software (StatSoft Corp., Tulsa, OK, USA) to assess the differences in soil, NDVI and yield data in the two-zone maps. The least significant difference (LSD) at p < 0.05 was used to indicate significant differences between the high and low zone.

3. Results

3.1. Crop Yields: Descriptive Statistics, Spatial Variability, and Spatial Maps

The descriptive statistics of the standardized crop yields in the seven years indicates a sizable variation around the mean yield for all the crops (

Table 4. Descriptive statistics of standardized crop yields over seven years (means = 100).

Source	Median	Min	Max	SD	Skewness	Kurtosis
BW_2005	102	43	125	12.6	−1.23	2.84
SU_2006	100	46	147	14	−0.11	0.56
BW_2007	101	57	140	12.8	−0.45	0.69
BW_2009	102	52	133	12.9	−0.8	0.75
SO_2010	100	66	208	10.4	1.21	12.5
MA_2018	101	79	107	4.4	−1.53	3.21
MA_2020	101	29	149	12.6	−0.88	2.38

Actual mean yield (t/ha); BW_2005 (6.91), SU_2006 (2.00), BW_2007 (3.50), BW_2009 (6.59), SO_2010 (2.71), MA_2018 (14.52), MA_2020 (24.5); BW, bread wheat; SU, sunflower; SO, sorghum; MA, maize; SD, standard deviation.

). Standard deviation (SD) was between 4.4 (MA_2018) and 14% (SU_2006) of the mean.

The pattern of crop yield spatial distribution indicates that, for five crops, the stable model was the best fitting, while for the remaining two crops, two different models (spherical and exponential) were the best fitting (

Table 5. Semi-variogram parameters of standardized crop yields over seven years.

Source	Model	C0	C0 + C	a (m)	C0/(C0 + C) %	SpD	R2	MAE	RMSE	ASE
BW_2005	Stable	27.8	140	29	20	S	0.96	0.028	9.75	10.5
SU_2006	Stable	40.6	330	146	12	S	0.95	−0.0035	10.1	12.5
BW_2007	Stable	0	283	52	0	S	0.90	0.0061	10.7	12.3
BW_2009	Spherical	24.8	164	10	15	S	0.94	0.021	7.87	7.78
SO_2010	Stable	0	212	87	0	S	0.92	0.098	6.69	10.3
MA_2018	Exp.	0	13	25	0	S	0.87	0.0004	1.46	1.92
MA_2020	Stable	30.7	189	22	16	S	0.87	−0.08	6.81	7.56

BW, bread wheat; SU, sunflower; SO, sorghum; MA, maize; C₀, nugget; C, partial sill; C₀ + C, sill; a, range; C₀/(C₀ + C), nugget-to-sill ratio; SpD, spatial dependence; S, strong; R², coefficient of determinants, MAE, mean absolute error; RMSE, root mean square error; ASE, absolute standard error.

). Based on the relationship between nugget and sill, a strong spatial dependence was evidenced in all cases up to the range, which varied from a minimum of 22 m (MA_2020) to a maximum of 146 m (SU_2006). The parameters of prediction accuracy indicate an overall good fitting (R²) and low error in terms of the three parameters, MAE, RMSE, and ASE (Table 5). The overall good fitting is reflected in all the fitted models closely interpolating measured points (Figure S1 in the Supplementary Materials).

The spatial maps reveal consistently low yielding areas near the SW and NE borders, whereas the high yielding areas are more variably dispersed across the field in the seven years (Figure 3a). The multiple crops map averages the single year patterns by indicating lower yielding areas near the SW and NE borders, and higher yielding areas quite evenly spread in the central part of the field (Figure 3b). Higher yields were associated with lower inter-annual variability (lower CV), with respect to lower yields (Figure 3b).

3.2. Soil Traits: Descriptive Statistics, Relationships with Yield, Spatial Variability, and Spatial Maps

On average, sand and silt particles predominated in the soil at 0–0.6 m of depth (

Table 6. Descriptive statistics of soil traits (depth 0–0.6 m).

Soil Variables	Mean	Median	Min	Max	SD	Skewness	Kurtosis
Sand (g/kg)	549	504	335	793	147	0.22	−1.32
Silt (g/kg)	354	379	153	529	116	−0.16	−1.33
Clay (g/kg)	97	110	35	156	37.9	−0.28	−1.39
pH	7.93	7.90	7.72	8.19	0.13	0.72	−0.1
CaCO3 (g/kg)	144	145	116	162	11.1	−0.71	0.66
C (g/kg)	13.5	14.0	8.7	17.1	2.33	−0.4	−0.54
N (g/kg)	1.62	1.68	1.08	2.06	0.28	−0.32	−0.7
C:N	8.34	8.40	7.96	8.81	0.23	0.07	−0.85
ECe (ds/m)	0.42	0.41	0.26	0.63	0.1	0.6	−0.39

Min, minimum; Max, maximum; SD, standard deviation.

), although all three particle size classes (sand, silt, and clay) staged a noticeable variability (SD ranging from 27 to 42% of the mean). The data of the two soil depths (0–0.3 and 0.3–0.6 m) are shown in Table S2 of the Supplementary Materials. The pH was always close to 8, which was consistent with a significant, quite stable content of total limestone (CaCO₃). The amount of total organic carbon and total nitrogen was within the range of the average data of this area. The C:N ratio, below 10, indicates a soil prone to organic matter mineralization. Lastly, soil salinity (ECe) was barely detectable, which is consistent with this soil being far away from sources of salinity, such as the seacoast or wells/streams supplying brackish water. The semi-variogram parameters of soil traits are reported in

Table 7. Semi-variogram parameters of soil traits.

Source	Model	C0	C0 + C	a (m)	C0/(C0 + C)%	SpD	R2	MAE	RMSE	ASE
Sand	Spherical	0	717	119	0	S	0.98	1.16	12.29	15.48
Silt	Exp.	0	531	125	0	S	0.96	−0.78	12.22	17.06
Clay	Exp.	0	39.5	120	0	S	0.95	−0.20	3.83	4.73
pH	Exp.	0	0.00003	90	0	S	0.92	−0.00011	0.0042	0.0045
CaCO3	Exp.	0	0.48	139	0	S	0.96	−0.012	0.40	0.49
C	Exp.	0	0.16	118	0	S	0.94	0.0034	0.26	0.30
N	Exp.	0	0.0024	130	0	S	0.94	0.00092	0.032	0.036
C:N	Spherical	0	0.0035	140	0	S	0.94	−0.0020	0.026	0.031

C₀, nugget; C, partial sill; C₀ + C, sill; a, range; C₀/(C₀ + C), nugget-to-sill ratio; SpD, spatial dependence; S, strong; R², coefficient of determinants, MAE, mean absolute error; RMSE, root mean square error; ASE, absolute standard error.

.

The relationships between soil traits and crop yields in the single and multiple crops point out some traits as having significant correlations with yield in many single years and in the seven-year combination (

Table 8. Pearson’s correlations between soil traits and crop yields (n = 99).

Soil Traits	BW_2005	SU_2006	BW_2007	BW_2009	SO_2010	MA_2018	MA_2020	Mean Yields
Sand	−0.69	−0.09	−0.19	−0.17	−0.08	−0.58	−0.07	−0.44
Silt	0.67	0.08	0.20	0.14	0.08	0.59	0.07	0.43
Clay	0.71	0.23	0.17	0.34	0.04	0.37	0.10	0.47
pH	−0.68	0.15	−0.08	−0.06	0.16	−0.53	0.24	−0.26
CaCO3	−0.18	0.29	0.38	0.37	0.64	0.14	0.43	0.36
Carbon	0.09	0.12	0.44	0.25	0.17	0.07	0.61	0.30
Nitrogen	0.15	0.11	0.47	0.23	0.20	0.16	0.60	0.34
C/N	−0.40	0.10	−0.12	0.15	−0.16	−0.56	0.08	−0.22
ECe	−0.01	0.10	0.36	0.17	0.10	−0.04	0.48	0.18

BW, bread wheat; SU, sunflower; SO, sorghum; MA, maize; bold letter correlations are significant at p < 0.05 according to Pearson correlation.

). This is the case for clay content, total limestone, and total nitrogen. Another trait, ECe, shows significant correlations with yield only for two years and is not significantly correlated with multiple, i.e., seven-year, crop yield. Based on this, ECe was dismissed in the subsequent assessment of spatial variability and the drawing of spatial maps.

The interpolated maps of soil traits exhibit a variable picture (Figure 4): sand shows a spatial distribution of high and low values, which is approximately opposed to that of silt and clay; carbon and nitrogen have the same spatial pattern; the rest of traits outline spatial distributions scarcely related to the previous traits.

3.3. NDVI Indices: Descriptive Statistics, Relationships with Yield, Spatial Variability, and Spatial Maps

Three or four NDVI dates were retained in each crop for the present work. They were mostly acquired in the central stages of the seven crops (

Table 9. Descriptive analysis of NDVI indices.

Crop & Year	Acquisition Date	Mean	Median	Min	Max	SD	Skewness	Kurtosis
BW 2005	15 April	0.69	0.70	0.53	0.73	0.03	−2.22	7.67
	23 April	0.47	0.49	0.22	0.52	0.06	−2.61	6.19
	1 May	0.76	0.76	0.7	0.79	0.02	−1.33	1.5
	2 June	0.5	0.50	0.41	0.55	0.03	−0.42	−0.81
SU 2006	26 April	0.16	0.16	0.13	0.28	0.03	2.04	5.01
	15 July	0.55	0.56	0.47	0.61	0.03	−0.4	0.51
	31 July	0.41	0.41	0.33	0.45	0.02	−0.96	1.12
	16 August	0.32	0.33	0.1	0.57	0.13	0.01	−1.25
BW 2007	4 March	0.72	0.73	0.58	0.75	0.03	−2.35	6.31
	5 April	0.84	0.84	0.78	0.86	0.02	−1.73	3.29
	29 April	0.67	0.67	0.61	0.68	0.01	−1.76	3.55
	15 May	0.61	0.61	0.58	0.64	0.01	−0.47	−0.02
BW 2009	21 February	0.54	0.55	0.44	0.61	0.03	−0.82	1.81
	18 April	0.68	0.67	0.57	0.73	0.03	−1.18	4.13
	20 May	0.61	0.62	0.56	0.64	0.01	−1.77	5.34
SO 2010	8 June	0.48	0.49	0.34	0.54	0.04	−1.14	1.79
	24 June	0.72	0.73	0.67	0.75	0.01	−1.46	2.58
	10 July	0.72	0.73	0.68	0.74	0.01	−2.47	8.54
MA 2018	24 July	0.59	0.60	0.31	0.66	0.06	−2.75	10.09
	9 August	0.64	0.65	0.41	0.69	0.05	−1.39	2.49
	17 August	0.78	0.81	0.41	0.84	0.07	−2.6	7.72
	10 September	0.56	0.59	0.31	0.63	0.07	−1.86	3.05
MA 2020	26 May	0.55	0.58	0.3	0.65	0.08	−1.55	1.67
	19 June	0.84	0.88	0.37	0.9	0.11	−2.93	8
	5 July	0.81	0.85	0.37	0.87	0.11	−2.75	7.13
	21 July	0.8	0.83	0.4	0.85	0.09	−2.84	7.64

BW, bread wheat; SU, sunflower; SO, sorghum; MA, maize; Min, minimum; Max, maximum; SD, standard deviation.

). The mean values are the most varied, depending on more or less favorable ambient conditions, as well as plant stage: crops approaching senescence displayed generally lower NDVI values. However, in SU_2006, NDVI values never attained the level of 0.60 during the whole crop cycle.

The relationships between NDVI data and single crop yields evidence several dates in specific crops where NDVI data are significantly correlated with yield (

Table 10. Pearson’s correlations (r values) between NDVI data in different acquisition dates and crop yields (n = 99).

BW_2005	SU_2006	BW_2007	BW_2009	SO_2010	MA_2018	MA_2020
0.79 (15 April)	−0.30 (26 April)	0.53 (4 March)	0.79 (21 February)	0.29 (8 June)	0.55 (24 July)	0.16 (26 May)
0.05 (23 April)	−0.17 (15 July)	0.57 (5 April)	0.48 (19 April)	0.53 (24 June)	0.36 (9 August)	0.18 (19 June)
0.71 (1 May)	0.06 (31 July)	0.58 (29 April)	0.54 (20 May)	0.31 (10 July)	0.32 (17 August)	0.20 (5 July)
0.58 (2 June)	0.18 (16 August)	0.31 (15 May)			0.24 (10 September)	0.22 (21 July)

BW, bread wheat; SU, sunflower; SO, sorghum; MA, maize; correlations in bold are significant at p < 0.05.

). However, in the case of SU_2006 and MA_2020, there is no NDVI data significantly correlated with the final yield. Based on insignificant correlations with yield, NDVI from these two crops and from 23-Apr in BW_2005 were dismissed in the subsequent assessment of spatial variability and the drawing of spatial maps.

In the remaining 17 cases, the spatial pattern of NDVI values was best described by an exponential model (

Table 11. Semi-variogram analysis of NDVI index and best fit exponential model performance.

Crop	Acquisition Date	C0	C0 + C	a (m)	C0/(C0 + C) %	SpD	R2	MAE	RMSE	ASE
BW 2005	15 April	0	0.87	66	0	S	0.99	0.001	0.017	0.012
	1 May	0.58	1.01	131	57	M	0.85	0.00049	0.012	0.0098
	2 June	0	1.14	102	0	S	0.90	0.0002	0.01	0.01
BW 2007	4 March	0.13	0.94	69	14	S	0.98	0.00086	0.017	0.015
	5 April	0.005	1.01	78	0	S	0.97	0.00069	0.0098	0.0085
	29 April	0	1.01	82	0	S	0.96	0.00056	0.01	0.0089
	15 May	0	1.11	92	0	S	0.95	0.00002	0.011	0.01
BW 2009	21 February	0	1.03	56	0	S	0.98	0.0007	0.017	0.014
	18 April	0	1.11	92	0	S	0.98	0.00003	0.011	0.0102
	20 May	0.57	0.96	56	59	M	0.87	0.00015	0.0083	0.0071
SO 2010	8 June	0	1.2	102	0	S	0.99	0.00006	0.02	0.0201
	24 June	0.36	0.92	56	39	M	0.88	0.00056	0.01	0.0081
	10 July	0.61	0.92	59	67	M	0.65	0.0003	0.0086	0.0069
MA 2018	24 July	0	1.09	130	0	S	0.99	0.0024	0.031	0.023
	9 August	0	1.01	110	0	S	0.99	0.001	0.019	0.016
	17 August	0.12	1	83	12	S	0.99	0.0033	0.038	0.029
	10 September	0	0.95	63	0	S	0.99	0.0019	0.026	0.022

BW, bread wheat; SU, sunflower; SO, sorghum; MA, maize; C₀, nugget; C, partial sill; C₀ + C, sill; a, range; C₀/(C₀ + C), nugget-to-sill ratio; SpD, spatial dependence; S, strong; M, moderate; R², coefficient of determinants, MAE, mean absolute error; RMSE, root mean square error; ASE, absolute standard error.

). The nugget-to-sill relationship indicates a spatial dependence that goes from moderate (three cases) to strong (14 cases). The range varies from a minimum of 56 m in SO_2010 to a maximum of 131 m in BW_2005. Prediction accuracy was assured by R² values that were generally high, apart from 10 July in SO_2010 (R² = 0.65); the error terms (MAE, RMSE; and ASE) were from low to very low. The good fitting of semi-variogram models, including the single case of modest R² values, can be visually checked by the interpolation of measured values (Figure S2 in the Supplementary Materials).

The maps of NDVI (actual values) outline a very differentiated behavior (Figure 5). Sometimes, the pale and dark green prevail in the map (e.g., 5 and 29 April 2007 ); some other times, orange and red prevail (e.g., 15 May 2007 and 24 July 2018); in the rest of the cases, a more balanced palette of colors is shown across the map.

3.4. Two-Zone Maps of Soil Properties, NDVI Data, Their Combination, and Agreement with Yield

Two-zone maps of the three soil properties and five NDVI, based on single property datasets, resulted in the best correlation with yield (Table 8 and Table 10, respectively), and are shown in Figure 6. Their composition in terms of data points, average value, and agreement with the two-zone multiple yield map are reported in

Table 12. Statistical differences between two-zone single soil properties and NDVI indices, and their pixel agreements (percentage) with two-zone multiple crops yield map.

Variables	Zones	Data Points	Mean ± SD	Agreement (%)
Clay (g/kg)	High clay Low clay	67	112 ±7.39 a	94
		32	79.3 ± 10.42 b
CaCO3 (g/kg)	High CaCO3	49	146 ± 1.31 a	76
	Low CaCO3	50	142 ± 1.05 b
Nitrogen (g/kg)	High Nitrogen	54	1.69 ± 0.08 a	81
	Low Nitrogen	45	1.48 ± 0.07 b
NDVI 15-Apr-2005	High NDVI	78	0.70 ± 0.012 a	95
	Low NDVI	21	0.65 ± 0.032 b
NDVI 29-Apr-2007	High NDVI	84	0.67 ± 0.006 a	89
	Low NDVI	15	0.64 ± 0.013 b
NDVI 21-Feb-2009	High NDVI	50	0.57 ± 0.014 a	77
	Low NDVI	49	0.52 ± 0.021 b
NDVI 24-Jun-2010	High NDVI	51	0.73 ± 0.004 a	78
	Low NDVI	48	0.71 ± 0.013 b
NDVI 24-Jul-2018	High NDVI	95	0.60 ± 0.033 a	78
	Low NDVI	4	0.38 ± 0.068 b

SD, standard deviation; means bearing different letters are significantly different within 2-zone classes at p < 0.05 according to the least-significant difference (LSD) test.

. Soil property maps evidence a variable pattern and a more balanced size between high and low zones; in contrast, NDVI maps show a more uniform pattern, associated with a prevalence of the high zone. The agreement with the multiple yield map goes from 76 (CaCO₃) to 95% (NDVI of 15 April in BW_2005). In each of the eight parameters, the high level is statistically differentiated from the low level. In addition, two-zone maps of soil traits and NDVI data exhibiting significant correlations with crop yields are also shown in the Supplementary Materials (Figures S3 and S4) with their statistical differences and pixel agreement with yields (Tables S3 and S4), respectively.

The two-zone soil and NDVI map, based on the combination of the eight above discussed parameters, and the two-zone yield map over seven years of multiple crops are shown in Figure 7. Their composition in terms of data points, average value, and agreement between pixels of the former and latter map are reported in

Table 13. Statistical differences between two-zone multiple soil properties and NDVI indices, and their pixel agreements (percentage) with two-zone multiple crops yield map.

Variables	Zones	Data Points	Mean ± SD	Agreement (%)
Clay (g/kg)	High values	89	103 ± 16.7 a	83%
	Low values	10	87.5 ± 19.8 b
CaCO3 (g/kg)	High values	89	144 ± 2.09 a
	Low values	10	143 ± 2.53 a
Nitrogen (g/kg)	High values	89	1.61 ± 0.13 a
	Low values	10	1.54 ± 0.14 a
NDVI 15-Apr-05	High values	89	0.70 ± 0.02 a
	Low values	10	0.63 ± 0.04 b
NDVI 29-Apr-07	High values	89	0.67 ± 0.01 a
	Low values	10	0.63 ± 0.01 b
NDVI 21-Feb-09	High values	89	0.55 ± 0.02 a
	Low values	10	0.50 ± 0.03 b
NDVI 24-Jun-10	High values	89	0.73 ± 0.01 a
	Low values	10	0.69 ± 0.01 b
NDVI 24-Jul-18	High values	89	0.60 ± 0.03 a
	Low values	10	0.47 ± 0.09 b
Multiple crops yield (%)	High	73	102 ± 1.72 a
	Low	26	94 ± 3.08 b

SD, standard deviation; means bearing different letters are significantly different within 2-zone classes at p < 0.05 according to the least-significant difference (LSD) test.

. The two maps are very similar, with a prevalence of the high zone (89 and 73 pixels in the two respective maps). Their agreement is high (83%). The high and low levels are statistically differentiated in six traits out of eight (CaCO₃ and pH being undifferentiated), as well as in the multiple crop yield.

4. Discussion

The seven annual crops surveyed in this work across the last two decades represent a good example of the georeferenced yield data that are increasingly being archived, as PA practices are extending to relatively small farmers. It does not often occur that the reliability and potential uses of such data are sufficiently understood by potential beneficiaries, who remain uncertain how to manage them. Therefore, experiments such as that carried out in this work, serve as good examples of the approach that may be followed in dealing with yield data, soil and vegetational properties, and the resulting outcomes.

The seven crops, intrinsically, do not have a high degree of similarity: three years of winter wheat; two years of maize for biomass production, one of which was cultivated as second crop in a season after a winter grass crop; the remaining two years with a summer oilseed (sunflower) and cereal (grain sorghum). They represent the variability that is normally met in mixed-cropping systems at intermediate latitudes, such as in the case of the investigated field. The missing years in the series are due to lack of equipment for georeferenced yield mapping, as in the case of the 2011–2017 gap, when a biogas plant built in the proximity commenced to be fed with biomass grown in this field, and a maize cutter-chopper provided with yield mapping was available only in 2018. However, such gaps may occur also in other cases; thus, the procedure we followed to gather the available years of data and use them all as a dataset, buffered by different ambient conditions and crop types, may be of general interest.

The mixed-cropping systems are seldom surveyed in the literature addressing the use of multiple-year yield data in spatial and temporal variability assessment, despite the fact that mixed-cropping systems are more common than mono-cropping systems, and they are supported by funding schemes in the European Union and elsewhere in the world. More often, multiple-year yield data of a specific crop are used in spatial and temporal variability assessment [4,40,41,42]. Alternatively, yield data from cropping systems with the same crop type as spring grain crops are also used [43]. When compared to these, the mixed-cropping systems involving different crop types in terms of harvested plant portions (grain vs. biomass/forage), environmental requirements (autumn- vs. spring-sown crops), or some other feature, represent a more challenging case. A recent work addressing such a case [44] envisaged the establishment of four classes of yield level and stability combined, using data from grain crops and pastures cultivated on the same fields in different years. The production stability index originating from this approach could discriminate field areas suited for different management, indicating, as in our case, the feasibility of mixed-cropping data fusion for interpreting field variability.

The seven crops staged a yield variation which was lower than observed in a previous study of ours’ [14]: average SD, which in standardized data (means = 100) corresponds to the coefficient of variation, was 11.4% across the seven years (Table 4), whereas it was 34% across five years in the cited work. This points out the strong diversity that may be observed between fields of similar size (~10 ha), shape, and topography (flat surface with underground pipe drainage) in the same area as the Po River floodplain in Northern Italy. This, in turn, was reflected in only two clusters as the optimum number in this work (Table 3), based on the principle of minimizing data dispersion. Therefore, no more than two different management zones would be needed to account for the existing variability, and this figure, based on seven-year average yield, was assumed as the reference number also for soil traits and remote NDVI data, in view of subsequent comparisons.

Despite the above-discussed modest variation, all the yield, soil, and NDVI datasets were fitted with good accuracy by semi-variogram models of various types: the stable model prevailed in yield fitting (Table 5 and Figure S1), while the exponential model prevailed in soil traits (Table 7 and Figure S5) and NDVI data (Table 11 and Figure S2). More to this, the fitted semi-variograms featured a generally strong spatial dependence, i.e., a nugget value that was modest, if not negligible, with respect to the total sill (Table 5, Table 7 and Table 11). This is a good premise for undertaking precision crop management, as the differences in, say, fertilizer dose that are associated with spatial distance, can be more reliably trusted [35]. For the same reason, a long range would be desirable, as it ensures that the control distance used for a site-specific operation, e.g., fertilizer spreading, falls in an area of relative uniformity, and the switch between different doses is sufficiently rapid to account for spatial variation in fertilizer needs. Based on the principle that the upper limit of cell size, i.e., the resolution, for a site-specific intervention should not exceed the length of half the range, to avoid excess heterogeneity within the cell [45,46], it results that the soil traits and NDVI data in our study have ranges long enough, in general, for a dose variation to be successfully implemented in real-time (Table 7 and Table 11). It is assumed that a vast majority of the ranges for soil and crop properties falls between 20 and 110 m [47]; therefore, comparing the range resulting from the semi-variogram analysis with the distance required by a specific implement to switch doses can be suggested as a preliminary check in view of variable crop input application.

In this respect, it is evidenced that the two-zone maps of soil properties, NDVI values, and multiple yields (Figure 6 and Figure 7) have pixels of sufficient size (30 m) for dose variations to be implemented when passing from the low to the high zone, and vice versa. This is associated with modest patchiness of the two-zone maps (Figure 6 and Figure 7): only some single NDVI maps (e.g., 24-06-2010 in Figure 6) denote a patchy distribution that could pose problems in the differential application of some factor doses. However, the multiple soil and NDVI map (Figure 7) was not affected by this drawback, as the eight single parameters contributing to the multiple soil and NDVI maps provided the buffering that is a desired characteristic in a map intended for future use. The multiple yield two-zone map (Figure 7) is an equally well-buffered map, as it originates from the multiple crops relative yield map (Figure 3b). The good agreement between the multiple soil and NDVI map, and multiple yield map (83% in Table 13), is a good premise for differential management of soil- and NDVI-based two zones in future cropping seasons.

5. Conclusions

Arable crop fields subjected for a number of years to georeferenced yield mapping obtain an amount of data deserving to be exploited in future crop management. The present work outlines a way to use these data in a field case study, by combining a temporal series of yield data from different crops into a standardized average yield, mapped through the field. This represents a robust indication of field spatial variability across varying weather and crop conditions. In parallel to this, spectral vegetation indices from remote sources and soil analyses could be obtained, and the soil and spectral vegetation data best correlated with yield could be combined into a multiple parameter mapped through the field. Geostatistical and clustering techniques were proven to be able to establish a suitable number of management zones across the field; in this case, it was only two zones, owing to the modest yield variability in the surveyed field. Lastly, a good agreement in the high and low zone between the multiple soil and spectral vegetation and multiple yield map, in this case 83 pixels out of 99 with agreement, is the ultimate proof that is needed to support the use of the soil- and vegetational-index-based map in the management of subsequent crops.