Estimating spatio-temporal patterns of agricultural productivity in fragmented landscapes using AVHRR NDVI time series

Remote Sensing of Environment 84 (2003) 367 – 384 www.elsevier.com/locate/rse

Estimating spatio-temporal patterns of agricultural productivity in fragmented landscapes using AVHRR NDVI time series Michael J. Hill a,*, Graham E. Donald b a

b

Bureau of Rural Sciences, PO Box E11, Kingston, ACT 2604, Australia CSIRO Livestock Industries, Leeuwin Centre for Remote Sensing Technologies, Floreat Park, WA 6014, Australia Received 2 May 2002; received in revised form 16 August 2002; accepted 24 August 2002

Abstract The characteristics of Normalized Difference Vegetation Index (NDVI) time series can be disaggregated into a set of quantitative metrics that may be used to derive information about vegetation phenology and land cover. In this paper, we examine the patterns observed in metrics calculated for a time series of 8 years over the southwest of Western Australia—an important crop and animal production area of Australia. Four analytical approaches were used; calculation of temporal mean and standard deviation layers for selected metrics showing significant spatial variability; classification based on temporal and spatial patterns of key NDVI metrics; metrics were analyzed for eight areas typical of climatic and production systems across the agricultural zone; and relationships between total production and productivity measured by dry sheep equivalents were developed with time integrated NDVI (TINDVI). Two metrics showed clear spatial patterns; the season duration based on the smooth curve produced seven zones based on increasing length of growing season; and TINDVI provided a set of classes characterized by differences in overall magnitude of response, and differences in response in particular years. Frequency histograms of TINDVI could be grouped on the basis of a simple shape classification: tall and narrow with high, medium or low mean indicating most land is responsive agricultural cover with uniform seasonal conditions; broad and short indicating that land is of mixed cover type or seasonal conditions are not spatially uniform. TINDVI showed a relationship to agricultural productivity that is dependent on the extent to which crop or total agricultural production was directly reduced by rainfall deficiency. TINDVI proved most sensitive to crop productivity for Statistical Local Areas (SLAs) having rainfall less than 600 mm, and in years when rainfall and crop production were highly correlated. It is concluded that metrics from standardized NDVI time series could be routinely and transparently used for retrospective assessment of seasonal conditions and changes in vegetation responses and cover. D 2002 Elsevier Science Inc. All rights reserved. Keywords: NDVI; Time series; Metrics; Agriculture; Spatial patterns; Temporal variation

1. Introduction Time series of National Oceanic and Atmospheric Administration Advanced Very High Resolution Radiometer (NOAA AVHRR) Normalized Difference Vegetation Index (NDVI) have been subjected to a wide variety of analytical approaches. There are two major categories for analysis; signal decomposition to extract major seasonal and interannual periodic responses (e.g. Li & Kafatos, 2000; Moody & Johnson, 2001; Roderick, Noble, & Cridland, 1999); and quantification of NDVI curve properties and characteristics in order to derive information about plant and land cover

* Corresponding author. E-mail address: [email protected] (M.J. Hill).

phenology (Reed et al., 1994; Reed, Swets, Bard, Brown, & Rowland, 2001). In each of these cases, spatial patterns of interannual or phenological responses are also of interest. NDVI time series may be analyzed to generate a set of metrics that summarize the phenology of vegetation (Lloyd, 1990; Malingreau, 1986; Reed et al., 1994). These metrics break the NDVI curve down into measures of the timing and magnitude of NDVI response. They have been used in the USA for the analysis of grassland land cover classes on the Prairies (Reed, Loveland, & Tieszen, 1996; Tieszen, Reed, Bliss, Wylie, & DeJong, 1997; Yang, Wylie, Tieszen, & Reed, 1998). One of the major metrics usually estimated is time integrated NDVI (TINDVI). Early research on monitoring and measuring rangeland biomass in Sahelian Africa used integrated NDVI as a predictor of end of season biomass (Diallo, Diouf, Hanan, Ndiaye, & Pre’vost, 1991; Prince, 1991).

0034-4257/02/$ - see front matter D 2002 Elsevier Science Inc. All rights reserved. PII: S 0 0 3 4 - 4 2 5 7 ( 0 2 ) 0 0 1 2 8 - 1

368

M.J. Hill, G.E. Donald / Remote Sensing of Environment 84 (2003) 367–384

Metrics represent one of a number of approaches to decomposition of NDVI time series. The metrics decompose the curve into a set of statistics—reducing the curve to its parts within an individual cycle. However, time series may also be unmixed or decomposed to separate temporal oscillations of different periodicities. These can involve traditional time series analysis to extract trend and seasonal components to identify the different greening rhythms of trees and grassland in a mixed woodland landscape (Roderick et al., 1999); using phase and amplitude of different harmonics from a Fourier transform to describe seasonal properties (Moody & Johnson, 2001); or using principal component analysis and wavelet transforms to extract from vegetation signatures the underlying climatic oscillations such as those driven by the El Nin˜o and La Nin˜a cycle (Li & Kafatos, 2000). With the advent of the TERRA satellite, the MODerate resolution Imaging Spectroradiometer (MODIS) sensor can provide many derived products for quantification of terrestrial responses such as leaf area index, net primary productivity and fraction of absorbed photosynthetically active radiation (Salomonson, Guenther, & Masuoka, 2001). However, there is still an essential role for NOAA AVHRR NDVI, and MODIS NDVI and Enhanced Vegetation Index (EVI; Huete, Liu, Batchily, & van Leeuwen, 1997) time series in description of the temporal and spatial patterns of vegetation behavior in response to climate and other agents of change (e.g. Hill, Vickery, Furnival, & Donald, 1999). Bimonthly composites of the NDVI at 1.1 km resolution have been constructed from NOAA AVHRR data for the whole of Western Australia (WA) since 1991. Such data are widely used in vegetation monitoring and have been the subject of substantial Australian research on calibration, and rangeland and agricultural monitoring (Cridland, Burnside, & Smith, 1994; Roderick, Smith, & Lodwick, 1997a,b; Smith, 1994; Smith, Adams, & Hick, 1994). In WA, the agro-climatic system is highly seasonal with winter rain, summer drought and a growing season defined by moisture availability for plant growth and determinate, obligate annual plant forms. Within this system, NDVI may be used to characterize the temporal extent of the growing season and productive potential of agricultural land on a regional basis (Roderick et al., 1997a,b). Initial analysis of NDVI metrics for the south west of WA found that on a Statistical Local Area basis, TINDVI was well correlated with wool production per hectare, but the correlation with wool production per sheep was poor (Hill, Donald, & Smith, 1998). However, the correlation really reflected the influence of long-term climate patterns on geographical distribution of sheep, and hence the resulting productivity of wool enterprises. Subsequently, we have improved the processing of the time series and adjusted the procedure for calculating NDVI metrics. NDVI metrics could provide the means to describe the seasonal pattern of agricultural production in some detail, and provide for retrospective analysis of trends in land use and production patterns. In this paper, we

use NDVI metrics to explore variation in seasonal vegetation response patterns for southwestern WA, and examine the relationship between spatial and temporal patterns of NDVIbased metrics and measures of agricultural land use and production in the south west of WA.

2. Methods 2.1. Data pre-processing Bimonthly composites of NOAA AVHRR NDVI were obtained from the WA Department of Land Administration for the period 1991 – 1999. Composites were calculated using a maximum value method (Holben, 1986). Data were supplied as integer images for each biweekly period scaled between 0 and 1000, but contained some maximum or minimum integer values associated with processing errors or clouds. The images were truncated to the range 0– 1000. Large negative or positive values were assigned to 1 or 1001, respectively. Cells with these values were replaced with cloud-free values from neighboring pixels using a linear interpolation of values over time between available data points where more than one time interval was missing. Cells with anomalously high NDVI values between 800 and 999, and anomalously low NDVI values between 1 and 100 were also replaced by the same method. Thresholds for high and low limits were set after examination of average maximum and minimum NDVI values for NDVI profiles not affected by data spikes or troughs. Extreme spikes or troughs within the normal data range of 100– 750 were filtered on the basis of exceeding a threshold percentage deviation from preceding and following data points. Time series were standardized to remove differences due to changes in satellite and sensor drift for each satellite using trends for regional mean data (Roderick, Smith, & Lodwick, 1996). 2.2. Agricultural production data The Australian Bureau of Statistics collects agricultural census data on a regular basis. The major unit for aggregation and reporting is the Statistical Local Area (SLA) which corresponds closely but not exactly to local government area boundaries, and which may change slightly between census dates particularly surrounding large rural towns. Agricultural production data were extracted from a long-term census data set for 1983 –1997 that was concorded for a single set of statistical local area boundaries. Data are reported at 31 March for each census year, and therefore statistics actually relate to production condition in the previous year, i.e. 1997 census data refer to production for the 1996 growing season and are shown as such. We chose four crops and three animal production systems to represent agricultural productivity for the Western Australian agricultural zone: wheat, barley, oats, lupins, beef cattle, sheep and dairy cows. Other agricultural production activities were

M.J. Hill, G.E. Donald / Remote Sensing of Environment 84 (2003) 367–384

either negligible for the study period, or occur, like grape growing, at a small spatial scale unsuited to analysis with AVHRR data. Agricultural production is represented by a variety of commodities. We needed to represent production from livestock and cropping by a common currency. A National Land and Water Resources Audit project (Walcott, personal communication) addressed this issue and chose dry sheep equivalents (DSE) to represent productivity. Productivity conversion factors were based on calculated metabolisable energy (ME) per kg of dry matter for each of the crops. These were then converted to DSEs based on a DSE consuming 2550 MJ ME per year using the following factors: wheat, 5.0 dse/t; barley, 4.7 dse/t; oats, 4.1 dse/t; and lupins, 4.7 dse/t. Cattle numbers were converted to DSE on the basis of 8 dse/beef animal and 12 dse/dairy cow. Productivity for the study region was therefore represented by the sum of the animal numbers converted to dse and the crop production in tones converted to dse using the above factors. 2.3. Calculation of NDVI metrics NDVI metrics were calculated using the method described by Reed et al. (1994). The time series was smoothed using a 5 interval running mean smoothing algorithm implemented in the Interactive Data Language (Research Systems Inc.) (Fig. 1). This provided the basic curve for metric estimation. In order to define the beginning and end of the period when NDVI is significantly greater than its minimum value, forward and backward lag curves are calculated. The period of lag should correspond approximately to the length of the non-growth period for the environment in question. Examination of historical seasonal

369

patterns for the agricultural zone of WA showed that in general the growing season commenced by the middle of May and concluded by the end of October, a period of 11 half-months. On the basis of this, the forwards and backwards lag periods were set to 13 intervals corresponding to non-growing period of 6.5 months (Reed et al., 1994; Fig. 1). The phenological metrics were calculated as shown in Fig. 2 with crossover points between the smoothed time series and lagged running means used to denote the beginning and end of the NDVI ‘‘green’’ period. In the first instance, we calculated the same basic set of metrics (Table 1) as described by Reed et al. (1994). However, we also devised some additional supplementary metrics to provide more sensitivity, and better seasonal definition for the environment of the Western Australian agricultural zone (Table 1). These were based on (a) the maximum and minimum values and times for the smooth curves (SMMaxV, SMMinV, SMMaxT, SMMinT); and (b) the value and time at which half of the range was achieved on the upwards and downwards path (HRanVO = half range value at onset, HRanVE = half range value at end, HRanTO = time of half range value at onset, HranTE = time of half range value at end); and (c) the time period between HranTO and HRanTE (HDurT). 2.4. Analysis We chose four approaches to analysis. (1) Characterization of spatial patterns in key NDVI metrics across the agricultural zone of WA. Spatial patterns of key metrics were characterized by calculating temporal mean and standard deviation layers for selected metrics showing significant spatial patterns. (2) Unsupervised classification of the agricultural zone based on temporal and spatial patterns of key NDVI metrics.

Fig. 1. Data for 1996 – 1998 for the Kojonup site giving a graphical depiction of the method for estimation of seasonal metrics. The graph shows the corrected raw NDVI, smoothed NDVI, and forward and backward lag curves used to define the beginning and end of the NDVI ‘‘green period’’. The bars indicate bimonthly rainfall for the same period showing the seasonal climate.

370

M.J. Hill, G.E. Donald / Remote Sensing of Environment 84 (2003) 367–384

Fig. 2. A diagram showing the derivation of the basic metrics described by Reed et al. (1994) as attributes of an NDVI profile and additional metrics associated with the smooth curve. In this figure, the concept is a modification of that published by Reed et al. (1994) in the Journal of Vegetation Science, however the curve is derived from the data used in this study, and the drawing is completely new.

Three specific metrics showed systematic spatial variation which indicated that temporal patterns may also be specific and diagnostic; TINDVI, HDurT and HRanTO (Table 1). Metric layers for 1992– 1998 were subjected to a clustering procedure followed by an unsupervised max-

imum likelihood classification in the ArcInfo GIS. Classes were aggregated to final groups using a dendrogram and relative distance statistics. Mean and standard deviations for each metric were extracted for each TINDVI temporal class.

Table 1 Description of NDVI metrics calculated for the agricultural zone of WA for the period 1991 – 1998 Abbreviation

Definition

Metric

Basic metrics OnT OnV EndT EndV MaxT MaxV DurT RanV RIN RDN TINDVI DurNT

Intersection of forward lag and smooth curve Value of NDVI at forwards intersection Intersection of backwards lag and smooth curve Value of NDVI at backwards intersection Time of maximum raw corrected NDVI Maximum value of corrected raw NDVI Time from forwards to backwards intersections Difference between minimum and maximum value of smooth curve Slope of line from forwards intersection to raw maximum Slope of line from raw maximum value to backwards intersection Integrated area under smooth NDVI curve Time from backwards to forwards intersection

Starting date of NDVI high period NDVI at start of high period End date of NDVI high period NDVI at end of high period Date of maximum NDVI Maximum NDVI Length of NDVI high period Amplitude of season Rate of NDVI increase Rate of NDVI decrease ‘‘Magnitude’’ of season Length of NDVI low period

Supplementary metrics RRINDN HRanTO HRanVO HRanTE HRanVE HDurT SMMaxT SMMaxV SMMinT SMMinV

Rate of increase/rate of decrease Time of half range value at onset—equals OnV+(RanV/2) when rising Half range value at onset—OnV+(RanV/2) Time of half range value at end—equals EndV+(RanV/2) when falling Half range value at end—EndV+(RanV/2) Duration of period from HRanTO to HRanTE Time of maximum smooth NDVI curve Maximum value of smooth NDVI curve Time of minimum smooth NDVI curve Minimum value of smooth NDVI curve

‘‘Quality’’ of season Start of active growing season NDVI at start of active growing season End of active growing season NDVI at end of active growing season Duration of active growing season Date of peak of season Value at peak of season Date of season minimum Value of season minimum

M.J. Hill, G.E. Donald / Remote Sensing of Environment 84 (2003) 367–384

(3) Analysis of NDVI metric properties of eight SLAs typical of climatic and production systems across the agricultural zone. The SLA is a useful common unit for relating agricultural production to land use and process information derived from remote sensing. We chose eight SLAs across the agricultural zone to represent the spatial variation in land use, predominant agricultural enterprise and climate (Fig. 3). The SLAs represent the northern coastal area (Northampton); the northern wheatbelt (Moora); the central wheatbelt (Merredin); the sheep/wheat zone (Brookton); the core wool zone (Kojonup); the dairy zone (Bussleton); the marginal-cropping zone (Jerramungup) and the unique environment of Esperance. The latter region with its south-facing coast line, and large distance from the Indian Ocean influence, experiences different climatic patterns to all of the other areas—it has a very rapid decline in average winter rainfall with distance from the south coast, and a tendency to receive significant late summer rainfall from decaying northern monsoon systems. Fig. 4 provides a characterization of the agricultural properties for these eight SLAs in terms of land use in 1994 and temporal trends between 1992 and 1997 in numbers of beef cattle and sheep, and areas of wheat and pasture sown. The data show a decline in sheep numbers and a rise in cattle numbers and area of wheat sown over the period. The SLAs show a mixture of land uses with remnant vegetation, sown pasture, cereal and legume (lupins) crops. However, extremes are represented by Merredin, predominantly cereal cropping; Bus-

371

selton, with cattle grazing only; and Northampton, where almost 70% of the area is non-agricultural. We hypothesized that metrical properties of these SLAs would be distinctive and able to be related to the agricultural production systems, land use and climate. (4) Examination of the relationship between TINDVI and agricultural productivity across the agricultural zone of WA. A long history of research has shown that TINDVI can be related to seasonal and/or end of season production (e.g. Prince & Tucker, 1986). We hypothesized that TINDVI would be an indicator of overall agricultural productivity particularly where productivity was highly dependent upon receipt of a minimum seasonal rainfall, and where this was not always received. The analysis looked at total production in DSEs, and productivity as DSE/unit area. The procedure involved forming regression relationships between production and area and the multiple of area and TINDVI, and between total and crop productivity in dse/unit area and average TINDVI per SLA for each year from 1992 to 1996, since census data were only available to coincide with seasonal metrics for this period.

3. Results 3.1. Spatial patterns An overview of the seasonal changes in spatial patterns of landscape greening can be seen in the time series of

Fig. 3. Map showing statistical local areas for the agricultural zone identifying eight SLAs selected to sample the temporal and spatial patterns of NDVI metrics.

372

M.J. Hill, G.E. Donald / Remote Sensing of Environment 84 (2003) 367–384

Fig. 4. Agricultural statistics, showing change between 1992 and 1997 in sheep and beef cattle numbers (1995 – 1997 only) and area of wheat and sown pasture; and land use characteristics based on the year 1994 for the eight sample SLAs showing the proportion allocated to different cropping systems and pasture.

M.J. Hill, G.E. Donald / Remote Sensing of Environment 84 (2003) 367–384

Fig. 5. Map of the agricultural zone of WA showing the variation in distribution and magnitude of TINDVI for growing seasons from 1992 to 1999.

373

374

M.J. Hill, G.E. Donald / Remote Sensing of Environment 84 (2003) 367–384

TINDVI for the southwest of WA for the years 1992 to 1999 (Fig. 5). These maps showed a high degree of spatial variation in TINDVI between years with 1993 standing out as a good season, and 1996 as a poor season. Data for 1999 were affected by cloud in the far eastern area and were not included in subsequent analysis. This interannual spatial variation in TINDVI is an indicator of the general regional vegetation productivity, since NDVI is a good surrogate for net primary productivity in landscape with relatively uniform plant functional types. However, further analysis is required to determine whether these patterns translate to an indicator of agricultural production. Several time-based metrics can provide insights into the temporal properties of the growing season. The end of the season in this area tends to result from an interaction between increasing temperature, reducing rainfall and determinate reproductive behavior from obligate annual pasture and crop plants. There is therefore less variation in the metrics associated with the end of the season. However, season start date (HRanTO), time of maximum smooth NDVI (SMMaxT) and duration of the season (HDurT) show distinct patterns (Fig. 6). The time of maximum smooth NDVI shows relatively early occurrence in the far north grading to a much later occurrence in the south western corner with relatively low temporal variability except along the western and southern coastline where cloud and coastal rainfall early and late in the season lead to greater variability in NDVI signal. This latitudinal trend is most likely related to the later onset of winter and higher winter temperatures in the northern areas. The starting date metric shows a clear demarcation between the coast and southwestern corner and the central inland. This boundary most probably relates to the transition from pastures to cropping as the predominant land use. This transition coincides with a rapid decline in seasonal rainfall. Interannual variability in this starting date metric is much higher in southern areas where rainfall sufficient to initiate the start of season may occur very early or be delayed until a more general break of season. By contrast, in the northern areas, interannual variability is much lower as significant rainfall does not usually penetrate into these areas until the winter patterns of frontal systems are well established. Season duration shows a similar pattern to starting date—this is expected since there is limited variation in ending date. The far eastern areas show slightly longer seasons and earlier starting dates than the central inland areas—this phenomenon reflects the influence of late summer rain from decaying cyclonic depressions crossing from the Indian Ocean to the Southern Ocean as discussed earlier, and some changes in soil type. 3.2. Classification Three metrics provided a sound basis of classification to provide generalized maps describing key seasonal properties of the agricultural zone. The timing of start of

season (HRanTO) and season duration based on the smooth curve (HDurT) provide relatively simple classifications with similar patterns (Fig. 7a,b). The zones of HRanTO from 1 to 10 represent increasing lateness in break of season in days from January 1 to 129.5 F 20.1 (insignificant in area); 135.6 F 8.7; 144.1 F 6.7; 150.9 F 7.7; 158.1 F 6.8; 164.0 F 9.9; 166.2 F 7.5; 172.9 F 6.1; 179.3 F 6.8; and 186.8 F 6.2. The starting date is almost entirely a function of the timing of significant rainfall in autumn reflecting the degree penetration of frontal rain inland. However, it is also affected by the lag in sowing and emergence of crops relative to emergence of pastures from soil seed stores. The zones of HDurT from 1 to 7 represent increasing average length of growing season in d a ys 1 2 7. 4 F 1 2 .0 ; 1 44 . 6 F 1 6. 8 ; 14 4 .5 F 15 . 4; 163.3 F 20.4; 177.2 F 34.1; 179.4 F 28.1; and 218.3 F 32.9. The variability increases with length of growing season, proximity to the coast and reduction in cropping as a land use, since crops provide phenological limits to the length of growing season. A classification based on TINDVI (Fig. 7c) reveals a set of classes that are characterized by differences in overall magnitude and differences in response in particular years (Table 2; Fig. 9). This latter property reflects the spatial variation depicted in Fig. 4. When averaged over all years, the classes are defined by increasing value of TINDVI with relatively similar spatial variability (data not given), but classes with similar average TINDVI are differentiated by degree of variability in temporal values (Table 2). The classes with their mix of different magnitudes and interannual variation in TINDVI (Fig. 7; Table 2) identify distinct interactions between landscape characteristics  production response, i.e. the mix of crops and pastures, the types of crops and pastures and the climate-induced patterns of production (Table 2). The NDVI metrics for these classes represent the components of the NDVI curve that go to make up the interannual and spatial pattern of TINDVI (Fig. 8). An examination of the average values of the other metrics for the TINDVI classes reveals both the major factors responsible for the TINDVI response and the limits and benefits of particular metrics. Firstly, in general, the time of maximum actual and smooth NDVI tends to be later, and the value of maximum NDVI tends to be higher for classes with high TINDVI than for classes with lower TINDVI. The starting and ending times and values for the season, and the rates of increase and decrease in NDVI vary with landscape properties, with the behavior associated with the eastern margins of the wheatbelt, classes 1, 2 and 3, providing an interesting study (Fig. 8). OnT is much the same for the three classes, however EndT becomes slightly later from class 1 to class 3. These metrics are not very sensitive to seasonal variation within a region having relatively similar response patterns, since they represent the time at which almost all trace of NDVI response above the minimum has disappeared. These met-

M.J. Hill, G.E. Donald / Remote Sensing of Environment 84 (2003) 367–384 Fig. 6. Means and standard deviations of metrics describing the temporal properties of the growing season. (a) Time of maximum smooth NDVI curve (SMMaxT)—date of peak of season; (b) standard deviation of SMMaxT; (c) time of half range value at onset (HranTO), —start of active growing season; (d) standard deviation of HranTO; (e) duration of period from HRanTO to HranTE (HdurT)—duration of active growing season; and (f) standard deviation of HDurT. 375

376

M.J. Hill, G.E. Donald / Remote Sensing of Environment 84 (2003) 367–384

Fig. 7. Map of the agricultural zone of WA showing a classification based on (a) TINDVI, (b) HDurT—duration of season in days and (c) HRanTO—start of season in days. Legends show class numbers and mean values for class.

Table 2 Class descriptions for the TINDVI classification (Fig. 7a) 1992

1993

1994

1995

1996

1997

1998

1999

Description

1

2280 F 1121

2633 F 1028

1916 F 657

1492 F 60

1384 F 516

1889 F 643

1723 F 662

1905 F 1232

2

2580 F 730

3306 F 579

2600 F 678

2461 F 517

2205 F 506

2656 F 621

2764 F 603

2880 F 958

3

2816 F 446

3697 F 428

2802 F 389

2777 F 376

2898 F 411

2788 F 480

3644 F 407

3445 F 530

4

3512 F 669

4174 F 694

1776 F 501

2542 F 642

2857 F 768

3243 F 585

3093 F 577

1859 F 1246

6 7

3480 F 494 2629 F 460

4116 F 427 3786 F 661

3536 F 400 4057 F 399

3329 F 420 3753 F 400

3124 F 524 2730 F 426

3584 F 424 3297 F 443

3510 F 439 3696 F 577

3495 F 845 3305 F 839

8

3080 F 522

3635 F 530

3497 F 453

3017 F 392

2032 F 438

4055 F 427

3635 F 453

3005 F 1019

9

3422 F 559

4186 F 612

4450 F 336

4261 F 344

3210 F 494

3975 F 481

4218 F 475

3887 F 545

5

2493 F 901

975 F 771

3560 F 929

3340 F 749

2500 F 726

3087 F 760

3318 F 982

2767 F 1651

10

2483 F 1046

NA

2752 F 939

2635 F 910

2273 F 848

2888 F 906

2874 F 1060

2654 F 1367

11

2623 F 1283

3552 F 2241

2843 F 981

2200 F 1092

1591 F 894

2533 F 1145

NA

2373 F 1823

Represents areas with consistently low TINDVI associated with rivers and streams, and mixtures of remnant vegetation and agricultural land throughout the area. Identifies the eastern margins of the central and northern wheatbelt with a consistent but moderately low TINDVI. Represents a transition zone between the eastern margins and the main wheat belt, with a consistently higher TINDVI than class 2. Identifies the majority of the non-coastal areas in Esperance, and a small area on the eastern edge of the main agricultural zone with high variability in TINDVI between years. Identifies the northern wheatbelt. Represents interannually variable areas in the core wool-growing region. Identifies the central and southern wheatbelt characterized by a large decline in TINDVI in the poor year of 1996 amongst generally high values for the other years. Represents the highest TINDVI characterizing the prime wool growing areas and some perennial pasture on the coastal sand plain north of Perth. The interannual pattern is also evident here. This class is insignificant and based on very low TINDVI for 1993 due to edge effects with heavily treed areas in the far south. Affected by very high values derived from cloud interference for 1993. Affected by very high values derived from cloud interference for 1998.

M.J. Hill, G.E. Donald / Remote Sensing of Environment 84 (2003) 367–384

Class

377

378

M.J. Hill, G.E. Donald / Remote Sensing of Environment 84 (2003) 367–384

Fig. 8. Histograms showing the means and standard deviations of 15 metrics for the main TINDVI classes (Fig. 5; Table 2). Insignificant classes 5, 10 and 11 are excluded for clarity. Bars indicate the temporal standard deviation.

rics may be more sensitive at continental or global scale where major differences in seasonal patterns are evident, such as changes from winter to summer growth (Malingreau, 1986). OnV and EndV increase from class 3 to class 1 reflecting the presence of more mixed vegetation at the edge of the agricultural zone and the effect of late summer rainfall described earlier. The area occupied by classes 1 –3 tends to receive the seasonal rainfall last but the area of class 3 misses out on late summer rain from decaying cyclonic depressions, as well as having a lower frequency of remnant and non-

agricultural vegetation which tend to elevate the baseline minimum NDVI. However, the rates of increase and decrease in NDVI become greater from class 1 to class three, as MaxV and SMMaxV are higher, EndV and OnV are lower, the range of increase is greater and HRanTO is later. In summary, class 3 represents cropland that is last to receive rain, but responds rapidly from a very low base NDVI, and where most of the NDVI signal is attributable to crop growth. Class 1 represents marginal land that sometime receives late summer rain but has more mixed land use

M.J. Hill, G.E. Donald / Remote Sensing of Environment 84 (2003) 367–384

meaning that the NDVI signal is a mixture of crop and other land use. The behaviour of the metrics for classes 6– 9, the core wheat and wool growing areas, is more straightforward. Class 6, the northern wheatbelt has a late HRanTO, an early MaxT and a high RIN reflecting the influence the later break of season in the north and warmer temperatures in late autumn and winter. Classes 7 and 9 in the core wool district, have an early HRanTO, a late MaxT, a moderate RIN, a high RDN (since MaxV is very high) and high EndV. They are only distinguished by higher variability for some metrics in class 7 than in class 9. Class 8, the southern wheatbelt, has a later HRanTO, a lower MaxT, a moderate RIN and high RDN. Classes 6 and 8 have lower EndV values than classes 7 and 9 probably reflecting the higher exposure of bare soil in the cropping areas compared with the pastured land of the core wool districts. The examination of these metrics reveals certain codependencies; if maximum values are high and minimum values are low then rates and half range values are high, unless starting dates are very early. The slope metric integrates these properties. Curves can be broad and tall, broad and short, narrow and tall or narrow and short. These four basic shapes encapsulate the major response classes— long growing season with mixed land cover; long growing season with mostly responsive agricultural land cover; shorter growing season with mixed land cover; and shorter growing season with mostly responsive agricultural land cover. These shape characteristics provide a basis for creation of a composite index that may be used to compare growing seasons with each other and the long-term average both in terms of mean response and spatial consistency of response. 3.3. Metrical properties of representative regions The characteristics of the SLAs selected for analysis are shown in Fig. 4. Most of these have a mixture of cropping and livestock production based on land use for 1994, however, Bussleton has virtually no cropping and Merredin has only a small proportion of pasture. The sheep numbers tend to decline to varying degrees for all SLAs due to the downturn in the wool industry in the early 1990s. This is accompanied by a gradual rise in wheat cropping in most areas. The metrics may be used to describe the seasonal production potential for these SLAs by looking at the frequency of occurrence of TINDVI by SLA and year— this provides a picture of the variation in production potential by year (Fig. 9). It can also give an indication of changing land use if the shape of the TINDVI frequency histogram changes permanently. These frequency histograms may be viewed using a simple shape classification: tall and narrow with high, medium or low mean indicating most land is responsive agricultural cover with uniform seasonal conditions; broad and short indicating that land is of mixed cover type or seasonal conditions are not

379

spatially uniform. The shapes can be quantified by simple analysis of the frequency histogram. If the individual years are standardized by subtracting the long-term mean, then mean of these distributions provides a relative rating for the season, and the skewness, relative to the long-term skewness, provides a relative rating of spatial uniformity of response, i.e., if the individual year has a markedly different skewness value to the long-term average, this indicates a large difference in the spatial distribution of response. Fig. 9 shows the frequency histogram for eight SLAs; the long-term mean and skewness value are shown on the first graph, and the deviation from the mean and the skewness value are shown for each individual year. This approach can be adopted to provide a relative rating for any unit of aggregation—in this case, the units of aggregation, SLAs, are different sizes, but the units of aggregation could take any form or shape which provides enough pixels to create an analyzable distribution. The SLAs in Fig. 9 show characteristic frequency histogram patterns and deviations from these patterns for particular years. Moora, Merredin, Brookton and Kojonup have generally tall and narrow histograms indicating predominantly responsive agricultural cover. These histograms move up and down the range of TINDVI according to the quality of the growing season. In 1993 at Kojonup and Busselton, there is distinct spatial variation in quality of season, with a skewness values of 2.22 and 2.0 compared to the long-term value of 0.58 and 0.48, respectively, indicating that vegetation response was both uneven across the SLA and distinctly different to the average. Examples where mean deviation was highly positive and spatial distribution was positively skewed compared to average are given by all SLAs in 1993 except Kojonup. Moora, Merredin, Brookton and Jerramungup in 1996 provide examples of poor responses, but only at Brookton is the skewness of the distribution very different to the average. The histograms for Northampton, Bussleton, Jerramungup and Esperance are generally broader than for the other locations. These histograms show a wider range of deviation in shape. The patterns for Bussleton and Northampton show a consistent lower end tail suggesting that at least parts of these areas have mixed land use contributing to a consistently low TINDVI. The metrics can be used to characterize the seasonal response of each SLA (Fig. 10a,b). The temporal metrics show that the four southern SLAs, Kojonup, Bussleton, Jeramungup and Esperance have longer growing seasons than the four northern SLAs, Brookton, Merredin, Moora and Northampton but more variation in starting date (Fig. 10a). The minimum NDVI is highest for Bussleton with dairying and substantial tree cover, next highest for Kojonup with predominantly pasture cover, and lowest for the Merredin, Moora and Northampton where cropping is widespread (Fig. 10b). These metric assessments of SLAs provide a general picture of the length and timing of seasons and the magnitude of vegetation response that can be expected.

380

M.J. Hill, G.E. Donald / Remote Sensing of Environment 84 (2003) 367–384

Fig. 9. Histograms showing the frequency of TINDVI values within each of eight SLAs for 8 years from 1992 to 1999 and showing characteristic distribution shapes and perturbation of these shapes across the years. Number at the far left indicates the long-term mean TINDVI and the skewness of the long-term mean frequency histogram. Number for each year indicate the difference in mean for that year and the long-term mean, and the skewness of the distribution for that year.

M.J. Hill, G.E. Donald / Remote Sensing of Environment 84 (2003) 367–384

381

Fig. 10. (a) Mean (point), temporal (left-hand bar) standard deviation and spatial (right-hand bar) standard deviation of time-based metrics for eight SLAs in SW Western Australia. (b) Mean (point), temporal (left-hand bar) standard deviation and spatial (right-hand bar) standard deviation of value-based metrics for eight SLAs in SW Western Australia. RIN and RDN are shown as bars and without standard deviations as these were all very similar.

3.4. TINDVI as an indicator of agricultural productivity Total agricultural production in Western Australia on an SLA basis is very well represented by area of SLA (Fig. 11). This is partly a function of the very strong area relationship with livestock numbers, and partly due to the lack of large differences in reliability of cropping, meaning that yields/ unit area are relatively similar, and therefore bigger areas produce bigger total yields. When cell numbers as a measure of area are multiplied by TINDVI as a measure of response, there is some small improvement in the relationship with total production in 1996 and 1997, but not in other years. In the dry year of 1996, TINDVI as a measure of response

appears to add significant information, however in the other years, production is basically a function of area. The key factor here is that in none of the years used in this study, were conditions poor enough to cause crop failure, nonsowing of crops or significant animal losses or transfer. Detection of differences in total production with TINDVI in this usually reliable seasonal climate depends on conditions being severe enough to leave significant areas unsown, or result in significant crop failures. Productivity per hectare provides more interesting and sensitive results (Fig. 11). The effectiveness of TINDVI as a measure of agricultural productivity was dependent on the nature of the growing season and the range of TINDVI

382

M.J. Hill, G.E. Donald / Remote Sensing of Environment 84 (2003) 367–384

Fig. 11. Graphs showing the relationship between area of SLA as represented by cell count and total agricultural production in dse; relationship between TINDVI.cell count and total agricultural production in dse; relationship between TINDVI and agricultural productivity in dse/ha; and relationship between TINDVI and crop productivity in dse/ha (based on wheat, oats, barley and lupins only).

across the agricultural zone. Total agricultural productivity is best represented in the dry year of 1996. The fitted curves are shown only for data points below a threshold rainfall value. This threshold represents the rainfall level above which there is no relationship with productivity in that year. The approximate thresholds were 600 mm for 1992, 1993 and 1994, 700 mm for 1995 and 500 mm for 1996. For 1992, 1994, and

1995, reasonable relationships could be formed between TINDVI and agricultural productivity (Fig. 11). In 1992, a dual response was evident with a negative relationship between productivity and TINDVI above 600 mm and a positive relationship between productivity and TINDVI below 600-mm rainfall. We have previously established that livestock distribution reflects long-term rainfall distribution.

M.J. Hill, G.E. Donald / Remote Sensing of Environment 84 (2003) 367–384

Since our agricultural productivity index is heavily weighted by livestock numbers, this has a large controlling influence on the pattern of values (Hill et al., 1998). When the pattern of TINDVI matches this pattern, agricultural productivity is apparently well represented; but in fact, general productivity associated with long-term average agro-climatic patterns is actually being reflected in both data sets. Analysis of crop productivity alone was restricted to those SLAs with greater than 30% of the area under cropping. In this case, productivity is more generally related to TINDVI except in 1996; however, the best relationship is still obtained for the driest growing season of 1995. A dual response is again evident in 1992 with a negative relationship for eastern areas and a positive relationship for western areas. The results suggest that if rainfall is low enough to influence crop production and the distribution and amount is reflected in TINDVI, then TINDVI can be a good measure of crop productivity. However, if rainfall is more plentiful, and the range of variation in TINDVI is reduced, then sensitivity to crop productivity is lower. In very good seasons, the range of TINDVI becomes compressed towards the top end of values and sensitivity is lost (e.g. Fig. 11— 1993). The main point arising from these analyses is that the relationship between productivity and TINDVI in a highly reliable and seasonal agricultural region is greatly influenced by a number of land cover and land use factors. These include the livestock –cropping mix, the long-term average rainfall, the change with rainfall in the annual/perennial vegetation mix, and the threshold of rainfall below which production is highly sensitive, and above which production is highly insensitive, to rainfall or ‘‘goodness’’ of season.

383

and climatic characteristics. The TINDVI for each zone is a composite of the other metrics defining the timing and magnitude of NDVI response, and can be summarized in four characteristic response forms. 3. The metrics can be used to quantitatively describe the growing season for individual Statistical Local Areas (or other forms of regionalization) indicating both the ‘‘normal’’ and possible range in values for starting dates, ending date, maximum value, etc. These metrics may be used as relative indices for comparison between regions and areas (Reed et al., 2001), however, they need a relationship with agricultural productivity or measures of landscape function. As is, they only represent a general index of non-specific response, but are very useful for comparison (Reed et al., 2001). 4. The properties of the frequency distribution of TINDVI for any region or area of aggregation may be used as a basis for assessment of the relative ranking of the individual years compared to a long-term average in terms of: (a) magnitude of response—the difference between year mean and long-term mean; and (b) distribution of response—the skewness of the distribution relative to the long-term skewness. 5. In the southwest of Western Australia, the TINDVI shows a relationship to agricultural productivity which is dependent on the extent to which crop or total agricultural production are directly reduced by rainfall deficiency. As this never occurred during the 7-year study period, TINDVI proved most sensitive to crop productivity for SLAs having rainfall less than a threshold level of rainfall, usually 600 mm but sometimes a higher or lower value, and in years when rainfall and crop production were highly correlated.

4. Discussion We have assessed the application of metric calculations on NDVI time series to the characterization of seasonal properties and agricultural productivity using the south west of Western Australia as a study area. The NDVI metrics are calculated using established methods (Reed et al., 1994) that extract magnitude and time-based values from the NDVI curves that may be used to summarize seasonal characteristics. The NDVI curves must have a significant amplitude and distinct response period for metric determination—in this study, metrics are not calculated where forests occur or where tree signal is sufficient to dampen the ‘‘greening’’ response to the point where modality is indeterminate. The results of the analysis have shown that: 1. Calculation of metrics for a series of seasons reveals large differences in spatial patterns of TINDVI between seasons, and characterizes regions based on the starting date, timing of maximum value and duration and the variability of these measures. 2. Classification of TINDVI for southwestern WA identifies a set of zones which can be readily related to land use

NDVI metrics for the USA are available via web-based software that enables analysis of interrelationships between seasonal metrics and a range of geographic data sources such as land cover type or administrative unit (Reed et al., 2001). In this context, metrics from standardized NDVI time series could be routinely and transparently used for retrospective assessment of seasonal conditions and changes in vegetation responses and cover. Full information extraction from NDVI time series may require application of a suite of approaches to curve definition and decomposition including classical time series (Roderick et al., 1999), metrics (Reed et al., 1994), time-based classification (Moody & Johnson, 2001) and Fourier and wavelet approaches (Li & Kafatos, 2000). These approaches each have strengths and limitations, and it would be valuable to compare the information yield from application of a range of such methods to the same data sets representing highly seasonal and aseasonal environments. There is a need for further work to compare these approaches on the same data sets, and formulate a framework for application and analytical outputs.

384

M.J. Hill, G.E. Donald / Remote Sensing of Environment 84 (2003) 367–384

Acknowledgements We are grateful to John Sims from the Bureau of Rural Sciences for helpful suggestions on approaches to analysis and presentation of the NDVI metrics for WA. We are grateful to Richard Smith from the WA Department of Land Administration for making NOAA AVHRR NDVI data available for this study.

References Cridland, S. W., Burnside, D. G., Smith, R. C. G. (1994). Use by managers in rangeland environments of near real-time satellite measurements of seasonal vegetation response. Proceedings of the 7th Australasian remote sensing conference, Melbourne, Australia ( pp. 1134 – 1141). Diallo, O., Diouf, A., Hanan, N. P., Ndiaye, A., Pre’vost, Y. (1991). AVHRR monitoring of savanna primary production in Senegal, West Africa: 1987 – 1988. International Journal of Remote Sensing, 12, 1259 – 1279. Hill, M. J., Donald, G. E., Smith, R. C. G. (1998). NDVI-based phenological indices; predictive potential for grazing systems. Proceedings of the 9th Australasian remote sensing and photogrammetry conference, UNSW, Sydney, July, 1998. Hill, M. J., Vickery, P. J., Furnival, E. P., Donald, G. E. (1999). Using of NOAA AVHRR NDVI and classified Landsat TM data to describe pastures in the temperate high rainfall zone (HRZ) of Eastern Australia. Remote Sensing of Environment, 67, 32 – 50. Holben, B. N. (1986). Characteristics of maximum-value composite images from temporal AVHRR data. International Journal of Remote Sensing, 11, 1417 – 1434. Huete, A. R., Liu, H. Q., Batchily, K., van Leeuwen, W. (1997). A comparison of vegetation indices over a global set of TM images for EOSMODIS. Remote Sensing of Environment, 59, 440 – 451. Li, Z., Kafatos, M. (2000). Interannual variability of vegetation in the United States and its relation to El Nino/Southern Oscillation. Remote Sensing of Environment, 71, 239 – 247. Lloyd, D. (1990). A phenological classification of terrestrial vegetation cover using shortwave vegetation index imagery. International Journal of Remote Sensing, 11, 2269 – 2279. Malingreau, J.-P. (1986). Global vegetation dynamics: satellite observations over Asia. International Journal of Remote Sensing, 9, 1121 – 1146. Moody, A., Johnson, D. M. (2001). Land-surface phenologies from AVHRR using the discrete fourier transform. Remote Sensing of Environment, 75, 305 – 323. Prince, S. D. (1991). Satellite remote sensing of primary production: com-

parison of results for Sahelian grasslands 1981 – 1988. International Journal of Remote Sensing, 12, 1301 – 1311. Prince, S. D., Tucker, C. J. (1986). Satellite remote sensing of rangelands in Botswana: II. NOAA AVHRR and herbaceous vegetation. International Journal of Remote Sensing, 7, 1555 – 1570. Reed, B. C., Brown, J. F., VanderZee, D., Loveland, T. R., Merchant, J. W., Ohlen, D. O. (1994). Measuring phenological variability from satellite imagery. Journal of Vegetation Science, 5, 703 – 714. Reed, B. C., Loveland, T. R., Tieszen, L. L. (1996). An approach for using AVHRR data to monitor U.S. Great Plains grasslands. Geocarto International, 11, 1 – 10. Reed, B. C., Swets, D., Bard, L., Brown, J., Rowland, J. (2001). Interactive visualisation of vegetation dynamics. In T. Stein (Ed.), 2001 IEEE international geoscience and remote sensing proceedings. Piscataway, NJ: IEEE. CDROM. Roderick, M. L., Noble, I. R., Cridland, S. W. (1999). Estimating woody and herbaceous vegetation cover from time series satellite observations. Global Ecology and Biogeography, 8, 501 – 508. Roderick, M. L., Smith, R. C. G., Lodwick, G. D. (1996). Calibrating long term AVHRR derived NDVI imagery. Remote Sensing of Environment, 58, 1 – 12. Roderick, M. L., Smith, R., Lodwick, G. (1997a). Annual growth dynamics over continental scale regions using satellite derived vegetation indices. In N. E. West (Ed.), Rangelands in a sustainable biosphere. Proceedings of the fifth international rangelands congress ( pp. 477 – 478). Denver: Society of Range Management. Roderick, M. L., Smith, R., Lodwick, G. (1997b). Rainfall/plant growth relationships using transfer functions and satellite observations. In N. E. West (Ed.), Rangelands in a sustainable biosphere. Proceedings of the fifth international rangelands congress ( pp. 479 – 480). Denver: Society of Range Management. Salmonson, V. V., Guenther, B., Masuoka, E. (2001). A summary of the status of the EOS terra mission Moderate Resolution Imaging Spectroradiometer (MODIS) and attendant data product development after one year of on-orbit performance. IEEE 2001 international geoscience and remote sensing symposium, 9 – 13 July 2001. CDROM. Smith, R. C. G. (1994). Australian Vegetation Watch, Vegetation Watch, Final Report, RIRDC reference No. DOL-1A. Smith, R. C. G., Adams, J., Hick, P. T. (1994). Forecasting wheat yields from satellite observations. Proceedings of the 7th Australasian remote sensing conference, Melbourne, Australia ( pp. 641 – 648). Tieszen, L. L., Reed, B. C., Bliss, N. B., Wylie, B. K., DeJong, D. D. (1997). NDVI, C3 and C4 production, and distributions in Great Plains grassland land cover classes. Ecological Applications, 7, 59 – 78. Yang, L., Wylie, B. K., Tieszen, L. L., Reed, B. C. (1998). An analysis of relationships among climate forcing and time integrated NDVI of grasslands over the US Northern and Central Great Plains. Remote Sensing of Environment, 65, 25 – 37.