Geostatistical Analyses of Field Spatial Variability of Cotton Yield ()
1. Introduction
Agricultural field experiments are traditionally conducted in small plots with treatment replicates. However, large strip plots have been increasingly utilized in field experiments during the recent years because farmers’ acceptance of research results from small plot experiments has dropped. A common belief among farmers is that small plots are inferior to large strip plots because large strip plots are more similar to the farmers’ field conditions.
In geostatistical analyses, autocorrelation has often been used to describe the spatial variabilities of soil properties and plant characteristics and the degree of dependencies among neighboring observations in a field experiment. This information is utilized to figure out the adequate sampling interval for which observations remain spatially correlated and for designing sampling protocols [1] [2] [3] . The Moran’s I statistic and scatter plot is used to express the degree of global spatial autocorrelation that presents among the plots of a field experiment [4] [5] [6] . This statistic calculates an autocorrelation index called Moran’s I among groups of paired plots separated by increasing distances. This index is calculated by dividing the spatial covariation by the total variation in the data. The resulting Moran’s I values are in the range from −1 to 1. Positive Moran’s I values represent positive spatial autocorrelation, which indicates that similar values (either high or low) are spatially clustered. Negative Moran’s I values mean negative spatial autocorrelation, which indicates that neighboring values are not similar. Moran’s I value of 0 indicates no spatial autocorrelation or spatial randomness. On the other hand, the Localized Indicators of Spatial Autocorrelation (LISA) cluster map is frequently used as a tool for estimating local spatial autocorrelation [7] . It contains information on only those locations that have significant spatial autocorrelation.
Spatial variabilities of soil properties and plant characteristics are central in precision agricultural research.
2. Materials and Methods
2.1. Site Description and Experimental Design
A strip plot cotton experiment was carried out on a Lexington silt loam soil on a private farm near Brazil in Gibson County, TN from 2009 to 2011. The plots used in 2009 were repeatedly utilized for 2010 and 2011. Cotton was grown on the field before the initiation of this experiment. Five side-dress N treatments of 0, 45, 90, 134, and 179 kg N ha−1 (0, 40, 80, 120, and 160 lb N a−1) were examined as urea and ammonium nitrate (UAN, 32% N) in strip plots under a RCB design with three replications. All the plots received 45 kg N ha−1 in the form of chicken litter as pre-plant N before cotton planting by broadcasting on the soil surface each year. The strip plot was 11.6 m wide by 244.0 m long, and split into eight sub-plots of 30.5-m length. Cotton was seeded in a 97-cm row spacing on May 8, 2009, May 14, 2010, and May 21, 2011, and managed with the recommended management practices except the N treatments for cotton in the region [8] . The specific dates of side dress N treatment implementation were presented in the article of Yin (2016) [9] . The coordinates for the four experiment corners were collected using a GPS hand held receiver on August 12, 2009.
2.2. Sampling and Measurements
The following sampling and measurements were taken from each individual sub-plot each season. A composite soil sample was collected at a 60-cm depth with a Concord hydraulic soil probe before initiation of the side-dress N treatments but after the pre-plant N application in early 2009. A composite leaf sample (10 blades + 10 petioles) was collected at the early square and early, mid, and late bloom growth stages each year (except no measurement at early square in 2009). Canopy NDVI readings were recorded three to four times each year at about the same dates when leaf samples were taken with the GreenSeeker® RT 200 Data Collection and Mapping System (NTech Industries, Inc., CA). The methods and dates for the above soil and plant sampling and measurements and the results of these measurements are shown in Yin (2016) [9] .
Harvest aids were applied to terminate the crop at approximately 10 to 20 days before cotton harvest. Cotton was harvested from the six center rows of each sub-plot with the farmer’s cotton picker on November 6, 2009, September 30, 2010, and October 1, 2011. Seed cotton was weighed, and a subsample was collected for ginning. Each subsample was ginned on a 10-saw laboratory gin to determine gin turnout and obtain lint yield for each sub-plot.
2.3. Statistical Analysis
Pearson product-moment correlation coefficients were calculated among lint field, company NDVI, and leaf N at each sampling date of each year with PROV CORR in SAS (SAS Institute, Cary, NC). The GPS coordinates at the four experiment corners were imported into Arc View GIS (Environmental Systems Research Institute, Redlands, CA). Spatial variations of lint yield within the strip plot experiment were visualized in GIS maps with Arc View v.9.3 each year. A quadratic regression of lint yield was conducted using the classic and spatial error models in GeoDa 0.9.5-i (Beta) (University of Chicago, Chicago, IL) with a weight matrix created using a 2nd order queen’s contiguity model that includes all lower contiguity orders. In order to evaluate the spatial dependence of lint yield relating to the characteristics of the test field (not to the N treatments), the effects of side dress N treatments on lint yield were removed from the lint yield data using the spatial error model, and then the residual lint yield data were used to make Moran’s I statistic and scatter plot and the LISA cluster map. Moran’s I statistic and scatter plot and the LISA cluster map of residual lint yield were created in GeoDa using the 2nd order queen’s contiguity model that includes all lower contiguity orders. Probability levels < 0.05 were designated as significant for all analysis.
3. Results and Discussion
3.1. Correlations of Lint Yield and Leaf N with Canopy NDVI
In 2009, the correlation of lint yield with NDVI was statistically significant or close to significance at early, mid, and late bloom stages, and became stronger as the season progressed from early to late bloom (Table 1). The correlation of lint yield with leaf N
![]()
Table 1. Correlations among lint yield, canopy NDVI, and leaf N in the N experiment at Gibson in 2009-2011.
was significant at early, mid, and late bloom stages, and became weaker as the season moved forward from early and mid bloom to late bloom (Table 1). Although the correlation of leaf N with NDVI was statistically significant at all growth stages, the determination coefficient (R2) was low; which suggests that NDVI may not be a strong indicator of plant N nutrition status during early to late bloom after a common rate of 45 kg N ha−1 from chicken litter was applied as pre-plant N before cotton planting.
The relationship of lint yield with NDVI was significant at early, mid, and late bloom stages in 2010 (Table 1). The association of lint yield with leaf N was significant at early square and mid and late bloom stages (Table 1). However, there was no significant correlation of leaf N with NDVI regardless of growth stage (Table 1).
In 2011, the correlation of lint yield with NDVI was significant at early square and early, mid, and late bloom stages (Table 1). The correlation of lint yield with leaf N was significant at mid and late bloom stages (Table 1). There was a significant correlation of leaf N with NDVI at mid and late bloom stages (Table 1). Overall, the R2 values for the above correlations in 2011 were similar to those in 2010, but lower than those in 2009; which suggests that the correlations of lint yield and leaf N with NDVI varied with year.
The University of Tennessee recommends 67 to 90 kg N ha−1 per year for cotton grown in Tennessee [10] . Application of chicken litter before cotton planting is a common practice to some producers in Tennessee and some other states in the Mid- South region. Application of 45 kg ha−1 of N to all the plots in the form of chicken litter as pre-plant N before cotton planting in this experiment accounted for a relatively high percentage of the annual N fertilizer recommendation for cotton, which probably weakened the correlations of lint yield and leaf N with NDVI because the application of pre-plant N increased N supply to the plant, and thus decreased the responses of lint yield, leaf N, and NDVI to the side dress N treatments. Our results suggest that the relationships of lint yield and leaf N with canopy NDVI may not be strong enough to be used to assess plant N nutrition status and predict cotton yield during the growing season for guiding in-season side dress N application. If we want to use the relationships of lint yield and leaf N with NDVI to guide in-season side dress N application, we may need to change the N management systems for cotton by reducing or even avoiding pre-plant N application but increasing side dress N application during the growing season.
3.2. Geostatistical Analyses of Field Spatial Variability of Lint Yield
The Arc View maps of lint yield are presented in Figures 1-3 for 2009, 2010, and 2011, respectively. The lint yield maps showed that although N application rate had impacts on lint yield, spatial variation in lint yield did exist within most strip plots each year, and it varied with year.
We conducted a quadratic regression of lint yield with side dress N rate treatments using the classic model in the GeoDa software in order to examine the spatial dependence of lint yield within the test field. Significant spatial dependence of lint yield
![]()
Figure 1. ArcView map of lint yields at harvest in the N experiment at Gibson in 2009.
![]()
Figure 2. ArcView map of lint yields at harvest in the N experiment at Gibson in 2010.
within the field was observed each year (data not presented). Then the spatial error model in GeoDa was used to perform the quadratic regression of lint yield with side dress N treatments; the output is presented in Table 2. It showed that the quadratic relationship of lint yield with side dress N rate treatments was significant on a sub-plot
![]()
Figure 3. ArcView map of lint yields at harvest in the N experiment at Gibson in 2011.
![]()
Table 2. Regression summary of output using spatial error model in the N experiment at Gibson in 2009-2011.
basis in 2009 and 2010, but not in 2011. Our results suggest that use of the spatial error model in the statistical analysis of N treatment effect on cotton yield is more appropriate in two out of the three years.
In order to visualize the spatial dependence of lint yield relating to the characteristics of the test field (not to the N treatments), we used the residual lint yield data (which were obtained in the spatial error model in GoeDa and in which N treatment effects on lint yield had been excluded) to make the Moran’s I statistic and scatter plot and LISA cluster map each year.
The Moran’s I and scatter plot is used to evaluate global spatial autocorrelation. Moran scatter plot provides a visual exploration of global spatial autocorrelation. The four quadrants in the Moran scatter plot provide a classification of four types of spatial autocorrelation: high-high and low-low for positive autocorrelation; low-high and high- low for negative autocorrelation. The value listed at the top of the graph is the Moran’s I index. In 2009, Figure 4 showed that the spatial autocorrelation of residual lint yield (N treatment effects on yield excluded) was not significant (p = 0.623) within the test field. However, in both 2010 and 2011, there was significant (p = 0.003 in 2010 & 0.001 in 2011) spatial autocorrelation of residual lint yield within the field (Figure 5 & Figure 6).
The LISA cluster map is a tool for estimating local spatial autocorrelation. It contains information on only those locations that have significant spatial autocorrelation. Four types of spatial autocorrelations are colored in four different colors: dark red for high-high, dark blue for low-low, pink for high-low, and light blue for low-high. In 2009, the LISA cluster map showed that there were some significant local clusters of residual lint yield (N treatment effects on yield excluded) within the test field (Figure 7). Specifically, there were six sub-plots with high residual yield surrounded by high residual yield neighbors, four high residual yield sub-plots were surrounded by low residual yield neighbors, and two sub-plots with low residual yields were surrounded by
![]()
Figure 4. Moran’s I and scatter plot of residual lint yield (N treatment effects on yields excluded) in the N experiment at Gibson in 2009.
![]()
Figure 5. Moran’s I and scatter plot of residual lint yield (N treatment effects on yields excluded) in the N experiment at Gibson in 2010.
![]()
Figure 6. Moran’s I and scatter plot of residual lint yield (N treatment effects on yields excluded) in the N experiment at Gibson in 2011.
high residual yield neighbors. In 2010, some significant local clusters of residual lint yield were also observed within the field in the LISA cluster map (Figure 8). There were six sub-plots with high residual yield surrounded by high residual yield neighbors, two low residual yield sub-plots were surrounded by low residual yield neighbors, seven sub-plots with low residual yield were surrounded by low residual yield neighbors, and two high residual yield sub-plots were surrounded by low residual yield neighbors. In 2011, the LISA cluster map showed that there were some significant local clusters of residual lint yield within the field (Figure 9). Specifically, there were eighteen sub-plots with high residual yield surrounded by high residual yield neighbors, sixteen low residual yield sub-plots were surrounded by low residual yield neighbors, four sub-plots with low residual yield were surrounded by high residual yield neighbors, and two high
![]()
Figure 7. LISA cluster map of lint yield (N treatment effects on yields excluded) in the N experiment at Gibson in 2009.
![]()
Figure 8. LISA cluster map of lint yield (N treatment effects on yields excluded) in the N experiment at Gi- bson in 2010.
![]()
Figure 9. LISA cluster map of lint yield (N treatment effects on yields excluded) in the N experiment at Gibson in 2011.
residual yield sub-plots were surrounded by low residual yield neighbors. Overall, our results suggest that both global and local spatial autocorrelation varied with year in this study.
Traditionally, a classical randomized complete block and similar designs were often used in agricultural field experiments with small plots. Replications of the treatment plots were used at a location as homogeneous as possible so as to avoid the influences of spatial variability. A widely accepted assumption was that existing spatial variability in the experiment could be compensated for by more replications of the treatment plots [11] . Because a strip plot covers a much larger area then a small plot, the spatial variation within a strip plot is likely greater than that within a small plot. Therefore, more replications may be required for experiments with strip plots than those with small plots to compensate for the same degree of spatial variability in the field. Another alternate way to cope with the possible larger spatial variation with the strip plots is to measure the spatial variation within each strip plot and then included it appropriately in the statistical analyses of the treatment effects for the strip plot study. Since this experiment was carried out on the identical plots during the three-year experimentation, the N treatment effect on cotton yield of the second and third years (2010 & 2011) was cumulative if there was any residual N treatment effect left from the previous year(s).
4. Conclusions
The relationship of cotton yield and leaf N with canopy NDVI was mostly statistically significant but not strong during early square to late bloom irrespective of year. The uniform application of 45 kg ha−1 of N to all the plots in the form of chicken litter as pre-plant N before cotton planting accounted for a high percentage of the annual N fertilizer recommendation for cotton, and thus probably weakened the correlations of lint yield and leaf N with NDVI. Our results suggest that the relationships of cotton yield and leaf N with canopy NDVI may not be strong enough to be used to assess plant N nutrition status and predict cotton yield during the growing season for guiding in-season side dress N applications after a common rate of N was applied as pre-plant N in the form of chicken litter before cotton planting. In order to use the relationships of cotton yield and leaf N with NDVI to guide in-season side dress N application, the N management systems for cotton need to be changed by reducing or even avoiding pre-plant N application but increasing side dress N application rate during the growing season.
Spatial variability was present in cotton yield at the strip plot experiment scale on a farmers’ field. However, the pattern and degree of the spatial variability varied with year. Our results suggest that the spatial variability of cotton yield is high, and needs to be taken into account appropriately in the statistical analyses of N treatment effects on cotton yield in strip plot field studies. On the other hand, our results indicate that cotton yield from farmers’ fields can be expected to have significant annual and within field spatial variation in the region, and this variation will significantly influence cotton yields. If the spatial variability is used appropriately in precision N management, producer profitability can be increased.
Acknowledgements
This project was partially supported by Cotton Incorporated Cooperative Agreement No. 09-497TN, managed by Dr. Bob Nichols. We appreciate the cooperative farmer: Jeff Dodd in Gibson County for allowing us to conduct the test on his farm. We also thank the technical cooperation of the Textile Service Laboratory of Cotton Inc. Technical assistance was provided by Robert Sharp, James Warren, Tracy Bush, Matt Ross, Dereck Eison, Ngowari Jaja, Will Goforth, and others.