Workflow to Establish Time-Specific Zones in Precision Agriculture by Spatiotemporal Integration of Plant and Soil Sensing Data

Abstract

Management zones (MZs) are used in precision agriculture to diversify agronomic management across a field. According to current common practices, MZs are often spatially static: they are developed once and used thereafter. However, the soil–plant relationship often varies over time and space, decreasing the efficiency of static MZ designs. Therefore, we propose a novel workflow for time-specific MZ delineation based on integration of plant and soil sensing data. The workflow includes four steps: (1) geospatial sensor measurements are used to describe soil spatial variability and in-season plant growth status; (2) moving-window regression modelling is used to characterize the sub-field changes of the soil–plant relationship; (3) soil information and sub-field indicator(s) of the soil–plant relationship (i.e., the local regression slope coefficient[s]) are used to delineate time-specific MZs using fuzzy cluster analysis; and (4) MZ delineation is evaluated and interpreted. We illustrate the workflow with an idealized, yet realistic, example using synthetic data and with an experimental example from a 21-ha maize field in Italy using two years of maize growth, soil apparent electrical conductivity and normalized difference vegetation index (NDVI) data. In both examples, the MZs were characterized by unique combinations of soil properties and soil–plant relationships. The proposed approach provides an opportunity to address the spatiotemporal nature of changes in crop genetics × environment × management interactions.

1. Introduction

Crop yields and resource use efficiency (e.g., nutrients and water) have a strong spatial component, which can be observed over a wide range of scales, from regional to subfield [1]. At the field scale, yield variability in uniformly managed fields is often related to the spatial variability of soil properties and their impact on plant growth [2,3,4,5,6]. This variability can be addressed using precision agriculture practices [7], such as variable rate management (VRM) [8]. According to VRM principles, efficiency or crop production can be increased by varying agronomic inputs over a field according to varying soil and crop conditions [9]. Within a field, areas with similar soil properties (e.g., texture, topography, water holding capacity) that can be managed uniformly are commonly called management zones (MZs) [10]. To justify the subdivision of a field into MZs, there should be sizeable difference in soil properties between MZs [9,11]. The use of MZs for VRM has been shown to increase productivity, decrease costs, and/or reduce environmental impacts of agronomic practices [12,13]. Several authors have shown that within-MZ management should change over time [14,15,16,17,18,19,20], such an approach is often referred to as “dynamic VRM” [14].

Shanahan et al. [21], Long et al. [22], and Quebrajo et al. [23] indicated that coupling soil and plant data should increase the efficiency of MZ-based precision agriculture. Several sensor measurements can be used as indicators of soil spatial variability, including: apparent electrical conductivity, gamma-ray spectrometry, visible and near-infrared reflectance, and penetrometry [24,25,26]. Throughout the growing season, field-scale crop canopy measurements from near-ground (e.g., tractor mounted sensors) and remote sensing (e.g., unmanned aerial vehicles—UAV, satellites) can be used to provide information about plant health [27]. In particular, visible, near-infrared, and thermal sensing can be used to infer crop water status, nitrogen deficiency, and the effect of other biotic/abiotic stressors [28,29,30,31] and as an indicator of potential yield [16,32,33,34,35].

According to current practices, MZ-designs derived from soil maps, crop sensing, and/or historical yield maps are generally designed once and then implemented year after year. Nevertheless, several authors have reported inconsistent benefits for this precision management strategy [21,22]. Spatial patterns in yield tend to change from year-to-year, mostly because of changes in meteorological conditions [12,36,37,38]. In other words, the spatial patterns of most soil properties are fairly stable in time, but at different times plants may be limited in different ways at the same location because of the influence of transient factors affecting the soil–plant relationship, such as meteorological factors and agronomical management [12,38,39].

Sadler et al. [40] indicated that there is a need for accurate and inexpensive systems to delineate dynamic management zones, obtained by sensing within-field variability in real time, so that agricultural management can be controlled adaptively. Recent research strongly suggested that MZ designs should change over time, both intra- and inter-seasonally [19,20,38,41,42,43,44]. Myers [45] formally justified the need for time-specific and spatially dynamic VRM through the “fundamental theorem of precision agriculture production” where Crop Yield is a function of Genetics × Environment × Management × Space × Time interactions. According to this theorem, the spatiotemporal variability of crop performance should be addressed by adjusting the agronomic prescriptions over time and space. Several examples of protocols and data analysis workflow for static soil and/or plant-based MZ delineation are present in the scientific literature [46,47,48,49]. To our best knowledge, time-specific MZ delineation based on soil and in-season plant information has not been commonly discussed in the literature. Particularly, there is a lack of protocols and analytical workflow that farm managers, agricultural consultants, and scientists can use to take advantage of free/inexpensive in-season crop information (e.g., from UAV, Sentinel 2 satellite) and high-resolution soil maps.

We aim to present a novel workflow for the selection of time-specific MZs according to in-season spatial measurements of crop growth status and its relationship with high-resolution soil spatial information. The MZ should identify areas with homogeneous and unique (within a single field) soil–plant relationships. The MZ-delineation workflow will be described in detail. We will also provide two examples on how to implement the workflow. The first example is based on synthetic data. The second example uses data from a maize (Zea mays L.) field in northeastern Italy.

2. Materials and Methods

2.1. Time-Specific MZ-Delineation Workflow

The time-specific MZ delineation can be implemented as follows:

• STEP 1. Soil and time-specific plant spatial information acquisition, pre-processing, interpretation, and interpolation
• STEP 2. Time-specific sub-field soil–plant modeling
• STEP 3. Time-specific MZ delineation with cluster analysis
• STEP 4. Evaluation/interpretation of time-specific MZ design

2.1.1. Soil and Time-Specific Plant Spatial Information

In STEP 1, high-resolution spatial measurements for target soil properties and in-season plant-canopy information, such as canopy reflectance, are acquired, pre-processed, interpreted, and interpolated. In-season plant-canopy reflectance is acquired and used as an indication of crop status. Soil spatial information is used to interpret the crop canopy measurements.

Soil sensor data acquisition should be carried out according to established protocols [50,51,52] to increase the accuracy and consistency of the survey across large areas. Attention should be paid to selecting those sensors which represent the spatial variability of soil properties known or believed to influence a crop at the site of interest [53,54]. For information on sensor data pre-processing (e.g., conversion of spatial coordinates, removals of outliers), readers are referred to the first protocol step of Córdoba, Bruno, Costa, Peralta, and Balzarini [49]. Subsequent to the sensor surveys, soil sampling should be carried out across the root zone (e.g., 0–1.2 m) to calibrate/interpret the sensor readings [50]. Sensor-directed sampling schemes can be used to minimize the number of sampling sites [55,56]. Exploratory analyses, such as correlation analysis and principal component analysis (PCA) should be carried out to investigate the relationships between soil sensor and laboratory soil analyses. The strength of the relationships between collocated soil sensor values and laboratory soil analyses should be investigated. If these relationships are moderately to very strong, the sensor data can be considered as an indicator of spatial variability of the target soil properties. Then, sensor data can be used to generate maps of the selected soil properties [57,58]. Soil maps should be generated only if acceptable prediction errors [59] are obtained. Alternatively, for weak to moderately strong relationships, the soil sensor maps should be used as qualitative indicator of soil spatial variability. Córdoba, Bruno, Costa, Peralta, and Balzarini [49] describe how to practically process interpolated data to obtain a raster of desired block support (e.g., of the same resolution chosen for the MZ design) in the second step of their MZ-delineation protocol.

Free (e.g., Sentinel 2 satellite) or inexpensive (e.g., from UAV) crop measurements are available throughout the growing season with moderately-high temporal resolution. For example, the Sentinel 2 satellite from the European Space Agency provides multi-spectral canopy reflectance at the 10 × 10-m resolution, with a 5-day revisit time, free of charge. Remote sensing of crop canopy data can be used to calculate vegetation indices [34]. Point measurements, such as those from tractor-mounted active spectrometers [30], should be pre-processed and interpolated similarly to soil sensing measurements. High-resolution raster data, such as that from satellite imagery, may need to be re-gridded and scaled to the selected block support.

2.1.2. Time-Specific Sub-Field Soil–Plant Modeling

In STEP 2, soil information from STEP 1 is used to interpret in-season (i.e., time-specific) measurements of crop status. The interpretation of crop canopy sensing measurements is not straightforward—as they are influenced by species × growth-stage × stress levels × soil background interactions [60,61,62].

Moving window spatial regression modeling, such as geographically-weighted regression (GWR) [63,64], can be used to understand the local (i.e., sub-field scale) variation in the plant–soil relationship. With GWR, a regression is run for each grid location, rather than for the whole study area [65]. Soil map(s) are used as the independent (explanatory) variable(s) and the in-season crop status maps are the dependent variable. When multiple explanatory variables are available, one should consider standardizing them. The GWR allows non-stationarity of the regression equation parameters a (e.g., intercept, slope), estimating their values at each location i. For a dependent variable y the equation reads:

$$y_i=a_{0i}+a_{1i}{x}_{1i}+\dots +a_{ki}{x}_{ki}+{\epsilon }_i$$

(1)

where ε, is a random error term, a₀ is the regression intercept, and a₁ to a_k are the regression coefficient(s) for each explanatory variable(s) x₁ to x_k. The spatial variability of the soil–plant relationship can be described with the maps of local a₁ to a_k coefficients (i.e., the regression slope[s]). Maps of GWR slope coefficient(s) show the spatial variability of the impact (i.e., sensitivity) of the explanatory variable on the regression [65].

In the GWR framework, spatial weighting is determined by incorporating all the dependent and explanatory variables falling within a geographical kernel of each target feature [65]. The values of the regression parameters and goodness-of-fit of the GWR depend on how the kernel size is chosen [66]. Maps of the estimated dependent variable, local coefficient of determination (R²), and local Pearson correlation coefficient r can be generated with the GWR [65]. The GWR is available in commercial GIS software platforms (e.g., ArcMap’s Spatial Statistics package [65], version 10.5.1; ESRI, Redlands, CA, USA), and freeware (e.g., spgwr package in R, version 3.4.1; the R Foundation for Statistical Computing Platform, Vienna, Austria).

2.1.3. Time-Specific MZ Delineation with Cluster Analysis

In STEP 3, the MZ delineation is carried out using the soil map(s) and the time-specific GWR slope maps from STEP 2 as ancillary variables. As indicated by Córdoba, Bruno, Costa, Peralta, and Balzarini [49], fuzzy c-means (also known as “k-means”) unsupervised clustering algorithms [67] can be employed to classify the data into MZs. ArcMap’s Grouping Analysis tool (e.g., [68]), the Management Zone Analyst software [69] or the EZZone online tool [70] can be used to delineate MZs. Several MZ designs can be tested. The optimum number of MZs can be identified using cluster validity functions, including: the Calinski–Harabasz criterion [71], the fuzziness performance index [67], the normalized classification entropy index [67,72], and the Jenks optimization method [73]. The Calinski–Harabasz criterion (CHC), also known as pseudo F-statistic, describes the ratio between within-MZ similarity and between-MZ differences. It is defined as:

$$CHC=\frac{BMZSS/\left(MZn-1\right)}{WMZSS/\left(N-MZn\right)}$$

(2)

where N is the number of data points, MZn is the number of considered MZs, BMZSS is the between-MZ sum of squares, and WMZSS is the within group sum of squares. The larger the value of the CHC the higher are the within-MZ homogeneity and between-MZ differences. Finally, one may consider smoothing the fuzzy c-means clustering results to reduce zone fragmentation [49,70].

2.1.4. Time-Specific MZ-Design Quality Control and Interpretation

In STEP 4, the quality of the MZ design from STEP 3 should be checked. Each MZ should identify a unique combination of soil and soil–plant relationship characteristics. To infer differences across MZs, parametric analysis of variance or the Kruskal–Wallis (KW) rank test [74] can be used. The KW test is a nonparametric analysis assessing if samples originate from the same distribution. The test can be used as alternative to the standard analysis of variance when assumptions for parametric testing are not met.

2.2. Synthetic Data Example

The proposed workflow can be demonstrated using synthetic data for soil clay content (percentage) and normalized difference vegetation index (NDVI) [75].

A spatial random field for clay content was generated using the RFsimulate function in the RandomFields package in R. The function estimated a Gaussian random field (RPgauss function) having an exponential covariance function (RMexp function with variance = 25 and scale = 300), a nugget effect (RMnugget function with variance = 5), and a pure trend model with covariance 0 (RMtrend function with mean = 15). The simulations were carried out over a square of size 250 × 250 cells. The random field was then converted into a 250 × 250-pixel raster with pixel size = 1 m using the function raster from the raster package in R. The clay content raster was then exported from R to a text file using the writeRaster function. We refer to this raster as the true clay content (TCC) map.

Next, we mimicked a field procedure for generating a soil clay content map, using TCC as the true, unknown, high resolution clay content. Accurate approximations of true soil properties can be obtained when high-resolution covariates, such as those from proximal-soil sensing [76], are available [77]. Soil sensor surveys can be calibrated to estimate the spatial variability of a target soil property by using laboratory measurements on collocated soil cores [58]. Heggemann et al. [78] reported that gamma-ray sensor readings can be calibrated to predict texture values with high accuracy (up to ~95 percent of observed variance in soil texture). We mimicked a typical sensor survey (e.g., Lesch [55]) in which 1870 data points (average nearest neighbor distance = 3.2 m) were spread across 24 nearly-parallel transects (average nearest neighbor distance = 9.9 m). Values from the TCC map were extracted at the sensor survey locations. A random error having mean = 0 and variance equal to 5 percent of the extracted TCC values was added to the sensor data. This simulated a realistic sensor calibration with goodness-of-fit with a R² close to 0.95. The spatial autocorrelation of the calibrated sensor measurements was described with a spherical semivariogram having range = 50.3 m, and nugget equal to 67% of the total sill. This spatial structure was similar to those reported by other authors in a 5-ha clay loam field in Italy [79] and in a cluster of fields with contrasting soil properties (from clay to gravelly) in Germany [80]. These calibrated sensor measurements were then interpolated using simple kriging in with ArcMap’s Geostatistical Analyst package. The resulting sensor-derived clay content (SCC) map was retained for further analyses.

NDVI is a vegetation index calculated from surface reflectance in the red (RED) and near-infrared (NIR) regions of the electromagnetic spectrum. It is defined as:

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

(3)

NDVI ranges from −1 to +1. NDVI for agricultural crops usually ranges from ~0.1 to ~0.9, with lush vegetation generally having high NDVI [81,82]. The NDVI information was simulated with the understanding that, often, the spatial relationship between remote sensing canopy measurement and collocated soil properties has both deterministic and spatial random components [83,84,85]. The deterministic component of the relationship can be a linear model between the soil property and the available remote sensing plant information [84]. The spatial random component is often equal to the field of spatially correlated residuals from the deterministic linear model [57,83]. We simulated the NDVI as follow:

$$\mathrm{NDVI}=S\times \left(TCC+SpERR\right)+O$$

(4)

where SpERR was a spatial error raster, S was a scaling factor = 0.01, and O was an offset coefficient = 0.4. SpERR was generated in R, using the RandomFields package. We estimated a Gaussian random field based on a model with exponential covariance function (variance = 7.5 and scale = 50) and a nugget effect (variance = 5).

A GWR with 30-m bandwidth size was used to model the spatial variability of the synthetic NDVI image using the SCC map as explanatory variable. Then, MZs were delineated with a fuzzy c-mean unsupervised cluster analysis. Clustering was carried out in ArcMap using the Grouping Analysis tool from the Spatial Statistics package. The SCC and the time-specific GWR slope maps were used for the cluster analysis. No spatial constraints were set for the cluster analysis: the covariate data points could be grouped according to their value, disregarding the values of their geographical coordinates (i.e., points did not need to be neighbors to be part of the same MZ). The classifications were carried out for scenarios outputting two, three, and four MZs. Then, the optimal number of MZs was identified using the CHC. Differences in clay (from the SCC map), NDVI, and GWR slope across MZs were investigated with the KW test using STATISTICA (version 12, StatSoft Inc., Tulsa, OK, USA).

2.3. Field Data Example

2.3.1. Site Description

Another example is given using experimental data from researchers at the University of Padua, Italy [86,87,88]. Here, we describe two time-specific MZ designs for a maize field based on remote sensing imagery acquired during the late kernel blister stage (R-2), in two consecutive growing seasons.

The study site (Figure 1) is a 21-ha field located at Chioggia, Venice, Italy, along the southern margin of the Venice Lagoon (45°10′7.0′′ N, 12°13′55.0′′ E). Previous studies [86,87,88,89] discussed the relationship between soil and maize yield at this site. The site lies below average sea level, and the groundwater level is kept fairly shallow (approximately between −0.6 to −1.5 m below ground level) by a pumping station [87] to promote sub-irrigation [88]. The soil in the area is classified as Molli-Gleyic Cambisol [90] with two coarsely-textured paleochannels crossing it with SW-NE direction (Figure 1). The field is affected by soil salinity due to its proximity to the Venice Lagoon [86].

Here, we analyzed the two maize growing seasons (April to September), 2010 and 2011. The two growing seasons differed greatly in terms of meteorology. Compared to the average April to September rainfall from 1993 to 2012, 2010 was rainy (in the upper third quartile, 540 mm) and 2011 was a drought year (in the first quartile, 200 mm) [88]. The daily average reference evapotranspiration was 4.01 mm day⁻¹ in 2010 and 4.42 mm day⁻¹ in 2011. For the entire growing season, the reference evapotranspiration was 569 mm in 2010 and 672 mm in 2011 [88]. Agronomic management (base-dressing of 64 kg N ha⁻¹, and 94 kg P₂O₅ ha⁻¹ and urea top-dressing of 184 kg N ha⁻¹) and maize hybrid (PR32P26, Pioneer Hi-Bred Italia, Gadesco-Pieve Delmona, Italy) were the same in the two years [86].

2.3.2. Soil Maps and Plant Spatiotemporal Information

Geospatial measurements of soil apparent electrical conductivity (EC_a) were used as an indicator of soil spatial variability at the study site. EC_a over the 0–1.5 m soil profile (EC_a Deep) was measured in April of 2011 across the site with a frequency-domain electromagnetic induction sensor (CMD-1, GF Instruments, Brno, Czech Republic) at 20,470 geo-referenced locations. Soil samples were taken at 91 locations in May 2010 (Figure 1) down to 1.2 m with 0.3-m increments. Here we discuss salinity and texture for the 0–1.2-m soil profile. Texture was measured with a Mastersizer 2000 (Malvern Instruments Ltd., Great Malvern, UK) and salinity was measured as EC_1:2 (i.e., the electrical conductivity of a 1:2 soil–water extract) [91]. Principal component analysis (PCA) carried out with STATISTICA was used to visually assess similarities and differences between the EC_a Deep dataset and the selected soil properties [92]. The point measurements were preprocessed with a procedure comparable to that of Córdoba, Bruno, Costa, Peralta, and Balzarini [49] and spatially interpolated with ordinary kriging following Scudiero, Teatini, Corwin, Deiana, Berti, and Morari [86] on a 10 × 10-m block support—the desired support for MZ at the site.

Two remotely sensed measures of crop status were available. On 31 July 31 2010 (10:39:40 UTC) and 9 July 2011 (10:47:30 UTC) WorldView-2 (DigitalGlobe, Westminster, CO, USA) satellite scenes were acquired over the study area. For both years, acquisition dates corresponded to late R-2 [88]. Pre-processing procedures, including sensor calibration, atmospheric correction, and radiometric normalization, were applied according to Vicente-Serrano et al. [93]. Radiometric calibration was required to convert digital numbers to top-of-atmosphere radiances [W m⁻² sr⁻¹ μm⁻¹], using the absolute radiometric calibration factors and effective bandwidths for each band, according to the satellite data provider [94]. No topographic correction was applied to the images because the study area is relatively flat and the solar incident angles (28.11° in 2010; 23.58° in 2011) were quite similar at the two acquisition times. The top-of-atmosphere radiance was then transformed to surface reflectance through the 6S code [95]. The values of aerosol optical thickness (AOT) at 550 nm collected from a nearby AERONET (Aerosol Robotic Network, https://aeronet.gsfc.nasa.gov/) [96] station (45°18′50.0″ N, 12°30′29.9″ E) were used as input for 6S. The AOT measurements were obtained simultaneously to the satellite overpasses. WorldView-2 reflectance has 2 × 2-m spatial resolution over eight spectral bands, at wavelengths spanning from 400 to 1040 nm. The red band (RED, 630–690 nm) and the near-infrared band at 770–895 nm where used to calculate NDVI according to Equation (3). The NDVI maps were aggregated (i.e., the coarsened pixel is the average of the pixels within the aggregated cell) to the 10 × 10-m cell size, as suggested by Córdoba, Bruno, Costa, Peralta, and Balzarini [49].

2.3.3. Time-Specific Spatial Soil–Plant Modeling

The time-specific relationship between soil properties and in-season NDVI were described using geographically weighted regressions (GWRs). The GWRs were carried out using the NDVI maps as the dependent variable and spatial soil information (i.e., EC_a Deep) as the independent variable. GWRs were carried out with ArcMap using an adaptive kernel (i.e., moving window) of 70 neighbors. For locations with significant GWR, the observed-estimated NDVI relationship was used as indication of goodness-of-fit of the GWR models. The GWR slope map was selected as indicator of soil–plant relationship type.

2.3.4. Delineating MZs with Cluster Analysis

Clustering was carried out using the Grouping Analysis tool in ArcMap. Spatial soil information (i.e., EC_a Deep) and the time-specific GWR slope map were used for the cluster analysis. No spatial constraints were set for the cluster analysis. The classifications were carried out for scenarios outputting three, four, and five MZs. Then, the optimal number of MZs was identified using the CHC. To compare the proposed time-specific soil–plant based MZ delineation with a more traditional, static MZ design, EC_a Deep alone was used to select a static MZ design.

2.3.5. Evaluation and Interpretation of MZ Design

Differences across MZs in the two years were investigated using the KW test in STATISTICA. The following variables were tested for differences across MZs: EC_a Deep, GWR slope, GWR local r, EC_1:2, and sand and clay contents. Additionally, in-season potential yield estimations were used to evaluate the MZ designs. Appendix A reports on how yield maps (obtained from Scudiero, Teatini, Corwin, Dal Ferro, Simonetti, and Morari [88]) and NDVI measurements were used to calculate yield prediction maps.

3. Results and Discussion

3.1. Synthetic Data Example

Every time the crop status is mapped during the season (e.g., every week), a new set of MZ is delineated (i.e., time-specific MZ design). Clearly, this applies only for crop growth stages when crop canopy plays a relevant role in determining reflectance readings. For soil tillage, sowing, and early vegetation stages when surface reflectance is mostly determined by soil, a static MZ delineation approach based on the spatial variability of soil properties (e.g., [9]) may be more adequate.

Figure 2 outlines the proposed analytical protocol to delineate time-specific MZs from soil information and in-season crop information using the synthetic dataset. In Step 1 of the protocol, the soil and plant information were acquired and pre-processed. The clay content (percentage) map had mean = 18.6%, standard deviation = 4.5%, minimum = 9%, and maximum = 30.1%. The NDVI map had mean = 0.58, standard deviation = 0.04, minimum = 0.45, and maximum = 0.69. The NDVI values are typical of vegetative stages for crops such as maize [97] and soybean (Glycine max (L.) Merr.) [98]. A simple linear model between the two maps had R² = 0.61, with an intercept of 0.44 and a positive slope of 0.007. Similar goodness-of-fit values between soil properties and vegetation reflectance were reported by Gomez et al. [99] for clay and by Gomez et al. [100] for soil organic carbon. In Step 2, the clay map was used as the explanatory variable in a GWR with NDVI as dependent variable. Through the GWR analysis, clay explained 94.5% of the observed variance of NDVI. Such high goodness-of-fit can be found in real-world data: Scudiero, Corwin, Wienhold, Bosley, Shanahan, and Johnson [54] observed R² ranging between 0.83 and 0.94 for their GWR analyses between soil maps and remotely sensed winter wheat (Triticum aestivum L.) canopy reflectance. The GWR slope had mean = 0.009, standard deviation = 0.008, minimum = −0.018, and maximum = 0.038. In Step 3, the clay and the GWR slope maps were used to delineate MZs via unsupervised fuzzy c-means clustering. The CHC indicated that the best number of MZs was three.

In Step 4, the MZ-design was evaluated and interpreted. The three MZs identified unique combinations of clay content and GWR slope values. MZ 1 was characterized by low NDVI, low clay content, and large GWR slope values (indicating high sensitivity of NDVI to clay). MZ 2 was characterized by high NDVI, high clay, and moderately low GWR slope values. MZ 3 was characterized by low NDVI, low clay content, and the smallest GWR slope values.

3.2. Field Data Example

3.2.1. Soil and Plant Information

Figure 3a shows an example of short-scale soil spatial variation at the site. Figure 3b depicts the frequency distribution of measured soil salinity. At this site, soils with EC_1:2 >0.79 dS m⁻¹ should be considered salt affected, with EC_1:2 >1.4 dS m⁻¹ indicating strongly saline soils [88]. Figure 3c reports the observed texture at the 91 soil sampling locations. Soil texture is predominantly in the loam and sandy loam classes. Few samples are in the clay loam and loamy sand class.

The EC_a Deep map is shown in Figure 4a. Descriptive statistics for the map are reported in

Table 1. Descriptive statistics for the soil apparent electrical conductivity map (EC_a Deep), and the 2010 and 2011 Normalized Difference Vegetation Index (NDVI) maps.

Variable	Mean	Median	Minimum	Maximum	Standard Deviation
ECa Deep (dS m−1)	0.96	0.86	0.32	2.49	0.36
2010 NDVI	0.790	0.80	0.470	0.867	0.046
2011 NDVI	0.829	0.83	0.664	0.858	0.018

. EC_a Deep showed a Pearson correlation coefficient of 0.89 with salinity, −0.51 with sand, and 0.39 with clay. The PCA extracted two factors (with eigenvalues of 2.61 and 0.89) explaining 65.3% (first component) and 22.3% (second component) of the total variance in the dataset. The bi-plot in Figure 3d—see Abdi and Williams [92] for interpretation guidelines—indicates that the first component contrasts the positive contribution of clay content and salinity on the EC_a Deep measurements with the negative correlation between sand content and EC_a Deep. In Figure 3d, clay content was not clustered with salinity and EC_a Deep, indicating that the three variables were not collinear, yet negatively correlated with sand content. According to the correlation and PCA analyses, high values of EC_a Deep were mainly interpreted as an indication of high salinity and fine soil texture. EC_a readings larger than 1–2 dS m⁻¹ are most likely due to high soil salinity rather than other edaphic factors contributing to soil conductivity [101,102]. Conversely, lower EC_a Deep was interpreted as an indication of coarser soils with low salinity. Note that Scudiero, Teatini, Corwin, Deiana, Berti, and Morari [86] reported that the spatial patterns of EC_a Deep remained stable over time at the study area.

The two in-season NDVI maps for the late R-2 growth stage are reported in Figure 4b (2010) and Figure 4e (2011). Table 1 reports the descriptive statistics for the two NDVI maps. The 2010 NDVI was more heterogeneous (mean = 0.79, standard deviation = 0.046) than that of 2011 (mean = 0.83, standard deviation = 0.018). During the R-2 growth stage, environmental stressors (e.g., drought, heat, severe nutrient deficiency) may prevent kernels from developing properly (e.g., kernel abortion at the ear tip) [103,104].

3.2.2. Plant–Soil Relationship at Different Time Points

Figure 4 summarizes the GWR analysis between soil EC_a Deep and in-season NDVI for 2010 and 2011. For 2010, 1149 cells (62.5% of total) were characterized by a significant (p < 0.05) local Pearson correlation coefficient r. For 2011, 56.8% of the cells had a significant local r. At the locations with significant local r, the observed-estimated NDVI relationship was characterized by R² = 0.73 in 2010 (Figure 3d) and R² = 0.62 in 2011 (Figure 3g). There was a strong negative relationship in both years (r = −0.36 in 2010 and r = −0.43 in 2011) between the significant local r and EC_a Deep maps. The slope of this relationship was steeper in 2011 (−0.46 with standard error = 0.027) than in 2010 (−0.38 with standard error = 0.026). This may indicate that low EC_a Deep (e.g., indicating coarse texture) was a greater constrain in 2011 than 2010. Scudiero, Teatini, Corwin, Dal Ferro, Simonetti, and Morari [88] showed that water stress in areas with coarser soil texture is particularly limiting in dry years, such as 2011.

In areas where the significant local r values were consistently negative in both years, the EC_a Deep averaged 1.17 dS m⁻¹ (standard deviation = 0.39 dS m⁻¹), whereas, in areas with negative local r only in 2010, it averaged 0.94 dS m⁻¹ (standard deviation = 0.28 dS m⁻¹). This suggested that only very high EC_a Deep, which indicates high salinity, consistently limited crop growth in the two years. As discussed by other authors [105,106,107], the effects of high salinity on crops are generally stable throughout different growing seasons.

The GWR slope maps (Figure 4c for 2010, Figure 4f for 2011) changed remarkably between the two years, indicating a change in the sensitivity of plant NDVI to changes in soil properties. Changes of crop growth and yield spatial patterns are expected between seasons having widely different meteorology, as discussed by Maestrini and Basso [38] and McBratney, Whelan, Ancev, and Bouma [37].

3.2.3. Time-Specific MZ Delineation

The GWR slope and EC_a Deep maps were used to drive the delineation of in-season site-specific MZs. The CHC index indicated that a five-MZ configuration was optimal for 2010. The CHC was 1901.3 for three MZs, 1979.7 for four MZs, and 2002.7 for five-MZs. In 2011, the CHC index indicated that a four-MZ configuration was optimal. CHC was 1248.9 for three MZs, 1383.4 for four MZs, and 1374.0 for five MZs. The two MZ delineations are shown in Figure 5a (for 2010) and Figure 5b (for 2011).

Overall, 45.1% of the cells were classified in the same MZ at both times. The consistency rate was 94.3% for MZ I, 75.6% for MZ II, 20.1% for MZ III, and 48.5% for MZ IV. Conversely, changes in local GWR slope led to 54.9% of the cells being assigned to a different time-specific MZ. As reported by Maestrini and Basso [38], it is reasonable to expect portions of a field to have stable (e.g., consistently high yields) and unstable crop outputs. Management in the unstable areas should be addressed within each growing season to meet desired agronomic output goals [18,38]. The proposed methodology helps characterizing time-specific changes in the crop output according to soil spatial variability.

3.2.4. Time-Specific MZ-Design Quality Control and Interpretation

In both years, each time-specific MZ was characterized by unique combinations of the variables used in MZ delineation. Moreover, the time-specific MZs differed greatly in terms of EC_a Deep, significant GWR r and slope, salinity, and texture (

Table 2. Count (n), mean, and standard deviation (Std. Dev.) of apparent electrical conductivity for the 0–1.5-m soil profile (EC_a Deep), local regression slope and Pearson correlation coefficient (r) obtained from significant geographically weighted regressions (GWR), soil salinity (EC_1:2), sand and clay contents, and potential yield estimations at each time-specific management zone (MZ) for 2010 and 2011 (n.a. = not available; NA = not assessed).

Variable (Unit)	MZ	2010				2011
		n	Mean	Std. Dev.	KW 1	n	Mean	Std. Dev.	KW 1
ECa Deep (dS m−1)	I	53	0.44	0.10	e	263	0.65	0.21	d
	II	779	0.68	0.12	d	1164	0.78	0.15	c
	III	402	1.54	0.19	a	146	1.24	0.45	a
	IV	513	1.09	0.12	b	574	1.38	0.23	b
	V	400	0.80	0.12	c	n.a.	n.a.	n.a.	n.a.
GWR slope	I	53 (100.0% 2)	0.790	0.316	a	221 (84.0% 2)	0.073	0.040	a
	II	371 (47.6% 2)	0.065	0.113	b	488 (41.9% 2)	−0.024	0.028	b
	III	334 (83.1% 2)	−0.069	0.069	c	146 (100.0% 2)	−0.132	0.050	c
	IV	294 (57.3% 2)	−0.048	0.055	c	364 (63.4% 2)	−0.018	0.029	b
	V	397 (99.3% 2)	−0.176	0.092	d	n.a.	n.a.	n.a.	n.a.
GWR local r	I	53 (100.0% 2)	0.587	0.102	a	221 (84.0% 2)	0.433	0.126	a
	II	371 (47.6% 2)	0.028	0.411	b	488 (41.9% 2)	−0.279	0.285	b
	III	334 (83.1% 2)	−0.439	0.239	c	146 (100.0% 2)	−0.597	0.190	d
	IV	294 (57.3% 2)	−0.399	0.281	c	364 (63.4% 2)	−0.308	0.349	c
	V	397 (99.3% 2)	−0.537	0.150	d	n.a.	n.a.	n.a.	n.a.
EC1:2 (dS m−1)	I	1	0.21	NA	ab	11	1.73	1.33	b
	II	32	1.38	0.89	b	49	1.49	1.01	b
	III	17	3.67	1.46	a	3	3.16	1.32	ab
	IV	26	2.28	1.18	a	28	3.03	1.45	a
	V	15	1.35	0.85	b	n.a.	n.a.	n.a.	n.a.
Sand (%)	I	1	71.1	NA	ab	11	54.1	15.6	ab
	II	32	50.6	12.6	b	49	53.1	12.4	a
	III	17	42.8	10.9	b	3	37.5	13.6	ab
	IV	26	42.7	8.8	b	28	42.1	9.2	b
	V	15	63.8	8.5	a	n.a.	n.a.	n.a.	n.a.
Clay (%)	I	1	9.0	NA	ab	11	12.6	6.6	ab
	II	32	13.5	5.1	a	49	12.6	4.8	b
	III	17	17.4	5.7	a	3	20.8	9.3	ab
	IV	26	16.8	4.7	a	28	17.3	4.7	a
	V	15	8.7	3.2	b	n.a.	n.a.	n.a.	n.a.
Potential Yield (Mg ha−1)	I	53	7.03	1.90	bc	263	10.0	1.7	b
	II	779	8.38	2.24	a	1164	10.5	1.2	a
	III	402	7.61	1.75	b	146	9.3	2.1	c
	IV	513	7.19	2.01	c	574	10.3	1.3	b
	IV	400	7.08	1.90	c	n.a.	n.a.	n.a.	n.a.

¹ Different letters are significantly different between MZs at the p < 0.05 level according to the Kruskal–Wallis (KW) test. ² Percentage of cells having significant GWR within the MZs.

). Note that the local GWR r and slope, which can be used to interpret the EC_a Deep–NDVI relationship, remained fairly consistent over time at the different MZs.

In both years, MZ I consisted of soils with low EC_a Deep, low salinity, and coarse soil texture. MZ I had the highest positive GWR slope values, meaning NDVI would decrease with increasing sand content. Scudiero, Teatini, Corwin, Dal Ferro, Simonetti, and Morari [88] monitored crop water and salt stress at five soil–plant–water monitoring stations at the site. One of their stations (“Station E”) was located in our MZ I. Scudiero, Teatini, Corwin, Dal Ferro, Simonetti, and Morari [88] indicated that at Station E maize was not stressed by salinity but was under water stress, particularly in 2011, when water stress was described as “severe”.

In contrast, MZ II in both years had GWR slope and r distributions overlapping with 0, indicating little to no influence of soil on crop NDVI variability. This is perhaps the reason why MZ II was consistently characterized by the highest potential yield estimations (Appendix A). The very saline MZ III (highest measured average EC_a Deep and EC_1:2) was characterized by lower potential yield estimations for both years in comparison with the other MZs (Table 2), as expected for (unmanaged) very saline portions of farmlands [106]. Similarly, the moderately saline MZ IV (second highest average EC_a Deep), was characterized by the second lowest yield predictions in 2010 and 2011, together with MZ I. MZ V was only selected in 2010. It grouped soils with moderately high EC_a Deep (third highest MZ average) together with the strongest negative GWR slope values and low potential yield estimations. In 2011, 89.5% of the 2010 MZ V locations were classified into MZ II. Previous research focusing on spatiotemporal variability of yield patterns [38], indicates that fields generally have areas of stable and unstable yield patterns. For unstable areas, in-season NDVI is suggested as the best predictor for yield spatial distribution [38]. Our proposed approach further refines the NDVI information by interpreting its spatial variability as a function of soil spatial variability.

A static MZ design based on EC_a Deep only (Figure 5c) would identify four areas with significantly different soil properties (

Table 3. Count (n), mean, and standard deviation (Std. Dev.) of apparent electrical conductivity for the 0–1.5-m soil profile (EC_a Deep), soil salinity (EC_1:2), and sand and clay contents for the static management zones (MZs) delineated using EC_a Deep only.

Variable (Unit)	Static MZ	n	Mean	Std. Dev.	KW 1
ECa Deep (dS m−1)	1	718	0.92	0.10	c
	2	181	1.71	0.17	a
	3	444	1.31	0.12	b
	4	804	0.63	0.10	d
EC1:2 (dS m−1)	1	28	1.98	1.04	b
	2	9	4.07	1.51	ab
	3	21	2.82	1.25	a
	4	33	1.07	0.68	c
Sand (%)	1	28	50.3	11.9	ab
	2	9	41.0	9.3	b
	3	21	41.1	9.8	b
	4	33	55.9	12.9	a
Clay (%)	1	28	13.3	5.0	bc
	2	9	18.3	5.7	ab
	3	21	17.9	5.0	a
	4	33	11.8	4.9	c

¹ Different letters are significantly different between MZs at the p < 0.05 level according to the Kruskal–Wallis (KW) test.

), but very heterogeneous and inconsistent in terms of GWR local r and slope (

Table 4. Mean, standard deviation (Std. Dev.), and count (n) for the correlation coefficient (r) and slope obtained from significant geographically weighted regressions (GWRs) for the static management zones (MZs) delineated using soil apparent electrical conductivity for the 0–1.5-m soil profile (EC_a Deep) only.

GWR Variable	Static MZ	2010			2011
		n	Mean	Std. Dev.	n	Mean	Std. Dev.
Local r	1	490 (68.2% 1)	−0.374	0.351	349 (48.6% 1)	−0.156	0.383
	2	163 (90.1% 1)	−0.474	0.158	148 (81.8% 1)	−0.542	0.158
	3	313 (70.5% 1)	−0.408	0.281	289 (65.1% 1)	−0.303	0.408
	4	483 (60.1% 1)	−0.098	0.472	433 (53.9% 1)	−0.041	0.420
Slope	1	490 (68.2% 1)	−0.086	0.148	349 (48.6% 1)	−0.009	0.039
	2	163 (90.1% 1)	−0.077	0.065	148 (81.8% 1)	−0.074	0.068
	3	313 (70.5% 1)	−0.058	0.065	289 (65.1% 1)	−0.023	0.045
	4	483 (60.1% 1)	0.066	0.304	433 (53.9% 1)	−0.002	0.079

¹ Percentage of cells having significant GWR within the MZ.

).

4. Conclusions

Characterizing the temporal variability of the soil–plant relationship within the growing season and over years is very relevant to developing best crop management practices with VRM. Soil influences plants according to complex interactions between factors that may change in time, such as meteorological conditions, nutrient availability, and water content. Static MZs are not ideal when spatial patterns of soil–plant relationship change in time—because of changing meteorology and/or other transient factors. We presented a time-specific MZ-delineation workflow to address such spatiotemporal variability.

The proposed approach for time-specific MZs should be further tested and evaluated in the field. Future work should focus on multi-field and multi-year comparisons of time-specific and static MZs and their relative impact on crop yields and resource use efficiency (e.g., water, nutrients).