ABSTRACT

Precision Agriculture (PA) recognizes and manages intra-field spatial variability to increase profitability and reduced environmental impact. Site Specific Crop Management (SSCM), a form of PA, subdivides a cropping field into uniformly manageable zones, based on quantitative measurement of yield limiting factors. In Mediterranean environments, the spatial and temporal yield variability of rain-fed cropping system is strongly influenced by the spatial variability of Plant Available Water-holding Capacity (PAWC) and its strong interaction with temporally variable seasonal rainfall. The successful adoption of SSCM depends on the understanding of both spatial and temporal variabilities in cropping fields. Remote sensing phenological metrics provide information about the biophysical growth conditions of crops across fields. In this paper, we examine the potential of phenological metrics to assess the spatial and temporal crop growth variability across a cropping field. The study was conducted at a wheat cropping field in Minnipa, South Australia. The field was classified into three management zones using prolonged observations including soil assessment and multiple year yield data. The main analytical steps followed in this study were: calculation of the phenological metrics using time series NDVI data from Moderate Resolution Imaging Spectroradiometer (MODIS) for 15 years (2001-2015); producing spatial trend and temporal variability maps of phenological metrics; and finally, assessment of association between the spatial patterns and temporal variability of the metrics with management zones of the cropping field. The spatial trend of the seasonal peak NDVI metric showed significant association with the management zone pattern. In terms of temporal variability, Time-integrated NDVI (TINDVI) showed higher variability in the “good” zone compared with the “poor” zone. This indicates that the magnitude of the seasonal peak is more sensitive to soil related factors across the field, whereas TINDVI is more sensitive to seasonal variability. The interpretation of the association between phenological metrics and the management zone site conditions was discussed in relation to soil-climate interaction. The results demonstrate the potential of the phenological metrics to assess the spatial and temporal growth variability across cropping fields and to understand the soil-climate interaction. The approach presented in this paper provides a pathway to utilize phenological metrics for precision agricultural management application.

Keywords:

Remote Sensing, Crop Phenology, Multi-Temporal Images, NDVI, Precision Agriculture, Spatio-Temporal Variability

1. Introduction

Precision agriculture (PA) recognizes and manages intra-field spatial variability with the desired outcome of increasing profitability and reduced environmental impact [1] . Site Specific Crop Management (SSCM) is a form of PA that focuses on managing Spatial crop yield variability through matching agricultural inputs with the site potential. SSCM subdivides a cropping field into uniformly manageable zones, based on quantitative measurement of yield limiting factors [2] [3] [4] .

Spatial variability in crop yield is the result of complex interaction of factors influencing crop growth that include soil (such as nutrients, soil water availability), topographical factors (such as elevation) and climatic factors (such as rainfall and temperature) [5] . In addition to their spatial variability, some of these factors such as climate have temporal variability, which causes the spatial pattern of crop yield to vary from season to season. The optimal choice of PA over uniform farm management requires a sound understanding of such temporal variability, and is well framed by Whelan et al. [6] as a “null hypothesis” for precision agriculture i.e. ‘‘Given the large temporal variation evident in crop yield relative to the scale of a single field, then the optimal risk aversion strategy is uniform management’’. The degree of variability, and whether PA can be technically and economically beneficial to manage the variability, are the most important issues to be considered [7] . Hence, SSCM requires a comprehensive understanding of both spatial and temporal variability in crop growth.

In Mediterranean environments like South Australia, yield variability of rain-fed crops is often controlled by soil water availability. The soil property Plant Available Water-holding Capacity (PAWC) explains a high degree of intra-field spatial variability [8] [9] [10] . Hence, it is suggested as a basis for management zone delineation [8] . However, the impact of PAWC on crop yield is highly determined by the seasonal rainfall. For example, high PAWC soils in dry seasons may have small differences in yield compared with low PAWC soils, as they rarely filled to their capacity. On the other hand in years with good opening rains, but with decreasing water availability during the growing season, large differences in yield may occur between high and low PAWC soils. In such way, soil PAWC interacts with the seasonal rainfall and controls the spatial and temporal crop yield variability.

While delineation of agricultural management zones vary in terms of the information used, it is generally based on soil and other yield determining factors [5] , with crop yield data often the primary source of information. However, interpreting the temporally variable pattern of multi-year yield maps and making use of them for PA purpose is a challenging task for farmers. Researchers have developed a number of different approaches to analyse multiple year yield data. The use of spatial and temporal variograms, where the semi variance is a function of spatial lag and temporal lag [11] [12] , and development of spatial trend and temporal stability maps [13] [14] are two examples of approaches for interpretation of multiyear yield data. However, the availability of multiple year yield data, representing variable climatic conditions, has also been identified as a limitation [14] [15] [16] . As an alternative approach, a few studies have used remote sensing technology to observe spatial crop yield variability [17] [18] .

Remote sensing vegetation indices have potential to assess crop growth variability by quantifying relative growth and health condition of the crop. The Normalized Difference Vegetation Index (NDVI) is one of the most widely used indices to quantify vegetation vigour from the spectral reflectance of vegetation. NDVI has been used for many crop-monitoring applications [17] [19] [20] [21] . In precision agriculture, NDVI has been used as a surrogate for crop yield for SSCM zone delineation [17] [22] [23] . While single images are useful for yield estimation, the inter annual growth variability resulting from soil-climate interaction can produce spurious results depending on the image date selection [24] . Multi-temporal NDVI, on the other hand, provides an additional temporal dimension to uncover the vegetation dynamics. Remote sensing phenology estimates phenological growth stages including the start of season and end of season from multi-temporal vegetation index data [25] [26] . The derived metrics may not necessarily correspond directly to conventional, ground-based phenological events, but they provide important information about the vegetation growth dynamics that can be associated with environmental factors such as soil properties [e.g. [24] [25] [27] ].

The aim of this paper is to explore the potential of remote sensing derived phenological metrics for assessing spatio-temporal variability in cropping fields and to understand the soil-climate interactions that strongly influence crop yield. The specific objective of the project was to examine the relationship between spatial and temporal trend of phenological metrics and predefined management zones, in a South Australian cropping field where intra-field variability is strongly attributed to soil PAWC. The approach presented in this paper provides a pathway for future studies in utilizing phenological metrics for management zone delineation.

2. Materials and Methods

2.1. Study Area

The study was conducted in the rain-fed cropping region of South Australia. The region experiences a Mediterranean climate, with hot summers and wet winters with average annual rainfall of approximately 325 mm and average growing season (April-October) rainfall of approximately 241 mm [28] . In this region, sowing starts in late March to May, following sufficient rainfall for seeding. Following seeding, the crop germinates, grows and progressively increases in cover to reach peak greenness in September. The crops ripen and reach senescence in October and harvested in November.

The study site was a field close to Minnipa town, in Upper Eyre Peninsula (Figure 1). The field is approximately 65 ha in area, predominantly cultivated for wheat interspersed with some years of barley and pastures. This field was chosen for its soil types, which are representative of those across the wider region. It has also been studied as a focus site for Eyre Peninsula Agricultural Research Foundation (EPARF), where numbers of agricultural research and developmental trials have been undertaken [29] . The field is characterized by sandy loam to sandy clay loam soils, with the land zones of sandy rises reported to perform well in dry years, and shallow flats, which rarely perform well regardless of crop or pasture choice [30] .

2.2. Data

2.2.1. NDVI

Derivation of crop phenological information from remote sensing imagery requires a high temporal frequency of images. Moderate Resolution Imaging Spectroradiometer (MODIS) is one of the widely used sensors for phenological

studies. For this study, we used the 16 days composite NDVI image product (MOD13Q1), which has 250 m spatial resolution. The MOD13Q1 product is derived from daily NDVI data using the Constrained View Angle Maximum Value Composite algorithm which extracts the maximum NDVI value for each pixel within the 16 day interval to create a cloud-free composite image [32] [33] [34] . For this study, a total of 345 images between 01/01/2001 and 18/12/2015 were downloaded from NASA Land Processes Distributed Active Archive Center (LP DAAC) website.

Overlaying the farm boundary on the MODIS NDVI image, using ArcGIS 10.2.1 software [35] , the area intersects 17 pixels to varying degrees. In order to minimize signal contamination only the six pixels with more than 75% of the pixel area lying within the field boundary were considered (Figure 2(b)).

2.2.2. Management Zone Map

The farm at Minnipa has been a focus site for low rainfall cropping research [29] . Management zones in this field are well understood and have been delineated on the basis of a long record of scientific observations. The field has been subdivided into three management zones: good, medium and poor, using historic yield data, soil Electro Magnetic survey (EM38) and elevation maps [30] (Figure 2(a)). The good zone is characterized by red light sandy clay loam with the maximum rooting depth of 80 cm and PAWC of 108 mm. The soil type in the medium zone is red loam with a constrained maximum rooting depth of 60 cm and PAWC of 74 mm. The poor zone, in contrast, has only 40 cm maximum rooting depth with soil type of red sandy clay loam (Minnipa heavy) and PAWC of 57 mm [30] [36] .

Intersecting the MODIS pixel grids with the management zone map, the pixels were assigned continuous management zone values between 1 (poor zone) and 3 (good zone), based on their relative coverage of the zones. For example,

a pixel covered 50% by good and 50% by medium zones was assigned the zone value of 2.5, (50/100 * 3 + 50/100 * 2). The resulting management zone values for the pixels range from 2.4 to 1.43. These values may not represent quantitative measure, rather they are intended to provide ordinal ranking of relative performance between the pixels. Figure 2(b) shows the selected six pixels. The management ranks for these pixels are presented in Table 1.

2.3. Analysis

2.3.1. Overview of the Approach

Figure 3 summarizes the steps followed in our analysis. Firstly, the phenological

metrics were derived from the MODIS NDVI data, their standard scores were calculated, and then the spatial trend and temporal variance of the metrics were calculated. Finally, the correlation between the management zone values and the spatial trend (temporal mean) and temporal variability (temporal variance) of each of the metrics were assessed.

2.3.2. Extraction of Phenological Metrics

The phenological metrics were extracted using “CropPhenology” package [37] , in the R software environment [38] . Figure 4 shows representation of the phenological metrics on the NDVI growth dynamic curve. The inferred physiological descriptions of the metrics are summarized in Table 2. The PhenoMetrics function of the CropPhenology package takes the time series of MODIS NDVI imagery and the boundary shape file of the study area and provides outputs of 15 phenological metrics in raster file format.

For the values of a given pixel to be comparable across time, we normalized the values relative to the mean of the field, using the standard score. The standard score value at the given pixel indicates how the pixel value deviates from the field mean in the measure of the standard deviation. Blackmore et al. [39] have used similar standardization to compare the relative yields among different crops. Accordingly, we calculated the pixel’s standard score for each phenological metrics, using Equation (1).

$\text{standard}{\text{score}}_{ijk}=\frac{{\text{Value}}_{ijk}-{\text{Mean}}_{jk}}{{\text{sd}}_{jk}}$ (1)

where Standard score_ijk is a standard score value of metric k, at pixel i for year j, Value_ijk is the value of metric k at pixel i for a year j, mean_jk and sd_jk is the mean and standard deviation of the metric k across all pixels of year j.

2.3.3. Spatial Trend and Temporal Variability of Phenological Metrics

Spatial trend map

The spatial trend at a pixel was calculated as the average of the pixel value over the year of interest (Equation (2)). It is adapted from Blackmoreet. Al. [13] . The spatial trend values of the phenological metrics from each pixel produced data for spatial trend map of that phenological metric, which shows the average value of the metric across observation years. A total of 15 spatial trend maps were created for the 15 metrics.

${\text{Spatialtrend}}_{ik}=\frac{1}{n}\left({\displaystyle \sum _{j=2015}^{j=2001}{\text{standardscore}}_{ijk}}\right)$ (2)

where Spatial trend_ik is the mean value of phenological metric k for the pixel i, standard Score_ijk is the standard score value for phenological metric k at of pixel i on year j, and n is the total number of observation year.

Temporal variability map

The temporal variability map was created using the temporal variance overthe observation years. It is adapted from Blackmoreet. Al. [13] . Temporal variability is classified into two categories: inter-year variability and relative temporal variability [13] . Inter-year variability is the temporal variability caused by the change in annual rainfall, which affects the overall productivity of the field. The second category described how each part of the field behaved relative to the other parts from year to year. In this analysis, we are interested in the second category, which assesses the relative variability of growth at a particular pixel, over the study period. This relative temporal variability was calculated using a temporal variance, a modified variance measured in a function of time [6] [13] . The temporal variance of all the metrics were calculated for each pixel across years, using Equation (3). And the temporal variance value of the metric from each pixels produced the temporal variability map of the phenological metric. Accordingly, 15 temporal variability maps were created for the 15 metrics. The temporal variance of each metric across the season revealed the temporal variability of the growth conditions at the various phenological stages. The calculated temporal variance has low values if the pixel is stable over time and often has values close to the mean of the observation seasons. On the other hand, if the metric is temporally unstable, with high values in some years and low values in others, its temporal variance will be high.

$${\text{Temporal Variance}}_{ik}=\frac{1}{n-1}\left({\displaystyle \sum _{j=2015}^{j=2001}{\text{standard score}}_{ijk}}-{\text{Spatial trend}}_i\right)$$ (3)

where Temporal Variance_ik is the temporal variance at pixel i of metric k, standard score_ijk is the standard score of metric k at pixel i for year j, Spatial trend_i is the Spatial trend of pixel i and n is the number of observation year.

2.3.4. Relationship between Management Zone and Trends of Phenological Metrics

In this study, the monotonic relationship between the management zone values and the trends of the 15 phenological metrics were assessed. Spearman rank correlation was used to test the relative direction and strength of the relationship between the management zone pixel values and phenological metrics, in R computing environment. Spearman correlation ranks the variables and provides the correlation coefficient, rho (ρ) that indicates the strength of the correlation [40] .

3. Results and Discussion

3.1. Relationship between Management Zone and Spatial trend of Phenological Metrics

The Spearman rank correlations between management zone and the spatial trend of phenological metrics are shown in Table 3. In this analysis, we have undertaken multiple tests on a single dataset. This can introduce the probability of making type I errors, family wise error rate (FWER). To address this problem, we used the Bonferroni correction that adjusts the significance level (p value) [41] , by dividing the p value by the number of comparisons made, which results a new p value of 0.0333.

The results indicate that the spatial trend of Max Value is significantly

correlated with the management zone values (ρ(13) = 0.82, p = 0.024). Although the associations are not significant, the GreenUpSlope and MaxT showed high correlation coefficients, with different directions of association. The Asymmetry, TINDVIBeforeMax and BeforeMaxT metrics also showed moderate association with the management zone values.

This result corroborates previous observations that have recognized the ability of NDVI to estimate crop yield [42] [43] [44] . Although strong associations between NDVI and yield have been reported at differing stages during the growth period, these dates are mostly around heading and flowering stages [17] [19] [20] [21] [45] [46] [47] [48] . As crops often attain their maximum NDVI close to heading and flowering stages [49] , the high association of MaxValue with the management zone can be directly linked to these observations.

The results also show trends of increasing Max Time and decreasing of GreenUpSlope, moving from the poor to good management zone. As the soil PAWC strongly relates with yield potential in South Australia, the observed relationship is likely be related with water availability. Numbers of researchers reported that water stress has a strong influence on crop growth stages with evident observation of faster growth to the next growth stages in water stressed crops. Under similar climatic conditions, plants in low PAWC soils will likely to face water stress between rainfall events, as there is less water stored in the soil profile compared with the high PAWC soils. The early Maximum (low MaxT) can then be related to the hasten growth of the crop to anthesis due to water stress in the poor zone [25] [50] [51] [52] . This is also reflected as faster rate of greenness, higher GreenUpSlope, in low PAWC soil.

Similar observations of low Max Time and high GreenUpSlope were observed in our previous study [25] . Unlike our current approach of averaging values across the years, the previous study considered paired comparison of pixels across number of years. The pared comparison study showed higher Max Values in low PAWC soils than the high PAWC soils. Such pattern of paired comparison observed in the previous study may average out thorough the observation years, in the current study, to reflect average higher Max Value for higher PAWC soils.

The positive moderate association of Asymmetry, TINDVIBeforeMax and BeforeMaxT metrics, indicated that the good zones have higher values of these metrics than the poor zone. Generally, this observation shows that the good zone experienced longer time before the NDVI peak, higher area under the curve between the Onset and Maximum NDVI, and higher difference between TINDVIBeforeMax and TINDVIAfetrMax than poor zone. These differences are related with the high and late NDVI peak in the good zone than the poor zone, which is explained above. Furthermore, the trend of TINDVIBeforeMax in our analysis agrees with pevious field studies that confirm the strong association of pre-anthesis growth with yield potential, as it is the time when the number of grains are determined [53] .

3.2. Relationship between Management Zone and Temporal Variability of Phenological Metrics

The seasons considered in this study, 2001-2015, were characterized by variable rainfall amount and seasonality. The temporal variance of the phenological metrics measures the relative variability of the metrics across the variable rainfall. The Spearman correlation analysis assessed the trend of temporal variability of phenological metrics across management zones (Table 4).

The results indicate that the temporal variability of TINDVI has a significant positive correlation with the management zone, indicating an increase of temporal variability from poor to good management zones. This variability demonstrates the interaction of the spatial soil variability with the temporal rainfall seasonality. In good seasons, the spatial variability of soil PAWC becomes significant, as the high PAWC soils store more water than the low PAWC soils. However, in dry seasons there will be less difference in stored water between the contrasting soils. This results in larger spatial variability in good seasons compared with dry seasons [8] [54] . Temporally, the good zone with high PAWC soil will have higher variability across variable seasons than the poor zone with low PAWC soils [54] .

Previous research has reported that TINDVI can be a good indicator for crop production [20] [44] [55] . However, the sensitivity of TINDVI for crop production depends on rainfall seasonality; it can be a good indicator for crop production when water is the limiting factor for crop growth [55] . Considering the profound correlation between TINDVI and yield in cropping fields, the high temporal variability in TINDVI in good zone than in poor zone coincide well with the yield data observation in the study area. Yield data from EPARF research shows higher difference in good zone (3.1 t/ha) between the dry year 2008 and the wet year 2010 than the poor zone (2.13 t/ha) [30] [56] .

The relationship between the temporal variability of MaxV and the management zone showed weak correlation. We have observed in the previous section that the spatial trend (temporal mean) of this metrics showed a strong correlation with management zone values, when the low correlation of temporal variability with the management zone can be interpreted as both the good and poor zones are variable to a similar degree with the higher mean values from good zone than the poor zone. Thus, the MaxV metrics is more sensitive to spatial yield variability than temporal variability; on the other hand, TINDVI better reflects the temporal yield variability.

3.3. Summary of Indicative Phenological Metrics for Crop Field Management

Generally, our analysis demonstrates the potential of phenological metrics to assess the intra-field crop growth variability and to provide a comprehensive understanding of soil-climate interactions. Single vegetation index images have been successfully utilized to recognize yield variability across a crop field [e.g. [17] ]. However, the single image approach has been criticized for lacking information on intra-seasonal growth dynamics [24] . On the other hand, multi- temporal images and the derived phenological metrics uncover the intra-annual biophysical properties of the crop across the field, as driven by soil-climate interaction. This potentially provides a better understanding of both inter-annual and intra-annual variability, which is a key factor to improve management practices.

Furthermore, the method presented in this analysis show the efficacy of phenological metrics to recognize crop growth variability. Considering the increasing availability of remote sensing imagery, the spatio-temporal variability estimation using phenological metrics can provide valuable information for PA suitability assessment in areas where there is limited availability of yield data.

Currently, MODIS is the most appropriate imagery for such analysis due to its high temporal resolution. However, its coarse spatial resolution limits the application of this method to broad acre cropping. The technology of integrating the higher spatial resolution of sensors such as Landsat or Sentinel with the temporal resolution of MODIS to generate high temporal/high spatial resolution datasets [57] may allow future application for smaller fields. Whilst these fused spatio-temporal datasets are presently not available to the public at broad extents, once available, may be utilized at a scale typical for yield mapping in precision agriculture today. In such a way, this paper provides a basis for further analysis to utilize phenological metrics for crop management support systems.

4. Conclusions

In this study, the spatial and temporal variability of crop growth was assessed using remote sensing phenological metrics in relation to the relative yield potential of sites in the field. Phenological metrics provide information about the temporal and spatial variability of plant growth across cropping fields. Our study demonstrates the potential of satellite based phenological metrics to provide information about growth variability in cropping fields, which provide relevant information for agricultural management. The results of our analysis show that Time-integrated NDVI reflects seasonal effect on crop growth, whereas the magnitude of NDVI peak strongly reflects the soil quality and showed spatial variability of long-term site conditions across the field.

Crop yield response to PAWC is strongly influenced by the amount and seasonal distribution of rainfall, which is temporally variable. The effect is stronger when there is sufficient rainfall to fill the soil profile at the start of the season and the plants are dependent on deep stored soil-water as the season progresses. Phenological metrics provide comprehensive insight of the spatio-temporal crop growth variability across the variable seasons, which advances our understanding of soil-climate interaction across the cropping field.

The method presented here provides a pathway towards better estimation of spatio-temporal crop growth variability, which is vital for PA success. It can be used to develop future, more detailed studies to fully utilize the potential of phenological indicators for site characterization through unravelling the complex spatio-temporal soil-climate interactions.

Acknowledgements

The Authors wish to thank the Eyre Peninsula Agricultural Research Foundation particularly Linden Masters for data support and the Faculty of Sciences, University of Adelaide for research funding.

Cite this paper

Araya, S., Ostendorf, B., Lyle, G. and Lewis, M. (2017) Remote Sensing Derived Phenological Metrics to Assess the Spatio-Temporal Growth Variability in Cropping Fields. Advances in Remote Sensing, 6, 212-228. https://doi.org/10.4236/ars.2017.63016

References