Published in Agron. J. 96:77-90 (2004).
© American Society of Agronomy
677 S. Segoe Rd., Madison, WI 53711 USA
RICE
Spatiotemporal Analysis of Rice Yield Variability in Two California Fields
Alvaro Roela,c and
Richard E. Plant*,b
a Graduate Group in Ecol., Univ. of California, Davis, CA 95616, USA
b Dep. of Agron. and Range Sci. and Dep. of Biol. and Agric. Eng., Univ. of California, Davis, CA 95616, USA
c Instituto Nacional de Investigaciones Agropecuarias, Treinta y Tres, Uruguay
* Corresponding author (replant{at}ucdavis.edu).
Received for publication November 19, 2002.
 |
ABSTRACT
|
|---|
Currently, little is known about the spatial and temporal variability of rice (Oryza sativa L.) yield patterns. This information is needed before implementing any site-specific management strategy. The objective of this research was to characterize the spatial and temporal yield variability of rice grown in commercial fields in California. Rice cultivars M-202 and Koshihikari were grown and managed by a cooperating farmer, who collected yield monitor data over a 4-yr period. Alternative methods of data quality analysis were applied to the data. To evaluate temporal variability, yields from different years must be placed on a common grid. The appropriate size for these grids was tested. Large-scale spatial structure was determined using median polish while small-scale spatial structure was evaluated by computing variograms of the yield residuals after subtracting the trends. Temporal variability was determined using two approaches: (i) computation of the variance among years of the original, trend, and residual yield values at fixed points in the field and (ii) cluster analysis of the standardized trend yield values. Results from the study showed that the grid density necessary to capture the spatial variability depended on site and year. Trend surface spatial behaviors depended on year, indicating a lack of temporal stability. Variograms showed strong spatial structure of yield residuals. Cluster analysis reduced the considerable complexity in a sequence of yield maps of these fields to a few general patterns of variations among years.
Abbreviations: GR, grid resolution MC, Moran coefficient RIC, relative information criterion
 |
INTRODUCTION
|
|---|
YIELD MONITORING and mapping technology that can measure, georeference, and record grain yields makes it possible to document the location and magnitude of yield variability with a spatial precision of meters. If the causes of this variability can be identified, then corrective action may be implemented to reduce costs, increase yield, and/or reduce environmental impacts by adopting site-specific management practices (Lowenberg-DeBoer and Erickson, 2000). Moreover, the availability of high-precision measurements may permit researchers to more efficiently test hypotheses by precisely measuring crop response to environmental conditions as these conditions vary in the field (Bhatti et al., 1991; Grondona and Cressie, 1991; Long, 1998). Both of these uses of precision-measurement technology require statistical methods that until now have more commonly been employed in ecological, epidemiological, and econometric research (Long, 1998; Griffith and Layne, 1999; Bongiovani and Lowenberg-DeBoer, 2001).
Yield map data quality is an important initial issue for farmers and even more so for scientists and engineers who wish to use these data in the course of a scientific study. Blackmore and Marshall (1996) identify six primary sources of error in yield map data: unknown crop width at the header, time lag in the harvester, GPS error, grain surge, grain losses, and sensor accuracy and calibration. Birrell et al. (1996), Blackmore and Marshall (1996), Doerge (1999), Colvin and Arslan (2001), and Haneklaus et al. (2001) provide discussions of the errors associated with yield maps.
It is a well-known property of spatial data that their statistical properties depend on the scale at which they are represented. This is often called the modifiable areal unit problem (Openshaw and Taylor, 1981; Wong, 1995).Yield monitor data are generally recorded periodically (e.g., once per second). To compare data from different sources (e.g., yield monitor data from different years, bulk soil electrical conductivity data, and aerial image data), these data must be converted to a common grid. In choosing the grid size, there is a trade-off between maintaining spatial precision by selecting a fine grid and reducing noise and making the data more manageable by selecting a coarser grid (Wong, 1995; Long, 1998). Since variability may be studied at any spatial scale, the choice of grid size depends on the aims of the investigation. In making this choice, the investigator is aided by knowledge of how much information is lost in moving from one scale to a larger one. One way of determining the effect of increasing grid size is to examine the experimental variograms of the measured quantities (Isaaks and Srivastava, 1989). Long (1998) examined the change in correlation coefficients between yield and a second variable as the grid size increased. Chen et al. (1995) developed the relative information criterion (RIC) to quantify efficiency of data representation. They defined the RIC as the correlation coefficient between block-kriged estimates based on the highest grid density and estimates obtained at the same points based on reduced grid densities.
One important application of yield map data is the study of the spatiotemporal properties of yield distribution. These properties should be understood before implementing any site-specific management strategy (Bakhsh et al., 2000). From a scientific perspective, understanding spatiotemporal variation in yield may help to determine the factors underlying yield variability. Several studies of spatiotemporal patterns in terrestrial crops have been performed. In corn (Zea mays L.) and soybean [Glycine max (L.) Merr.], Jaynes and Colvin (1997) reported that yields displayed substantial spatial and temporal variability. This variability may be due to interactions among climatic growing conditions, soil properties, and the crop (Porter et al., 1998).
There have been few studies on the spatiotemporal variability of rice yields (Dobermann et al., 1995). Rice is one of the world's most important staple crops. The development of effective site-specific management techniques for rice production could, if it increases production efficiency, contribute to increased yields (Dobermann et al., 2002). The study of spatiotemporal variability in rice production also provides scientifically useful information as a comparison with terrestrial crops. Eghball and Power (1995) found that rice yields displayed less year-to-year variability than terrestrial crops. Since rice is grown as a monoculture in California, this production system provides a particularly good environment to analyze the evolution of yields in space and time.
Yield variability is often defined in terms of summary statistics such as temporal variance and spatial variance (Whelan and McBratney, 2000). To fully characterize the spatiotemporal behavior of a field, however, one must also understand the tendency of yield patterns to persist season after season, that is, the temporal stability of the spatial pattern. Correlation coefficients between years are often relatively small (Jaynes and Colvin, 1997). Moreover, determining the significance of a correlation coefficient is complicated due to the effect of spatial autocorrelation of the data (Cliff and Ord, 1981). One way to evaluate temporal variability is to compute the temporal variances of yield values at fixed points in the field (Whelan and McBratney, 2000). By comparing the variances among seasons (temporal) with those within a season (spatial), the importance of both components can be estimated. Cluster analysis of annual yields can also be used to assess temporal patterns (Lark and Stafford, 1997; Lark et al., 1999; Perez-Quezada et al., 2003). Cluster analysis may provide an objective quantification of the spatial structure of yield patterns as well as an indication of the consistency of these patterns from year to year (Perez-Quezada et al., 2003).
Our study was performed to perform a spatial and temporal analysis of 4 yr of rice yield monitor data collected in two California fields. The first objective of this study was to examine different methods for evaluating the quality of the data set. The second objective was to determine the grid resolution (GR) that captures enough information to represent yield spatial variability at a scale appropriate to management. The third objective was to assess the usefulness of summary statistics in quantification of the spatiotemporal variability. Finally, the fourth objective was to examine the use of clustering to delineate areas in the field with different spatiotemporal yield behaviors.
 |
MATERIALS AND METHODS
|
|---|
The study was performed from 1998 through 2001 in two rice fields approximately 2 km apart, one of 38 ha (denoted Field 1) and one of 52 ha (denoted Field 2), located near Marysville, CA (UTM Zone 10 coordinates: E: 627102, N: 4340769; and E: 624970, N: 4341076 for Field 1 and 2, respectively). The soils of the study fields consist of approximately 45% Kimball loam (fine, mixed, active, thermic Mollic Palexerafls), 30% San Joaquin loam (fine, mixed, active, thermic, Abruptic Durixeralfs), and 25% Bruella loam (fine-loamy, mixed, Ultic Palexeralfs). Medium-grain rice cultivar M-202 and short-grain cultivar Koshihikari were grown and managed by the cooperator in Fields 1 and 2, respectively, using standard practices for the area (Hill et al., 1992). The fields were selected based on conversations with the cooperating grower, with the object that one field would tend to have spatially uniform properties and the second would be very heterogeneous. Figure 1 shows gray-scale renditions of false color infrared aerial images taken of the two fields during the first year of the experiment. These images were taken using Kodak 2443 infrared film on 18 Aug. 1998 and are shown before georeferencing. The dark areas in Field 2 correspond to areas of very sparse vegetation. The field had recently been laser-leveled and brought into production, and these dark areas in the image were areas where considerable topsoil had been cut off in the leveling process. The fields were managed without regard to spatial variability with one exception: The consistently poor-yielding portion of Field 2 was observed by the grower to mature earlier than the rest of the field, and so the grower harvested this area first. In 2001, the grower harvested this portion of the field 5 d earlier than the rest of the field.

View larger version (104K):
[in this window]
[in a new window]
|
Fig. 1. Gray-scale renditions of false-color infrared aerial images taken on 8 Aug. 1998 of the two study fields. (a) Field 1. Note that in the image the field becomes darker as one moves from west to east. This corresponds to higher near-infrared values in the original image. (b) Field 2. The dark areas in this image correspond to areas of very sparse vegetation in the field. The field had recently been laser-leveled and brought into production, and these dark areas in the image were areas where considerable topsoil had been cut off in the leveling process.
|
|
Rice was harvested using a stripper harvester combine equipped with a John Deere Green Star yield-mapping system with real-time differential global positioning system (DGPS) receiver in all years. The harvester followed a back-and-forth, north-to-south harvest pattern in Field 1 and a series of concentric patterns in Field 2. Combine speed ranged from 1.1 to 1.25 m s1, header width ranged from 5.5 to 6.7 m, grain flows were recorded once per second, and moisture content was recorded once per 15 s. Yield map data files (yield, grain, moisture, longitude, and latitude) were collected and imported into the ArcView (ESRI, Redlands, CA) geographic information system for analysis. The yield monitor was calibrated at the start of each harvest season by comparing with a sample of known weight, and yield measurements were adjusted based on this calibration. Yield monitor data points for a distance of approximately 50 m from field edges were deleted from the data set to remove border effects and end-of-field yield monitor errors. Aerial images were taken approximately mid-July and mid-August of each year, using Kodak 2443 infrared film in the first year and a four-band multispectral digital camera in each other year, and were georeferenced to boundary files made using a Trimble Ag 132 GPS.
Data Quality Analysis
Yield monitor data were imported into ArcView GIS as point shapefiles (each point representing the grain flow and moisture measurements in an area of size approximately 1 by 6 m). The GPS accuracy and consistent header width were verified by visual inspection of the records. The John Deere Green Star yield monitor contains built-in proprietary software to correct for time lag in the harvester (D. Goebel, John Deere Co., personal communication, 2003), so we did not attempt to further correct for this. We could not measure errors due to grain surge or losses, so this left sensor calibration errors as the remaining source of inaccuracy we could detect.
Grain flow rates and moisture content were converted to yield at constant 14% moisture. Considering the sensor record as a time series, we tested for two ways in which sensor calibration could change during the course of the harvest of a single field: a gradual trend and a sudden change. A convenient method of trend analysis in time series data is by studying the differences between successive records (Kendall and Ord, 1990). The data sets consisted of between approximately 50000 and 90000 values. We first removed outliers by visual inspection of the data record. To make the data more manageable and remove the statistical problems associated with very large data sets (Matloff, 1991), we next selected every 10th point from the data records for time series analysis.
The resulting data sets each consisted of a sequence of time series in which discontinuities in the GPS clock time indicated points in which the harvest had been stopped and then restarted. Yield data were plotted against GPS clock time. The resulting plot was then inspected first for sudden changes in calibration and second for evidence of a trend (Haneklaus et al., 2001). We assumed that sudden changes in calibration would occur at discontinuity in the GPS clock time (indicating that the harvest had been stopped and restarted). Because in Field 2 the cooperating grower harvested the low-yielding areas first and then the better yielding areas, an abrupt change in yield at a gap in GPS time in Fig. 2 does not necessarily represent a calibration artifact. To determine which if any of the records showed an evident change in yield monitor calibration the yield records were visually compared with false color infrared aerial images of the field, with a change in calibration being indicated if the change in yield did not correspond to a change in vegetation intensity in the aerial image.

View larger version (42K):
[in this window]
[in a new window]
|
Fig. 2. Yield monitor data (every 10th data point) of Field 1 plotted in the order in which records were entered into the database. Abscissa is the difference between the GPS clock time for that record and the lowest GPS time of the data set. Ordinate is yield (kg ha1) converted to common moisture content. Arrows placed at gaps in the GPS time record indicate that the harvest had been stopped and restarted.
|
|
The comparison was performed as follows. In each year and for each field, the raw data values were displayed as points in the GIS. The display was examined for evidence of abrupt changes in yield trend that occurred at discontinuities in the GPS time. A false color aerial image of the field was then examined to determine if there was a corresponding change in infrared reflectance at this location. Because of the difficulties in estimating yield from aerial images (Plant et al., 2001), this process was done visually rather than statistically. If there was no change in reflectance properties corresponding to a change in yield, the change was assumed to be caused by a sudden change in calibration of the yield monitor. In this case, the yield data after the jump in GPS time was adjusted by multiplying by a constant value to bring the yield trends before and after the change into visual alignment.
To test for gradual drift or trend in the data, the yield records were differenced, and the differences were tested against the null hypothesis of zero mean using the sign test. The third data quality test was performed in Field 1 to take advantage of the back-and-forth harvest pattern. The test was performed at 10 randomly selected locations in the field. It consisted of comparing the mean of a sequence of 10 yield values with the mean of the sequence of 10 values at points immediately to the left of the first 10. Confidence intervals were computed for the means based on the standard deviation of the 10 sequential yield values. These confidence intervals are not exact due to the autocorrelation of the data.
Spatial Resolution Analysis
Yield point shapefile data were interpolated to a fixed 5- by 5-m grid using inverse-distanceweighted interpolation with power 2 and number of neighbors 12. This grid provided a set of locations of the yield data that was consistent over the 4 yr and approximated the spacing of the original yield data. These interpolated grids were used as the starting point for the analyses. Our primary interest was in studying variability on a scale of tens of meters to eliminate very short-range effects while still maintaining the ability to distinguish important patterns of spatial variability. Therefore, three different regularly spaced square point grids of size 90, 60, and 30 m were used to extract the values from the interpolated yield maps in both fields. These three grids form resolutions with numbers of grid points: n = 36, n = 76, and n = 292 and n = 42, n = 93, and n = 402 for the 90-, 60-, and 30-m grid in Field 1 and 2, respectively. Each larger grid was made up of a subset of the smaller grids.
Large-scale variability (trend) and small-scale variability were separated by detrending the data, using the median polish technique (Cressie, 1991; Jaynes and Colvin 1997; Bakhsh et al., 2000). Median polish by rows and columns may not capture the entire large-scale trend as the trend orientation is not known a priori (Cressie, 1991). Therefore, an additional term was included in the median polish equation to detect any further trend in the polished data as described by Cressie (1991) and Jaynes and Colvin (1997). No significant further trends were detected in any case in which this procedure was used in our data sets.
For the characterization of the small-scale variation, experimental variograms were calculated using the residual yields from the median polish. Variography analyses assume the data have a Gaussian distribution. Therefore, the distributions of the yield residuals for the 4 yr were checked using stem and leaf plots and normal probability plots. Outliers, indicated by stem and leaf plots, were removed and replaced by kriged estimates following the analysis of Bakhsh et al. (2000). In Field 1, 15, 22, 11, and 12 values for the 30-m grid and 9, 4, 0, and 0 values in the 60-m grid were removed in the 1998, 1999, 2000, and 2001 yield data sets, respectively. In Field 2, 11, 13, 9, and 22 values for the 30-m grid were removed in the 1998, 1999, 2000, and 2001 yield data sets, respectively. In the 60-m grid, 4 values in 1998 and 13 in 2001 yield data sets were removed. Experimental variograms were calculated from the yield values. All sample variogram computations and model fitting were performed using GS+ v. 5.0 software (Gamma Design Software, Plainwell, MI) assuming isotropic conditions. The isotropic assumptions were verified by variography analysis in different orientations. Only lags less than half the total grid length were computed. A theoretical variogram model was fit to each experimental variogram. Different models were tested for fitting the data (Isaaks and Srivastava, 1989), and the isotropic spherical model provided a good fit in each case. One method for estimating the appropriate resolution for spatial analysis was to base it on the ranges of the fitted variograms.
A second method for determining the appropriate cell size for analysis is the RIC (Chen et al., 1995). In this paper, we use the square of the RIC as defined by Chen et al. (1995). The RIC2 is the square of the product moment correlation coefficient of kriged residuals from the same locations in the field, namely:
In this equation, r2(kj,ki) denotes the square of the correlation coefficient of the block-kriged yield residual estimates using GR i and j. The RIC2 as we use it gives the fraction of the sample variation using GR j that is explained if it is related to GR i using a simple linear regression.
Block kriging was used to interpolate the residuals from the different grids. Following Chen et al. (1995), block kriging with a block size equal to a 4- by 4-m plot was used to produce block means of yield residuals, from which distribution maps of yield residuals were developed. The block size was chosen to approximate as closely as possible the original cell size of the interpolated grid, which could not be fit exactly due to limitations of the software. The search radius in the kriging procedure was 350 m in both fields. In this procedure, we first estimated the variogram as described above and block-kriged using all locations in both fields, corresponding to a GR of 30 m (i = 3). This procedure was repeated at GR of 60 and 90 m (i = 2 and 1).
Computation of Summary Statistics
All temporal analysis was performed on the 30-m grid. The four seasons worth of data do not provide the opportunity to study large- vs. small-scale temporal variability; any long-term variability would only be evident with a longer temporal record (Eghball and Power, 1995). We therefore focused on the stability of the spatial structure over the 4 yr. We computed the mean temporal variance and the mean spatial variance of each field. These statistics are defined as follows. Let N be the total number of grid cells and T the total number of years. Let Y(n,t) be the yield in grid cell n in year t, and let
T
be the mean yield of cell n over all years. The temporal yield variance
2T
of cell n is the variance in yield over the T years,
 | [1] |
and the mean temporal variance is the average of these grid cell variances over all cells,
 | [2] |
The spatial yield variance is the yield variance over all grid cells,
 | [3] |
where
N
is the mean yield in year t over the N cells, and the mean spatial variance is the average of these spatial variances over all years, given by
 | [4] |
Finally, it is of interest to determine the extent to which temporal variability is distributed over a field. The summary statistic for this is the standard deviation of the temporal variance, given by
 | [5] |
Cluster Analysis
We used the k-means clustering algorithm (Jain and Dubes, 1984). This algorithm forms a fixed number k of cluster groups by selecting those clusters that minimize within-group variances and maximize the between-group variance. Cluster analysis was performed on the trend data obtained by median polish. Yields were standardized by subtracting the mean and dividing by the standard deviation to avoid emphasis on data from higher-yielding years. Cluster analysis was performed using Statistica (StatSoft, Tulsa, OK).
 |
RESULTS AND DISCUSSION
|
|---|
Data Quality Analysis
Figures 2 and 3 show the yield monitor records presented as time series, or more precisely, as connected sequences of time series, plotted in the order in which the data were entered into the record. The ordinate of each axis is the difference between the GPS time in seconds of the corresponding record and the lowest GPS time of all the records in the data set. Arrows mark the points at which there is an abrupt change in the GPS time, indicating that the harvest had been stopped and resumed. Visual inspection of the figures indicates no obvious anomalies in any of the data sets except that of Field 2 in 2001. This data set shows evidence of poor quality, including apparently truncated values in the early part of the data set. Based on this evidence, we elected to remove the Field 2 2001 data set from most of the analysis.

View larger version (47K):
[in this window]
[in a new window]
|
Fig. 3. Yield monitor data (every 10th data point) of Field 2 plotted in the order in which records were entered into the database. Abscissa is the difference between the GPS clock time for that record and the lowest GPS time of the data set. Ordinate is yield (kg ha1) converted to common moisture content. Arrows placed at gaps in the GPS time record indicate that the harvest had been stopped and restarted.
|
|
Based on the analysis for abrupt changes in calibration, we concluded that the abrupt changes in Field 2 all corresponded to actual abrupt changes in yield but that the abrupt change in yield in Field 1 in 1998 at the second GPS time gap did represent a calibration change. This record was adjusted accordingly. Specifically, all data values occurring after the second GPS time gap were multiplied by 1.37, which brought the data set into visual alignment.
The second test for data quality was the check for gradual drift of the calibration. Tests of the null hypothesis m = 0, where m is the median of the yield differences yt+1 yt, were not significant (p > 0.05) for Field 1 in any year. They were significant (p < 0.05) for Field 2 in 1999 in 2001. The 2001 data set had been removed from analysis already. Comparison of 1999 data with the corresponding aerial image led us to conclude that the trend was due to actual conditions in the field rather than sensor drift.
Figure 4 shows the means and confidence intervals for the years with the most (five in 2000) and fewest (two in 1998) overlapping confidence intervals in 10 randomly selected sequences of 10 points. In general, the agreement between subsequent passes of the yield monitor at the same location was poor. This same high level of local spatial variability in the data may be observed in the data presented by Searcy et al. (1989) for grain yield and was reported by Searcy et al. (2001) for cotton (Gossypium hirsutum L.) yield and mentioned by Perez-Munoz and Colvin (1996) for grain yield. It was also observed in data recorded for a variety of crops by Perez-Quezada et al. (2003) and indicates that yield monitor data may be an imprecise predictor of yield at specific locations in the field.

View larger version (14K):
[in this window]
[in a new window]
|
Fig. 4. Relationship between the means of sequences of 10 yield monitor data points recorded next to each other in Field 1. Error bars show 95% confidence intervals based on standard deviation of the 10 data values.
|
|
Grid Resolution
Figures 5 and 6 show the sequences of yield maps of the 4 yr for Field 1 and 2, respectively, with data corrected as described in the Materials and Methods section and interpolated to a 30-m grid. Yield in Field 1 has no consistent visually apparent pattern over the years while in Field 2, there are persistent low- and high-yielding areas. The low-yielding areas correspond to cut areas of the laser-leveling process, which are the dark areas in Fig. 1b. In both fields, the 30-m grid appeared to visually correspond to patterns of spatial variability observed in infrared aerial images of the field.

View larger version (59K):
[in this window]
[in a new window]
|
Fig. 5. Sequence of maps of Field 1 data as analyzed for each of the years 1998 through 2001. Yield data are in kg ha1. Outliers have been deleted, and data have been resampled to a common grid representing 30-m square cells (Grid 1).
|
|

View larger version (61K):
[in this window]
[in a new window]
|
Fig. 6. Sequence of maps of Field 2 data as analyzed for each of the years 1998 through 2001. Yield data are in kg ha1. Outliers have been deleted, and data have been resampled to a common grid representing 30-m square cells (Grid 1).
|
|
Figure 7 shows the variograms of the median polish residuals for the 30- and 60-m GR of the 1998 data sets in Field 1. The variograms for Field 2 are almost identical. These are consistent with the variograms for the other years, which are not shown. To facilitate comparison, the sill and nugget are scaled to the sample variance for each variogram. The variogram for the 90-m grid, which consisted of only three points, displayed no structure (pure nugget). The 90-m grid was therefore not considered adequate for representing the spatial structure of yield variability.

View larger version (11K):
[in this window]
[in a new window]
|
Fig. 7. Experimental and fitted variograms of yield data in Field 1. Values of the fitted variogram sill, co + c1; range, Ao; and model fitness, r2; are given in the figures. (a) 30- and (b) 60-m sampling grid. Variograms for Field 2 are very similar.
|
|
Table 1 gives the parameters of the variogram fit to the median polish residuals using the spherical model for each of the years for the 30- and 60-m GRs for Fields 1 and 2, respectively. Table 2 gives the RIC2 values for each field. In both fields, RIC2 values vary from year to year in a manner that corresponds to the variogram ranges displayed in Table 1. The values in Table 2 may be interpreted as the portion of the variance of the data in the 30-m grid explained by a simple linear regression of 60 m on the 30-m grid. Comparison of the data in Table 1 with the data of Table 2 indicates that efficiency of representing the grid is directly related to the spatial continuity of the variable, i.e., that the 60-m grid captures more of the information of the finer grid in cases where the range of the variogram is smaller. Our results agree with those of Chen et al. (1995). Since RIC2 analysis indicated that a substantial amount of information may be lost in some years by reducing the grid density, all further analysis were done only with the 30-m grid.
View this table:
[in this window]
[in a new window]
|
Table 1. Range, sill, and nugget of variograms fit using the spherical model to 30- and 60-m grid resolutions of the yield median polish residuals for both fields.
|
|
View this table:
[in this window]
[in a new window]
|
Table 2. Squared relative information criteria (RIC2) between the 30- and 60-m grid resolutions for Fields 1 and 2 for each of the 4 yr.
|
|
Summary Statistics
Table 3 contains the spatial variances
2s
of the total yield, the median polish, residuals, and the covariance between median polish and residual values. The variances of the residuals, which are a combination of sampling error and short-range effects, are of the same magnitude as those of the median polish. There is a negative relationship between large- and small-scale values in Field 2, reflecting a tendency of the lower-yielding areas to have smaller short-range variability. Table 4 gives the mean temporal variance (Eq. [2]), standard deviation of the temporal variance (Eq. [5]), and the ratio between the mean temporal and the mean spatial variance (Eq. [4]) for the original, trend, and residual yield values. Computations for Field 2 were performed using the years 19982000 only. Mean temporal variance values over the seasons are larger than the mean spatial variances for the original and trend (large-scale) yield values in Field 1. Field 2 presented both higher spatial and temporal variance values.
View this table:
[in this window]
[in a new window]
|
Table 3. Spatial variance in total yield, yield trend obtained through median polish, and small-scale yield variation obtained as median polish residuals, and covariance between median polish and residuals. Units are (Mg ha1)2.
|
|
View this table:
[in this window]
[in a new window]
|
Table 4. Mean temporal variance (Eq. [2]), standard deviation (std. dev.) (Eq. [5]), and the ratio between the mean temporal and the mean spatial variance (Eq. [4]) for the years 19982001 for Field 1 and 19982000 for Field 2.
|
|
Figures 8a and 8b show the spatial distribution of the temporal variance (Eq. [1]) computed with the trend yield values in Field 1 and 2, respectively. To visualize different variance behaviors, the magnitude of the variances were classified in three different categories: values less than 1, between 1 and 2, and greater than 2 (Mg ha1)2. These may be considered as representing low, intermediate, and high variability behaviors, respectively. To aid in visualization, the figures are drawn using Thiessen polygons surrounding each grid point. Visual inspection of the figures shows that the number of locations with high temporal variability was very small in Field 1 and that most of the points with high temporal variability were concentrated in the southwest corner of the field. Most of the rest of the field (184 of the 292 polygons, or 63%) had a small temporal variance (<1 Mg ha1). Figure 8b indicates that Field 2 displays very different behavior from Field 1. In this field, most of the polygons (296 of the 402, or 73%) had a large temporal variance (>2 Mg ha1).

View larger version (40K):
[in this window]
[in a new window]
|
Fig. 8. Spatial distribution of the temporal variance (Mg ha1)2 for (a) Field 1 and (b) Field 2. Each figure consists of Thiessen polygons surrounding each sample point.
|
|
Cluster Analysis
Figure 9 shows the results of the k-means cluster analysis, with k = 3. The figure gives the values of the standardized yields in each of the clusters. In Field 1, Cluster 1 is composed of locations of the field with yield performance above field yield average in 1998, around field yield average in 1999, and below field yield average in 2000 and 2001; Cluster 2 is composed of locations in the field that presented low yield performance in 1998 and 1999, average performance in 2000 and above-average performance in 2001; and Cluster 3 is composed of locations in the field with yield performances similar to or above the annual field average in all years. The three clusters in Field 2 consist of a consistently high-yielding area, a consistently moderate-yielding area, and a consistently low-yielding area (Fig. 9b).

View larger version (15K):
[in this window]
[in a new window]
|
Fig. 9. (a) Cluster behaviors defined by standardized yields in 4 yr in Field 1. (b) Cluster behaviors defined by standardized yields over the three years 19982000 in Field 2.
|
|
Figures 10a and 10b show the spatial location of the three different cluster groups in Fig. 9a and 9b, respectively. These figures allow the visualization of locations in the field classified in each cluster and give an estimation of the importance of the spatial configuration of these clusters. The clustering procedure does not use any information regarding the spatial location of the data sets. Nevertheless, Figure 10 shows that the three different cluster groups tend to be spatially grouped. The Moran coefficient (MC) can be used to test the null hypothesis of zero spatial autocorrelation (Long, 1998). The distribution of cluster values was significantly different from random in both fields (MC = 0.673, Z value = 72.0, and p
0.0029 in Field 1 and MC = 0.620, Z value = 69.0, and p
0.0031 in Field 2).

View larger version (48K):
[in this window]
[in a new window]
|
Fig. 10. (a) Thiessen polygons interpolated clusters in Field 1 (19982001). (b) Thiessen polygons interpolated clusters in Field 2 (19982000).
|
|
Comparing Fig. 5, 9a, and 10a shows that although Field 1 appears to have highly disorganized yield behavior, there are yield patterns that emerge from the cluster analysis. Comparing Fig. 6, 9b, and 10b indicates that the organization of clusters in this field closely matches the consistent pattern observed in the original data. This indicated that yield patterns in Field 1 were associated with transient phenomena while yield patterns in Field 2 were associated with more permanent factors.
 |
SUMMARY AND CONCLUSIONS
|
|---|
The study of spatial yield variation of field crops in commercial fields may be addressed at several scales. Yield data collected with commercial yield monitor clearly cannot be used to study variation at a scale smaller than the resolution of the harvesting equipment. The results of our study, in agreement with those reported by others, indicate that commercial yield monitors may not provide useful information at the scale of one or two multiples of the width of the harvester. Nevertheless, our results provide an indication that at a somewhat larger scale (30 m in our case), the commercial yield monitor can provide useful information about patterns of relative yield. In Field 2, the pattern of consistently low yield in certain areas matched observations of aerial infrared photographs and was explainable in terms of the laser-leveling history of the field. The behavior in Field 1 was more complex and dynamic, but cluster analyses indicated a strong spatial pattern.
Our analysis was based on interpolating yield monitor data into a regular grid and then examining the properties of this grid. This is a commonly used procedure (Birrell et al., 1996; Swindell, 1997; Long 1998; Perez-Quezada et al., 2003). Although this is technically a geostatistical procedure, and has been analyzed by Long (1998) as such, the dense pattern of the data also justify considering it from a lattice perspective. In this context, each data point may be considered as an estimate of the mean of the yield distribution within a lattice cell whose width is the width of the harvester and whose length is the distance traveled. The interpolation process then becomes one of changing the scale and extent of the lattice cells to a regular grid. Openshaw and Taylor (1981) and Long (1998) have investigated the effects of such changes and have shown that they profoundly affect the autocovariance structure. Thus, conclusions about spatial autocorrelation of such resampled data apply only to aggregation at the scale in which the study was performed.
Variogram and RIC2 analysis indicated that the grid density required to capture the spatial variability was different for different years. The low RIC2 values in both fields coincide with the smaller variogram ranges observed in those years. Visual inspection of the yield maps of the two fields (Fig. 5 and 6) supports the idea that Field 1 had a less consistent spatial structure than Field 2. Trend surfaces, which followed the same patterns as Fig. 5 and 6, were very dissimilar among years for Field 1. For Field 2 trend surfaces were quite similar from 1998 to 2000. The estimates of temporal variance over the four seasons are larger than or similar to the mean spatial variances for the original and trend (large-scale) yield values in both fields. This indicates that the magnitude of both components of the spatiotemporal variance can be of equal importance in these fields. Thus, the temporal component of the spatial variability must be taken into account to make effective site-specific management decisions. Although both fields presented similar temporal/spatial ratios, there were differences between them. Field 2 presented both higher spatial and temporal variance values than Field 1 although it had an obviously persistent spatial pattern. This indicates that simple summary statistics may not capture the spatiotemporal structure of the field with sufficient descriptiveness to make them useful by themselves.
The visually persistent spatiotemporal pattern displayed by Field 2 was captured by a k-means cluster analysis. Moreover, the cluster analysis also captured a significant spatial structure in Field 1 even though this field did not display a corresponding persistent pattern. This indicates that cluster analysis may play a useful role in organizing and simplifying the complex spatiotemporal behavior of yield trends over time. This may aid in identifying the underlying causes of this behavior.
 |
ACKNOWLEDGMENTS
|
|---|
We are very grateful to the cooperating farmers, Charley Mathews and Charley Mathews, Jr., for allowing us to carry out research on their farm. We are grateful to David Clay and to an anonymous reviewer for many helpful comments. This research was supported by the California Rice Research Board, by the Instituto Nacional de Investigaciones Agropecurarias of Uruguay, and by a William F. Golden Fellowship to A. Roel.
 |
REFERENCES
|
|---|
- Bakhsh, A., D.B. Jaynes, T.S. Colvin, and R.S. Kanwar. 2000. Spatio-temporal analysis of yield variability for a cornsoybean field in Iowa. Trans. ASAE 43:3138.
- Bhatti, A.U., D.J. Mulla, F.E. Koehler, and A.H. Gurmani. 1991. Identifying and removing spatial correlation from yield experiments. Soil Sci. Soc. Am. J. 55:15231528.[Abstract/Free Full Text]
- Birrell, S.J., K.A. Sudduth, and S.C. Borgelt. 1996. Comparison of sensors and techniques for crop yield mapping. Comput. Electron. Agric. 14:215233.
- Blackmore, B.S., and C.J. Marshall. 1996. Yield mapping: Errors and algorithms. p. 403416. In P.C. Robert, R.H. Rust, and W.E. Larson (ed.) Precision agriculture. Proc. Int. Conf., 3rd, Minneapolis, MN. 2326 June 1996. ASA, CSSA, and SSSA, Madison, WI.
- Bongiovani, R., and J. Lowenberg-DeBoer. 2001. Nitrogen management in corn using site-specific crop response estimates from a spatial regression model. In P.C. Robert et al. (ed.) Precision agriculture [CD-ROM]. Proc. Int. Conf., 5th, Minneapolis, MN. 1619 July 2000. ASA, CSSA, and SSSA, Madison, WI.
- Chen, J., J.W. Hopmans, and G.E. Fogg. 1995. Sampling design for soil moisture measurements in large field trails. Soil Sci. 159:155161.
- Cliff, A.D., and J.K. Ord. 1981. Spatial processes: Models and applications. Pion, London.
- Colvin, T.S., and S. Arslan. 2001. A review of yield reconstruction and sources of errors in yield maps. In P.C. Robert et al. (ed.) Precision agriculture [CD-ROM]. Proc. Int. Conf., 5th, Minneapolis, MN. 1619 July 2000. ASA, CSSA, and SSSA, Madison, WI.
- Cressie, N.A.C. 1991. Statistics for spatial data. Wiley, New York.
- Dobermann, A., M.F. Pampolino, and H.U. Neue. 1995. Spatial and temporal variability of transplanted rice at the field scale. Agron. J. 87:712720.[Abstract/Free Full Text]
- Dobermann, A., C. Witt, D. Dawe, S. Abdulrachman, H.C. Gines, R. Nagarajan, S. Satawathananont, T.T. Son, P.S. Tan, G.H. Wang, N.V. Chien, V.T.K. Thoa, C.V. Phung, P. Stalin, P. Muthukrishnam, V. Ravi, M. Babu, S. Chatuporn, J. Sookthongsa, Q. Sun, R. Fu, G.C. Simbahan, and M.A.A. Adviento. 2002. Site-specific nutrient management for intensive rice cropping systems in Asia. Field Crops Res. 74:3766.
- Doerge, T.A. 1999. Yield map interpretation. J. Prod. Agric. 12:5461.
- Eghball, B., and J.F. Power. 1995. Fractal description of temporal yield variability of 10 crops in the United States. Agron. J. 87:152156.[Abstract/Free Full Text]
- Griffith, D.A., and L.J. Layne. 1999. A casebook for spatial statistical data analysis. Oxford University Press, Oxford, UK.
- Grondona, M.O., and N. Cressie. 1991. Using spatial considerations in the analysis of experiments. Technometrics 33:381392.
- Haneklaus, S., H. Lilienthal, E. Schnug, K. Panten, and E. Haveresch. 2001. Routines for efficient yield mapping. In P.C. Robert et al. (ed.) Precision agriculture [CD-ROM]. Proc. Int. Conf., 5th, Minneapolis, MN. 1619 July 2000. ASA, CSSA, and SSSA, Madison, WI.
- Hill, J.E., S.R. Roberts, D.M. Brandon, S.C. Scardaci, J.F. Williams, C.M. Wick, W.M. Canevari, and B.L. Weir. 1992. Rice production in California. Publ. 21498. Univ. of California Div. of Agric. and Nat. Resour., Oakland.
- Isaaks, E.H., and R.M. Srivastava. 1989. An introduction to applied geostatistics. Oxford University Press, Oxford, UK.
- Jain, A.K., and R.C. Dubes. 1984. Algorithms for clustering data. Prentice-Hall, Upper Saddle River, NJ.
- Jaynes, D.B., and T.S. Colvin. 1997. Spatiotemporal variability of corn and soybean yield. Agron. J. 89:3037.[Abstract/Free Full Text]
- Kendall, M., and J.K. Ord. 1990. Time series. Edward Arnold Publ., Kent, Great Britain.
- Lark, R.M., H.C. Bolam, T. Mayr, R.I. Bradley, R.G.O. Burton, and P.M.R. Dampney. 1999. Analysis of yield maps in support of field investigations of soil variation. p. 151161. In J.V. Stafford (ed.) Precision agriculture '99. Sheffield Academic Press, Sheffield, UK.
- Lark, R.M., and J.V. Stafford. 1997. Classification as a first step in the interpretation of temporal and spatial variability of crop yield. Ann. Appl. Biol. 130:111121.[Web of Science]
- Long, D.S. 1998. Spatial autoregression modeling of site-specific wheat yield. Geoderma 85:181197.
- Lowenberg-DeBoer, J., and K. Erickson. 2000. Precision farming profitability. Purdue Univ., West Lafayette, IN.
- Matloff, N. 1991. Statistical hypothesis testing: Problems and alternatives. Environ. Entomol. 20:12461250.
- Openshaw, S., and P. Taylor. 1981. The modifiable areal unit problem. p. 6069. In N. Wrigley and R.J. Bennett (ed.) Quantitative geography: A British view. Routledge and Kegan, London.
- Perez-Munoz, F., and T.S. Colvin. 1996. Continuous grain yield monitoring. Trans. ASAE 39:775783.
- Perez-Quezada, J.F., G.S. Pettygrove, and R.E. Plant. 2003. Spatial-temporal analysis of yield and the influence of soil factors in two four-crop-rotation fields in the Sacramento Valley, California. Agron. J. 95:676687.[Abstract/Free Full Text]
- Plant, R.E., D.S. Munk, B.A. Roberts, R.N. Vargas, D.W. Rains, R.L. Travis, and R.B. Hutmacher. 2001. Application of remote sensing to strategic cotton management. J. Cotton Sci. 5:3041.
- Porter, P.M., J.G. Lauer, D.R. Huggins, E.S. Oplinger, and R.K. Crookston. 1998. Assessing spatial and temporal variability of corn and soybean yields. J. Prod. Agric. 11:359363.
- Searcy, S.W., A.D. Beck, and J.P. Roades. 2001. Cotton yield mapping: Texas experiences in 2000. p. 342344. In Proc. 2001 Beltwide Cotton Conf., Anaheim, CA. 913 Jan. 2001. Natl. Cotton Counc., Memphis, TN.
- Searcy, S.W., J.K. Schueller, Y.H. Bae, S.C. Borgetl, and B.A. Stout. 1989. Mapping of spatially variable yield during grain combining. Trans. ASAE 32:826829.
- Swindell, J.E.G. 1997. Mapping the spatial variability in the yield potential of arable land through GIS analysis of sequential yield maps. p. 827834. In J.V. Stafford (ed.) Precision agriculture '97. Proc. Eur. Conf. on Precision Agric., 1st, Warwick Univ., UK. 810 Sept. 1997. BIOS Sci. Publ., Oxford, UK.
- Whelan, B.M., and A.B. McBratney. 2000. The "null hypothesis" of precision agriculture management. Precis. Agric. 2:265279.
- Wong, D.W.S. 1995. Aggregation effects in geo-referenced data. p. 83106. In S.L. Arlinghaus et al. (ed.) Practical handbook of spatial statistics. CRC Press, Boca Raton, FL.
This article has been cited by other articles:

|
 |

|
 |
 
S. B. Blanche and S. D. Linscombe
Stability of Rice Grain and Whole Kernel Milling Yield is Affected by Cultivar and Date of Planting
Agron. J.,
April 3, 2009;
101(3):
522 - 528.
[Abstract]
[Full Text]
[PDF]
|
 |
|

|
 |

|
 |
 
R. E. Plant
Comparison of Means of Spatial Data in Unreplicated Field Trials
Agron. J.,
March 12, 2007;
99(2):
481 - 488.
[Abstract]
[Full Text]
[PDF]
|
 |
|

|
 |

|
 |
 
A. N. Kravchenko, G. P. Robertson, K. D. Thelen, and R. R. Harwood
Management, Topographical, and Weather Effects on Spatial Variability of Crop Grain Yields
Agron. J.,
March 1, 2005;
97(2):
514 - 523.
[Abstract]
[Full Text]
[PDF]
|
 |
|

|
 |

|
 |
 
A. Roel and R. E. Plant
Factors Underlying Yield Variability in Two California Rice Fields
Agron. J.,
September 1, 2004;
96(5):
1481 - 1494.
[Abstract]
[Full Text]
[PDF]
|
 |
|