- Email: [email protected]
Crop identification using harmonic analysis of time-series AVHRR NDVI data
Computers and Electronics in Agriculture 37 (2002) 127 /139 www.elsevier.com/locate/compag
Crop identification using harmonic analysis of time-series AVHRR NDVI data Mark E. Jakubauskas a, , David R. Legates b, Jude H. Kastens a a
Kansas Applied Remote Sensing (KARS) Program, 2335 Irving Hill Road, University of Kansas, Lawrence, KS 66045, USA b Center for Climatic Research, University of Delaware, Newark, DE, USA
Abstract Harmonic analysis of a time series of National Oceanic and Atmospheric Administration (NOAA) advanced very high resolution radiometer normalized difference vegetation index (NDVI) data was used to develop an innovative technique for crop type identification based on temporal changes in NDVI values. Different crops (corn, soybeans, alfalfa) exhibit distinctive seasonal patterns of NDVI variation that have strong periodic characteristics. Harmonic analysis, or Fourier analysis, decomposes a time-dependent periodic phenomenon into a series of constituent sinusoidal functions, or terms, each defined by a unique amplitude and phase value. Amplitude and phase angle images were produced by analysis of the timeseries NDVI data and used within a discriminant analysis to develop a methodology for crop type identification. For crops that have a single distinct growing season and period of peak greenness, such as corn, the majority of the variance was captured by the first and additive terms, while winter wheat exhibited a bimodal NDVI periodicity with the majority of the variance accounted for by the second harmonic term. # 2002 Elsevier Science B.V. All rights reserved. Keywords: Harmonic analysis; Normalized difference vegetation index; Advanced Very High Resolution Radiometer
Corresponding author. Tel.: /1-785-864-7316; fax: /1-785-864-0392 E-mail address: [email protected] (M.E. Jakubauskas). 0168-1699/02/$ - see front matter # 2002 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 8 - 1 6 9 9 ( 0 2 ) 0 0 1 1 6 - 3
128
M.E. Jakubauskas et al. / Computers and Electronics in Agriculture 37 (2002) 127 /139
1. Introduction Identification and mapping of crop types using coarse-resolution satellite imagery is critical to many applications, from simple estimation of cultivated areas (Loveland et al. 1995, 1991) to stratification for crop yield models (Lee, 1999; Kastens et al., 1998; Doraiswamy and Cook, 1995) to a critical component of mesoscale storm prediction (Gutman and Ignatov, 1998; Gutman et al., 1989) and hydrologic models (Yin and Williams, 1997). Characterization and mapping of crop types using the National Oceanic and Atmospheric Administration (NOAA) Advanced Very HighResolution Radiometer (AVHRR) 1.1 km resolution data has typically taken one of two approaches: temporal profiles of crop phenology as manifested in the normalized difference vegetation index (NDVI) (DeFries et al., 1995; Reed et al., 1994; Lloyd, 1990); and classification of multitemporal data using maximum likelihood or other classification algorithms (Loveland et al., 1995; Brown et al., 1993; Loveland et al., 1991). To date, few researchers in remote sensing have applied time-series analysis techniques (harmonic or Fourier analysis) developed for other disciplines, such as electrical engineering, hydrology, or climatology (Legates and Willmott, 1990; Davis, 1986), to the NDVI time series. This paper describes the application of harmonic analysis to a single year (1992) of NOAA-AVHRR data (26 periods) for identification of several common crop types occurring within the western Great Plains. The products of the harmonic analysis (additive, amplitude, and phase terms) are used within a discriminant analysis for mapping crop types within Finney County in southwestern Kansas. Briefly defined, harmonic (Fourier) analysis permits a complex curve to be expressed as the sum of a series of cosine waves (terms) and an additive term (Davis, 1986; Rayner, 1971). Each wave is defined by a unique amplitude and a phase angle, where the amplitude value is half the height of a wave, and the phase angle (or simply, phase) defines the offset between the origin and the peak of the wave over the range 0/2p (Fig. 1a). Each term designates the number of complete cycles completed by a wave over the defined interval (e.g., the second term completes two cycles, Fig. 1b). Successive harmonic terms are added to produce a complex curve (Fig. 1c), and each component curve, or term, accounts for a percentage of the total variance in the original time-series data set.
2. Methods 2.1. Study area Finney County, the second-largest county in Kansas at approximately 1300 mile2, is located in southwestern Kansas in the High Plains at approximately 1008W longitude. Cropland dominates the relatively flat landscape of Finney County, comprising over 76% of the county. The Arkansas River extends east /west through the central portion of the county (Fig. 2), and irrigated agriculture (predominantly corn, milo, and alfalfa) is clustered along the river, drawing on the Ogalalla Aquifer
M.E. Jakubauskas et al. / Computers and Electronics in Agriculture 37 (2002) 127 /139
129
Fig. 1. (a) Simple cosine curve representative of the first harmonic; (b) curves for harmonic terms 1, 2, and 3; (c) curve produced by addition of curves in (b).
to supply center-pivot irrigation systems. Dryland agriculture dominates the northern and southern parts of the county, with winter wheat as the principal crop. Some native shortgrass prairie remains in the northeast where the topography is too rugged for agriculture, and the land is used for grazing. 2.2. AVHRR times-series NDVI data The NDVI is a well-established and commonly used vegetation index in studies using remote sensing data because it is roughly correlated with green plant biomass and vegetation cover (Box et al., 1989; Tucker, 1979). The NDVI is based on the
130
M.E. Jakubauskas et al. / Computers and Electronics in Agriculture 37 (2002) 127 /139
Fig. 2. Finney County in southwest Kansas.
relative reflectance values in the red and near infrared (NIR) wavelengths: NDVI / (NIR /Red)/(NIR/Red). Vegetation indices are commonly used to reduce effects of atmospheric conditions or different soil backgrounds on spectral reflectance values. The amount of red solar energy reflected by vegetation cover depends primarily on chlorophyll content, whereas the amount of NIR energy reflected by vegetation is affected by the amount and condition of green biomass, leaf tissue structure, and water content (Jensen, 1996). NOAA AVHRR NDVI biweekly composites for 1992 were acquired from the US Geological Survey EROS Data Center in Sioux Falls, South Dakota. The AVHRR system acquires images on a daily basis, but cloud cover often obscures the land surface. In order to produce a cloud-free image of vegetation conditions during a 2week period, a set of 14 daily images is used to create a single biweekly composite image. Each biweekly NDVI composite consists of the maximum NDVI value within a defined 2-week period for each pixel (Eidenshink, 1992). Vegetation index data are rescaled by EROS during processing from a range of /1.0 to /1.0, to 0 to 200. Values less than 100 typically represent snow, ice, water, and other non-vegetated earth surfaces. Data for Finney County were downloaded from CD-ROMs, coregistered, and images for missing biweekly composite periods were created by averaging image data for periods bracketing the missing period (e.g., the previous and succeeding periods).
M.E. Jakubauskas et al. / Computers and Electronics in Agriculture 37 (2002) 127 /139
131
2.3. Harmonic analysis Images of the additive term, and amplitude and phase angle for each term to the seventh harmonic were produced on a per-pixel basis for each pixel in the NDVI data set, using the procedures described by Jakubauskas et al. (2001), modified from Davis (1986). Similar to principal component analysis, the majority of the variance in a data set is contained in the first few terms (components). Over 90% of the variance in the original time series was captured in the additive and the first two harmonic terms. The additive term, and the phase and amplitude for the first two harmonic terms were extracted and used for further analysis. 2.4. Generation of test and training samples Because farm fields in the Great Plains are smaller than the resolution of the AVHRR sensor, most AVHRR pixels are typically mixtures of two or more different land use/cover types. To identify representative pixels where the phenological response for a pixel as recorded by the AVHRR NDVI time series resulted from predominantly one land use/cover type, the 30-meter resolution 1992 land use/land cover map produced by Egbert et al. (1995) was used to calculate the percent cover of each land use/land cover type within each AVHRR 1000-meter pixel for Finney County (Fig. 3). The 1992 land use/land cover map was created by supervised and
Fig. 3. Finney County 1992 land use/cover map from Egbert et al. (1995), resampled to 1000-meter resolution by majority rule.
132
M.E. Jakubauskas et al. / Computers and Electronics in Agriculture 37 (2002) 127 /139
unsupervised classification of three Landsat 5 TM scenes from May, July, and September 1992. Overall accuracy of this map was 98%, and accuracies for individual classes ranged from a low of 84% (corn) to 99% (wheat). For the discriminant analysis, a set of training samples was generated using the percent cover map to identify pixels with greater than 75% cover of a single land use/cover class (e.g. pixels where a single cover type predominated). Stepwise discriminant analysis was used to determine the probability of a pixel belonging to one of the crop types, using the additive term and the phase and amplitude for the first two terms as predictor variables. Discriminant functions were calculated in SPSSTM (SPSS, Chicago, IL) from the set of training samples, and applied to all pixels in the data set, excluding ‘background’ pixels (values of 0). The predicted class membership for each pixel was then retransformed from an SPSS variable back to an ERDASTM (ERDAS, Atlanta, GA) .img image. For comparison to the discriminant analysis map, the Egbert et al. (1995) 30-meter land use/land cover map was degraded to a 1000-meter resolution map. Degradation was performed by assigning to an output 1000-meter pixel the value of the majority of the 30-meter pixels occurring within the co-registered 1000-meter pixel. Classification accuracy was determined by crosstabulation between the discriminant analysis map and the 1000-meter majority class map, calculating overall accuracy and producers and users accuracies for each class.
3. Results and discussion 3.1. Harmonic characteristics of crop types Irrigated corn and alfalfa in Finney County possess a single distinct growing season, attaining peak greenness during midsummer. This is manifested in the harmonic analysis by a strongly unimodal periodic pattern, with a high amplitude value in the first term and low amplitude values in successive terms (Fig. 4a). The majority of the total variance in seasonal NDVI for corn is contained in the first harmonic term. Alfalfa, although it also exhibits a strong summer-peak greenness pattern (Fig. 4b), differs from the strongly unimodal pattern of corn and grass in that it has relatively high second- and third-term amplitudes resulting from cultivation practices. Alfalfa is harvested repeatedly during the summer, producing a secondary and tertiary periodicity in greenness with the growth /cutting/growth /cutting cycle. Irrigated alfalfa is typically harvested four times in western Kansas during a typical season, beginning in May and extending until September. Grassland exhibits a unimodal phenological pattern similar to corn (Fig. 4c), with high first-term amplitude values and the majority of the variance (/90%) in the first harmonic term. First-term phase values indicate that peak greenness for grasslands occurs during July, consistent with the warm-season nature of the grass species. Winter wheat has strikingly different phenological patterns in comparison to other crop types described above in that it exhibits a bimodal temporal NDVI curve. Winter wheat is planted in the fall (late September/early October), sprouts, and goes
M.E. Jakubauskas et al. / Computers and Electronics in Agriculture 37 (2002) 127 /139
133
Fig. 4. First three harmonic curves for: (a) corn; (b) alfalfa; (c) grassland; (d) winter wheat.
dormant over winter and may be covered by snow. In the spring, the wheat greens up and is harvested by May, followed by tilling or fallow. The second harmonic term, representing the spring growth period, therefore has the highest amplitude and contains the majority of the variance (Fig. 4d). 3.2. Harmonic additive, amplitude, and phase images Within Finney County, patterns of the additive term generally follow irrigated/ non-irrigated land use patterns (Fig. 5). The additive, or zero term, is the arithmetic mean of NDVI over the time series (26 periods). High values, in particular irrigated cropland along the Arkansas River, indicate high mean NDVI over a season. Low additive term values are manifested by less vegetated urban areas and the sandsage prairie south of Garden City. The first amplitude term image also closely matches land use/cover patterns within the study area. High first term amplitude values indicate a unimodal temporal NDVI pattern, where a land use/cover has a wide annual range in NDVI values. The highest first term amplitude values (white areas on Fig. 6) occur on areas planted in corn and alfalfa, while the native grasslands in
134
M.E. Jakubauskas et al. / Computers and Electronics in Agriculture 37 (2002) 127 /139
Fig. 5. Additive term image.
Fig. 6. First harmonic amplitude image.
M.E. Jakubauskas et al. / Computers and Electronics in Agriculture 37 (2002) 127 /139
135
the northeastern part of the county appear in light grey. Dark tones on the first term amplitude generally appear white or light grey on the second amplitude image, indicating semiannual peaks in greenness, a phenomena exhibited in this region by winter wheat. Phase indicates the time of year at which the peak value for a term occurs. For land cover types with unimodal NDVI phenological profiles that peak in midsummer (such as irrigated corn or grasslands), phase values typically are near p, or 3.14, over a possible range of 0 /2p. Although grasslands and irrigated cropland in Finney County exhibit different additive term values and first term amplitude values, first term phase values are similar, indicating the midsummer peak greenness timing of these cover types (grey areas, Fig. 7). Winter wheat, which dominates the northwestern and central portions of Finney County, appears dark grey to black on the phase image, indicating a first-term amplitude peak earlier in the year consistent with the spring maturing of the winter wheat crop.
3.3. Crop type mapping General patterns of agricultural land use/cover, as mapped by the discriminant analysis of the harmonic components (Fig. 8), match those of the degraded 1000-m Egbert et al. (1995) map (Fig. 3), although results from the accuracy assessment crosstabulation produced poor classification accuracies (Table 1). Overall class accuracy was 52%, comparable to accuracies produced for cultivated land mapping using AVHRR data in other land use/cover mapping studies (e.g., DeFries et al.,
Fig. 7. First harmonic phase image.
136
M.E. Jakubauskas et al. / Computers and Electronics in Agriculture 37 (2002) 127 /139
Fig. 8. Finney County 1992 crop type map produced by discriminant analysis of harmonic components derived from AVHRR NDVI data.
1995). Irrigated crops (milo, corn, alfalfa) and dryland agriculture (wheat and fallow) typically were misclassified within each broad category of agricultural land use, e.g., alfalfa misclassified with corn, and corn misclassified as milo (Table 1). Wheat and fallow lands exhibited a high degree of cross-classification, with 25% of wheat pixels misclassified as fallow (300/1155), and a comparable percentage (22%) of fallow misclassified as wheat (107/486). Fallow lands (producer’s accuracy /50%) were also misclassified as grasslands (90/486 /18%). Grassland, with the highest individual producer’s accuracy at 72%, misclassified with both dryland and irrigated crops. Results should be considered, however, in light of the difficulty of mapping and performing an accuracy assessment for coarse-resolution data such as the AVHRR. Crop type mapping is often confounded by the spatial resolution of the fields being smaller than the resolution of the scanner, a problem recognized by other researchers (Atkinson et al., 1997; Kremer and Running, 1993; Loveland et al., 1991). Fields that are smaller than the resolution of the sensor produce a composite NDVI response if multiple land use/cover types exist within the resolution of the AVHRR pixel, producing misclassification. In this study, highest classification accuracies were attained for grassland, which in Finney County is concentrated primarily in three large polygons (Fig. 3). Conversely, milo had the lowest classification accuracy (27%, producer’s accuracy) and occurs in small, isolated fields distributed across the county. This suggests that accuracies may be greatly enhanced by applying the
Reference class
Predicted class Wheat Milo Corn Alfalfa Fallow Grass Total Producer’s accuracy (%)
Wheat
Milo
Corn
Alfalfa
Fallow
Grass
Total
User’s accuracy (%)
471 75 97 33 300 179 1155 41
33 33 18 5 10 22 121 27
37 49 120 15 10 43 274 44
12 3 38 57 0 9 119 48
107 31 15 2 241 90 486 50
49 20 86 29 86 687 957 72
709 211 374 141 647 1030 3112
66 16 32 40 37 67 52
M.E. Jakubauskas et al. / Computers and Electronics in Agriculture 37 (2002) 127 /139
Table 1 Accuracy assessment for crop type mapping by discriminant analysis of harmonic components
137
138
M.E. Jakubauskas et al. / Computers and Electronics in Agriculture 37 (2002) 127 /139
methods used here to the newly-available Terra Moderate-Resolution Imaging Spectrometer data, which collects daily reflectance data at a much finer spatial resolution of 250 meters, as opposed to the 1000-m resolution of the AVHRR. Beyond crop type identification and mapping, harmonic analysis of time-series data offers considerable promise as a tool for monitoring vegetation change. Given that different land use/cover types produce different phase and amplitude values in the harmonic analysis, results suggest that harmonic analysis of time-series remotely sensed data could be used to detect changes in land use/land cover by examining changes in annual values of amplitude, phase, or the additive term over a period of years. Changes in seasonal amplitude, with phase unchanged, could indicate changes in crop type or vegetation condition. Changes in phase, with amplitude unchanged, could indicate changes in time of maximum greenness, which in turn may occur as a result of changes in the onset of greenness, changes in planting time, or changes in the time of harvest. Changes in both amplitude and phase values over a period of years to decades would be indicative of major and significant changes in land surface condition. Such radical change could result from changes in land management (crop rotation, enrollment in conservation programs), loss of vegetation by natural or anthropogenic disturbance, or changes in regional climate itself driving changes in vegetation.
Acknowledgements This project was conducted at the Kansas Applied Remote Sensing (KARS) Program (Edward A. Martinko, Director). The research described in this paper was funded by the National Institute for Global Environmental Change (NIGEC) (South Central Regional Center, Tulane University, David Sailor, Director), through the US Department of Energy (Cooperative Agreement No. DE-FC03- 90ER61010). Any opinions, findings, and conclusions or recommendations expressed in this publication are those of the authors and do not necessarily reflect the views of the DOE. Dr Stephen L. Egbert graciously permitted the use of his 1992 high-resolution land use/land cover classification map.
References Atkinson, P.M., Cutler, M.E., Lewis, H., 1997. Mapping sub-pixel proportional land cover with AVHRR imagery. International Journal of Remote Sensing 18 (4), 917 /935. Box, E., Holben, B., Kalb, V., 1989. Accuracy of the AVHRR vegetation index as a predictor of biomass, primary productivity and net CO2 flux. Vegetation 80, 71 /89. Brown, J.F., Loveland, T.R., Merchant, J.W., Reed, B.C., Ohlen, D.O., 1993. Using multi-source data in global land cover characterization: concepts, requirements and methods. Photogrammetric Engineering and Remote Sensing 59, 977 /987. Davis, J.C., 1986. Statistics and data analysis in geology, second ed. Wiley, New York. DeFries, R., Hansen, M., Townshend, J., 1995. Global discrimination of land cover types from metrics derived from AVHRR Pathfinder data. Remote Sensing of Environment 54, 209 /222.
M.E. Jakubauskas et al. / Computers and Electronics in Agriculture 37 (2002) 127 /139
139
Doraiswamy, P.C., Cook, P.W., 1995. Spring wheat yield assessment using NOAA AVHRR data. Canadian Journal of Remote Sensing 21 (1), 43 /51. Egbert, S.E., Price, K.P., Nellis, M.D., Lee, R. 1995. Developing a land cover modeling protocol for the High Plains using multiseasonal Thematic Mapper imagery. Proceedings of the 1995 ASPRS/ACSM Annual Meeting, Charlotte, N.C. pp. 836 /845. Eidenshink, J.C., 1992. The 1990 conterminous US AVHRR data set. Photogrammetric Engineering and Remote Sensing 58 (6), 809 /813. Gutman, G., Ignatov, A., 1998. Derivation of green vegetation fraction from NOAA-AVHRR for use in numerical weather prediction models. International Journal of Remote Sensing 19, 1533. Gutman, G., Ohring, G., Tarpley, D., Ambroziak, R., 1989. Albedo of the US Great Plains as determined from NOAA-9 AVHRR data. Journal of Climate 2, 608 /617. Jakubauskas, M.E., Legates, D.R., Kastens, J., 2001. Harmonic analysis of time-series AVHRR NDVI data. Photogrammetric Engineering and Remote Sensing 67 (4), 461 /470. Jensen, J., 1996. Introductory Digital Image Processing. Prentice Hall, Englewood Cliffs, NJ, p. 316. Kastens, D.L., Price, K.P., Martinko, E.A., 1998. Assessing wheat conditions in Kansas using biweekly AVHRR datasets and crop phenological indices. First International Conference on Geospatial Information in Agriculture and Forestry, Orlando, FL, June 1 /3, 1998. Kremer, R.G., Running, S.W., 1993. Community type differentiation using NOAA/AVHRR data within a sagebrush-steppe ecosystem. Remote Sensing of Environment 46, 311 /318. Lee, R. 1999. Modeling Corn Yields in Iowa Using Time-series Analysis of AVHRR Data and Vegetation Phenological Metrics. PhD Dissertation, Department of Geography, University of Kansas, Lawrence, KS. Legates, D.R., Willmott, C.J., 1990. Mean seasonal and spatial variability in global surface air temperature. Theoretical and Applied Climatology 41 (1), 11 /21. Lloyd, D., 1990. A phenological classification of terrestrial vegetation cover using shortwave vegetation index imagery. International Journal of Remote Sensing 11 (12), 2269 /2279. Loveland, T.R., Merchant, J.W., Ohlen, D.O., Brown, J.F., 1991. Development of a land-cover characteristics database for the conterminous US. Photogrammetric Engineering and Remote Sensing 57, 1453 /1463. Loveland, T.R., Merchant, J.W., Brown, J.F., Ohlen, D.O., Reed, B.C., Olson, P., Hutchinson, J., 1995. Seasonal land-cover regions of the United States. Annals of the Association of American Geographers 85, 339 /355. Rayner, J.N., 1971. An Introduction to Spectral Analysis. Pion Ltd, London, p. 174. 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. Tucker, C.J., 1979. Red and photographic infrared linear combinations for monitoring vegetation. Remote Sensing of Environment 8, 127 /150. Yin, Z., Williams, T.H.L., 1997. Obtaining spatial and temporal vegetation data from Landsat MSS and AVHRR/NOAA satellite images for a hydrologic model. Photogrammetric Engineering and Remote Sensing 63 (1), 69 /77.
Crop identification using harmonic analysis of time-series AVHRR NDVI data Mark E. Jakubauskas a, , David R. Legates b, Jude H. Kastens a a
Kansas Applied Remote Sensing (KARS) Program, 2335 Irving Hill Road, University of Kansas, Lawrence, KS 66045, USA b Center for Climatic Research, University of Delaware, Newark, DE, USA
Abstract Harmonic analysis of a time series of National Oceanic and Atmospheric Administration (NOAA) advanced very high resolution radiometer normalized difference vegetation index (NDVI) data was used to develop an innovative technique for crop type identification based on temporal changes in NDVI values. Different crops (corn, soybeans, alfalfa) exhibit distinctive seasonal patterns of NDVI variation that have strong periodic characteristics. Harmonic analysis, or Fourier analysis, decomposes a time-dependent periodic phenomenon into a series of constituent sinusoidal functions, or terms, each defined by a unique amplitude and phase value. Amplitude and phase angle images were produced by analysis of the timeseries NDVI data and used within a discriminant analysis to develop a methodology for crop type identification. For crops that have a single distinct growing season and period of peak greenness, such as corn, the majority of the variance was captured by the first and additive terms, while winter wheat exhibited a bimodal NDVI periodicity with the majority of the variance accounted for by the second harmonic term. # 2002 Elsevier Science B.V. All rights reserved. Keywords: Harmonic analysis; Normalized difference vegetation index; Advanced Very High Resolution Radiometer
Corresponding author. Tel.: /1-785-864-7316; fax: /1-785-864-0392 E-mail address: [email protected] (M.E. Jakubauskas). 0168-1699/02/$ - see front matter # 2002 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 8 - 1 6 9 9 ( 0 2 ) 0 0 1 1 6 - 3
128
M.E. Jakubauskas et al. / Computers and Electronics in Agriculture 37 (2002) 127 /139
1. Introduction Identification and mapping of crop types using coarse-resolution satellite imagery is critical to many applications, from simple estimation of cultivated areas (Loveland et al. 1995, 1991) to stratification for crop yield models (Lee, 1999; Kastens et al., 1998; Doraiswamy and Cook, 1995) to a critical component of mesoscale storm prediction (Gutman and Ignatov, 1998; Gutman et al., 1989) and hydrologic models (Yin and Williams, 1997). Characterization and mapping of crop types using the National Oceanic and Atmospheric Administration (NOAA) Advanced Very HighResolution Radiometer (AVHRR) 1.1 km resolution data has typically taken one of two approaches: temporal profiles of crop phenology as manifested in the normalized difference vegetation index (NDVI) (DeFries et al., 1995; Reed et al., 1994; Lloyd, 1990); and classification of multitemporal data using maximum likelihood or other classification algorithms (Loveland et al., 1995; Brown et al., 1993; Loveland et al., 1991). To date, few researchers in remote sensing have applied time-series analysis techniques (harmonic or Fourier analysis) developed for other disciplines, such as electrical engineering, hydrology, or climatology (Legates and Willmott, 1990; Davis, 1986), to the NDVI time series. This paper describes the application of harmonic analysis to a single year (1992) of NOAA-AVHRR data (26 periods) for identification of several common crop types occurring within the western Great Plains. The products of the harmonic analysis (additive, amplitude, and phase terms) are used within a discriminant analysis for mapping crop types within Finney County in southwestern Kansas. Briefly defined, harmonic (Fourier) analysis permits a complex curve to be expressed as the sum of a series of cosine waves (terms) and an additive term (Davis, 1986; Rayner, 1971). Each wave is defined by a unique amplitude and a phase angle, where the amplitude value is half the height of a wave, and the phase angle (or simply, phase) defines the offset between the origin and the peak of the wave over the range 0/2p (Fig. 1a). Each term designates the number of complete cycles completed by a wave over the defined interval (e.g., the second term completes two cycles, Fig. 1b). Successive harmonic terms are added to produce a complex curve (Fig. 1c), and each component curve, or term, accounts for a percentage of the total variance in the original time-series data set.
2. Methods 2.1. Study area Finney County, the second-largest county in Kansas at approximately 1300 mile2, is located in southwestern Kansas in the High Plains at approximately 1008W longitude. Cropland dominates the relatively flat landscape of Finney County, comprising over 76% of the county. The Arkansas River extends east /west through the central portion of the county (Fig. 2), and irrigated agriculture (predominantly corn, milo, and alfalfa) is clustered along the river, drawing on the Ogalalla Aquifer
M.E. Jakubauskas et al. / Computers and Electronics in Agriculture 37 (2002) 127 /139
129
Fig. 1. (a) Simple cosine curve representative of the first harmonic; (b) curves for harmonic terms 1, 2, and 3; (c) curve produced by addition of curves in (b).
to supply center-pivot irrigation systems. Dryland agriculture dominates the northern and southern parts of the county, with winter wheat as the principal crop. Some native shortgrass prairie remains in the northeast where the topography is too rugged for agriculture, and the land is used for grazing. 2.2. AVHRR times-series NDVI data The NDVI is a well-established and commonly used vegetation index in studies using remote sensing data because it is roughly correlated with green plant biomass and vegetation cover (Box et al., 1989; Tucker, 1979). The NDVI is based on the
130
M.E. Jakubauskas et al. / Computers and Electronics in Agriculture 37 (2002) 127 /139
Fig. 2. Finney County in southwest Kansas.
relative reflectance values in the red and near infrared (NIR) wavelengths: NDVI / (NIR /Red)/(NIR/Red). Vegetation indices are commonly used to reduce effects of atmospheric conditions or different soil backgrounds on spectral reflectance values. The amount of red solar energy reflected by vegetation cover depends primarily on chlorophyll content, whereas the amount of NIR energy reflected by vegetation is affected by the amount and condition of green biomass, leaf tissue structure, and water content (Jensen, 1996). NOAA AVHRR NDVI biweekly composites for 1992 were acquired from the US Geological Survey EROS Data Center in Sioux Falls, South Dakota. The AVHRR system acquires images on a daily basis, but cloud cover often obscures the land surface. In order to produce a cloud-free image of vegetation conditions during a 2week period, a set of 14 daily images is used to create a single biweekly composite image. Each biweekly NDVI composite consists of the maximum NDVI value within a defined 2-week period for each pixel (Eidenshink, 1992). Vegetation index data are rescaled by EROS during processing from a range of /1.0 to /1.0, to 0 to 200. Values less than 100 typically represent snow, ice, water, and other non-vegetated earth surfaces. Data for Finney County were downloaded from CD-ROMs, coregistered, and images for missing biweekly composite periods were created by averaging image data for periods bracketing the missing period (e.g., the previous and succeeding periods).
M.E. Jakubauskas et al. / Computers and Electronics in Agriculture 37 (2002) 127 /139
131
2.3. Harmonic analysis Images of the additive term, and amplitude and phase angle for each term to the seventh harmonic were produced on a per-pixel basis for each pixel in the NDVI data set, using the procedures described by Jakubauskas et al. (2001), modified from Davis (1986). Similar to principal component analysis, the majority of the variance in a data set is contained in the first few terms (components). Over 90% of the variance in the original time series was captured in the additive and the first two harmonic terms. The additive term, and the phase and amplitude for the first two harmonic terms were extracted and used for further analysis. 2.4. Generation of test and training samples Because farm fields in the Great Plains are smaller than the resolution of the AVHRR sensor, most AVHRR pixels are typically mixtures of two or more different land use/cover types. To identify representative pixels where the phenological response for a pixel as recorded by the AVHRR NDVI time series resulted from predominantly one land use/cover type, the 30-meter resolution 1992 land use/land cover map produced by Egbert et al. (1995) was used to calculate the percent cover of each land use/land cover type within each AVHRR 1000-meter pixel for Finney County (Fig. 3). The 1992 land use/land cover map was created by supervised and
Fig. 3. Finney County 1992 land use/cover map from Egbert et al. (1995), resampled to 1000-meter resolution by majority rule.
132
M.E. Jakubauskas et al. / Computers and Electronics in Agriculture 37 (2002) 127 /139
unsupervised classification of three Landsat 5 TM scenes from May, July, and September 1992. Overall accuracy of this map was 98%, and accuracies for individual classes ranged from a low of 84% (corn) to 99% (wheat). For the discriminant analysis, a set of training samples was generated using the percent cover map to identify pixels with greater than 75% cover of a single land use/cover class (e.g. pixels where a single cover type predominated). Stepwise discriminant analysis was used to determine the probability of a pixel belonging to one of the crop types, using the additive term and the phase and amplitude for the first two terms as predictor variables. Discriminant functions were calculated in SPSSTM (SPSS, Chicago, IL) from the set of training samples, and applied to all pixels in the data set, excluding ‘background’ pixels (values of 0). The predicted class membership for each pixel was then retransformed from an SPSS variable back to an ERDASTM (ERDAS, Atlanta, GA) .img image. For comparison to the discriminant analysis map, the Egbert et al. (1995) 30-meter land use/land cover map was degraded to a 1000-meter resolution map. Degradation was performed by assigning to an output 1000-meter pixel the value of the majority of the 30-meter pixels occurring within the co-registered 1000-meter pixel. Classification accuracy was determined by crosstabulation between the discriminant analysis map and the 1000-meter majority class map, calculating overall accuracy and producers and users accuracies for each class.
3. Results and discussion 3.1. Harmonic characteristics of crop types Irrigated corn and alfalfa in Finney County possess a single distinct growing season, attaining peak greenness during midsummer. This is manifested in the harmonic analysis by a strongly unimodal periodic pattern, with a high amplitude value in the first term and low amplitude values in successive terms (Fig. 4a). The majority of the total variance in seasonal NDVI for corn is contained in the first harmonic term. Alfalfa, although it also exhibits a strong summer-peak greenness pattern (Fig. 4b), differs from the strongly unimodal pattern of corn and grass in that it has relatively high second- and third-term amplitudes resulting from cultivation practices. Alfalfa is harvested repeatedly during the summer, producing a secondary and tertiary periodicity in greenness with the growth /cutting/growth /cutting cycle. Irrigated alfalfa is typically harvested four times in western Kansas during a typical season, beginning in May and extending until September. Grassland exhibits a unimodal phenological pattern similar to corn (Fig. 4c), with high first-term amplitude values and the majority of the variance (/90%) in the first harmonic term. First-term phase values indicate that peak greenness for grasslands occurs during July, consistent with the warm-season nature of the grass species. Winter wheat has strikingly different phenological patterns in comparison to other crop types described above in that it exhibits a bimodal temporal NDVI curve. Winter wheat is planted in the fall (late September/early October), sprouts, and goes
M.E. Jakubauskas et al. / Computers and Electronics in Agriculture 37 (2002) 127 /139
133
Fig. 4. First three harmonic curves for: (a) corn; (b) alfalfa; (c) grassland; (d) winter wheat.
dormant over winter and may be covered by snow. In the spring, the wheat greens up and is harvested by May, followed by tilling or fallow. The second harmonic term, representing the spring growth period, therefore has the highest amplitude and contains the majority of the variance (Fig. 4d). 3.2. Harmonic additive, amplitude, and phase images Within Finney County, patterns of the additive term generally follow irrigated/ non-irrigated land use patterns (Fig. 5). The additive, or zero term, is the arithmetic mean of NDVI over the time series (26 periods). High values, in particular irrigated cropland along the Arkansas River, indicate high mean NDVI over a season. Low additive term values are manifested by less vegetated urban areas and the sandsage prairie south of Garden City. The first amplitude term image also closely matches land use/cover patterns within the study area. High first term amplitude values indicate a unimodal temporal NDVI pattern, where a land use/cover has a wide annual range in NDVI values. The highest first term amplitude values (white areas on Fig. 6) occur on areas planted in corn and alfalfa, while the native grasslands in
134
M.E. Jakubauskas et al. / Computers and Electronics in Agriculture 37 (2002) 127 /139
Fig. 5. Additive term image.
Fig. 6. First harmonic amplitude image.
M.E. Jakubauskas et al. / Computers and Electronics in Agriculture 37 (2002) 127 /139
135
the northeastern part of the county appear in light grey. Dark tones on the first term amplitude generally appear white or light grey on the second amplitude image, indicating semiannual peaks in greenness, a phenomena exhibited in this region by winter wheat. Phase indicates the time of year at which the peak value for a term occurs. For land cover types with unimodal NDVI phenological profiles that peak in midsummer (such as irrigated corn or grasslands), phase values typically are near p, or 3.14, over a possible range of 0 /2p. Although grasslands and irrigated cropland in Finney County exhibit different additive term values and first term amplitude values, first term phase values are similar, indicating the midsummer peak greenness timing of these cover types (grey areas, Fig. 7). Winter wheat, which dominates the northwestern and central portions of Finney County, appears dark grey to black on the phase image, indicating a first-term amplitude peak earlier in the year consistent with the spring maturing of the winter wheat crop.
3.3. Crop type mapping General patterns of agricultural land use/cover, as mapped by the discriminant analysis of the harmonic components (Fig. 8), match those of the degraded 1000-m Egbert et al. (1995) map (Fig. 3), although results from the accuracy assessment crosstabulation produced poor classification accuracies (Table 1). Overall class accuracy was 52%, comparable to accuracies produced for cultivated land mapping using AVHRR data in other land use/cover mapping studies (e.g., DeFries et al.,
Fig. 7. First harmonic phase image.
136
M.E. Jakubauskas et al. / Computers and Electronics in Agriculture 37 (2002) 127 /139
Fig. 8. Finney County 1992 crop type map produced by discriminant analysis of harmonic components derived from AVHRR NDVI data.
1995). Irrigated crops (milo, corn, alfalfa) and dryland agriculture (wheat and fallow) typically were misclassified within each broad category of agricultural land use, e.g., alfalfa misclassified with corn, and corn misclassified as milo (Table 1). Wheat and fallow lands exhibited a high degree of cross-classification, with 25% of wheat pixels misclassified as fallow (300/1155), and a comparable percentage (22%) of fallow misclassified as wheat (107/486). Fallow lands (producer’s accuracy /50%) were also misclassified as grasslands (90/486 /18%). Grassland, with the highest individual producer’s accuracy at 72%, misclassified with both dryland and irrigated crops. Results should be considered, however, in light of the difficulty of mapping and performing an accuracy assessment for coarse-resolution data such as the AVHRR. Crop type mapping is often confounded by the spatial resolution of the fields being smaller than the resolution of the scanner, a problem recognized by other researchers (Atkinson et al., 1997; Kremer and Running, 1993; Loveland et al., 1991). Fields that are smaller than the resolution of the sensor produce a composite NDVI response if multiple land use/cover types exist within the resolution of the AVHRR pixel, producing misclassification. In this study, highest classification accuracies were attained for grassland, which in Finney County is concentrated primarily in three large polygons (Fig. 3). Conversely, milo had the lowest classification accuracy (27%, producer’s accuracy) and occurs in small, isolated fields distributed across the county. This suggests that accuracies may be greatly enhanced by applying the
Reference class
Predicted class Wheat Milo Corn Alfalfa Fallow Grass Total Producer’s accuracy (%)
Wheat
Milo
Corn
Alfalfa
Fallow
Grass
Total
User’s accuracy (%)
471 75 97 33 300 179 1155 41
33 33 18 5 10 22 121 27
37 49 120 15 10 43 274 44
12 3 38 57 0 9 119 48
107 31 15 2 241 90 486 50
49 20 86 29 86 687 957 72
709 211 374 141 647 1030 3112
66 16 32 40 37 67 52
M.E. Jakubauskas et al. / Computers and Electronics in Agriculture 37 (2002) 127 /139
Table 1 Accuracy assessment for crop type mapping by discriminant analysis of harmonic components
137
138
M.E. Jakubauskas et al. / Computers and Electronics in Agriculture 37 (2002) 127 /139
methods used here to the newly-available Terra Moderate-Resolution Imaging Spectrometer data, which collects daily reflectance data at a much finer spatial resolution of 250 meters, as opposed to the 1000-m resolution of the AVHRR. Beyond crop type identification and mapping, harmonic analysis of time-series data offers considerable promise as a tool for monitoring vegetation change. Given that different land use/cover types produce different phase and amplitude values in the harmonic analysis, results suggest that harmonic analysis of time-series remotely sensed data could be used to detect changes in land use/land cover by examining changes in annual values of amplitude, phase, or the additive term over a period of years. Changes in seasonal amplitude, with phase unchanged, could indicate changes in crop type or vegetation condition. Changes in phase, with amplitude unchanged, could indicate changes in time of maximum greenness, which in turn may occur as a result of changes in the onset of greenness, changes in planting time, or changes in the time of harvest. Changes in both amplitude and phase values over a period of years to decades would be indicative of major and significant changes in land surface condition. Such radical change could result from changes in land management (crop rotation, enrollment in conservation programs), loss of vegetation by natural or anthropogenic disturbance, or changes in regional climate itself driving changes in vegetation.
Acknowledgements This project was conducted at the Kansas Applied Remote Sensing (KARS) Program (Edward A. Martinko, Director). The research described in this paper was funded by the National Institute for Global Environmental Change (NIGEC) (South Central Regional Center, Tulane University, David Sailor, Director), through the US Department of Energy (Cooperative Agreement No. DE-FC03- 90ER61010). Any opinions, findings, and conclusions or recommendations expressed in this publication are those of the authors and do not necessarily reflect the views of the DOE. Dr Stephen L. Egbert graciously permitted the use of his 1992 high-resolution land use/land cover classification map.
References Atkinson, P.M., Cutler, M.E., Lewis, H., 1997. Mapping sub-pixel proportional land cover with AVHRR imagery. International Journal of Remote Sensing 18 (4), 917 /935. Box, E., Holben, B., Kalb, V., 1989. Accuracy of the AVHRR vegetation index as a predictor of biomass, primary productivity and net CO2 flux. Vegetation 80, 71 /89. Brown, J.F., Loveland, T.R., Merchant, J.W., Reed, B.C., Ohlen, D.O., 1993. Using multi-source data in global land cover characterization: concepts, requirements and methods. Photogrammetric Engineering and Remote Sensing 59, 977 /987. Davis, J.C., 1986. Statistics and data analysis in geology, second ed. Wiley, New York. DeFries, R., Hansen, M., Townshend, J., 1995. Global discrimination of land cover types from metrics derived from AVHRR Pathfinder data. Remote Sensing of Environment 54, 209 /222.
M.E. Jakubauskas et al. / Computers and Electronics in Agriculture 37 (2002) 127 /139
139
Doraiswamy, P.C., Cook, P.W., 1995. Spring wheat yield assessment using NOAA AVHRR data. Canadian Journal of Remote Sensing 21 (1), 43 /51. Egbert, S.E., Price, K.P., Nellis, M.D., Lee, R. 1995. Developing a land cover modeling protocol for the High Plains using multiseasonal Thematic Mapper imagery. Proceedings of the 1995 ASPRS/ACSM Annual Meeting, Charlotte, N.C. pp. 836 /845. Eidenshink, J.C., 1992. The 1990 conterminous US AVHRR data set. Photogrammetric Engineering and Remote Sensing 58 (6), 809 /813. Gutman, G., Ignatov, A., 1998. Derivation of green vegetation fraction from NOAA-AVHRR for use in numerical weather prediction models. International Journal of Remote Sensing 19, 1533. Gutman, G., Ohring, G., Tarpley, D., Ambroziak, R., 1989. Albedo of the US Great Plains as determined from NOAA-9 AVHRR data. Journal of Climate 2, 608 /617. Jakubauskas, M.E., Legates, D.R., Kastens, J., 2001. Harmonic analysis of time-series AVHRR NDVI data. Photogrammetric Engineering and Remote Sensing 67 (4), 461 /470. Jensen, J., 1996. Introductory Digital Image Processing. Prentice Hall, Englewood Cliffs, NJ, p. 316. Kastens, D.L., Price, K.P., Martinko, E.A., 1998. Assessing wheat conditions in Kansas using biweekly AVHRR datasets and crop phenological indices. First International Conference on Geospatial Information in Agriculture and Forestry, Orlando, FL, June 1 /3, 1998. Kremer, R.G., Running, S.W., 1993. Community type differentiation using NOAA/AVHRR data within a sagebrush-steppe ecosystem. Remote Sensing of Environment 46, 311 /318. Lee, R. 1999. Modeling Corn Yields in Iowa Using Time-series Analysis of AVHRR Data and Vegetation Phenological Metrics. PhD Dissertation, Department of Geography, University of Kansas, Lawrence, KS. Legates, D.R., Willmott, C.J., 1990. Mean seasonal and spatial variability in global surface air temperature. Theoretical and Applied Climatology 41 (1), 11 /21. Lloyd, D., 1990. A phenological classification of terrestrial vegetation cover using shortwave vegetation index imagery. International Journal of Remote Sensing 11 (12), 2269 /2279. Loveland, T.R., Merchant, J.W., Ohlen, D.O., Brown, J.F., 1991. Development of a land-cover characteristics database for the conterminous US. Photogrammetric Engineering and Remote Sensing 57, 1453 /1463. Loveland, T.R., Merchant, J.W., Brown, J.F., Ohlen, D.O., Reed, B.C., Olson, P., Hutchinson, J., 1995. Seasonal land-cover regions of the United States. Annals of the Association of American Geographers 85, 339 /355. Rayner, J.N., 1971. An Introduction to Spectral Analysis. Pion Ltd, London, p. 174. 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. Tucker, C.J., 1979. Red and photographic infrared linear combinations for monitoring vegetation. Remote Sensing of Environment 8, 127 /150. Yin, Z., Williams, T.H.L., 1997. Obtaining spatial and temporal vegetation data from Landsat MSS and AVHRR/NOAA satellite images for a hydrologic model. Photogrammetric Engineering and Remote Sensing 63 (1), 69 /77.