 |
|
|
||||
| USGS Earth Resources Observation and Science Center (EROS) |
J. Eidenshink, D. Hollaren, R. Vandersnick, C. VanBeek, D. Steinwand, C. Wivell, D. Meyer
The science requirements for the global land 1-km AVHRR data set are documented by IGBP, CEOS, NASA, CEC and other international research organizations. Among the requirements is the need for higher level products that are derived from the raw 1-km AVHRR data, such as vegetation indices and periodic temporal composites. The composites must include the minimum ten bands. The products must be generated from radiometrically calibrated, atmospherically corrected, and geometrically registered data. As a result, the typical user, who is mainly interested in the data for assessment of vegetation condition, derivation of biophysical parameters, mapping of land cover or monitoring land surface characteristics, is freed of the burden of preprocessing the raw data. The assumption here is that raw data have been processed using widely accepted, well defined and documented processing standards.
The definition of the processing standards for the prototype products has been reached by consensus as a result of the efforts of the IGBP-DIS land cover working group (Townshend et al.), CEOS, AVHRR Pathfinder land science working group (James and Kalluri) and l'Observatoire Sahara Sahel (Faizoun, 1992). The standards for the primary processing steps of radiometric calibration, atmospheric correction, geometric registration, and compositing have been documented and guidelines for other functions such as map projection, compositing period, and product format were provided. Overall, the recommendations of the IGBP were the most comprehensive and are being used to create the prototype products.
In a project of this magnitude, data processing efficiency is very important. Each day, the amount of data to be processed increases incrementally. The goal for any processing scenario is the ability to process a days worth of data in at least a day. Existing operational programs, such as the greenness mapping program at the EDC meet the requirement for regional and subcontinental data sets of the conterminous U.S., Eurasia, and Alaska. The AVHRR Pathfinder Program has developed the capability to process daily 4-km Global Area Coverage (GAC) data (James and Kalluri). But the capability to process daily 1-km AVHRR data on a global basis in a timely fashion does not currently exist at the EDC. Therefore, new processing technology will be developed and existing technologies will be refined to meet the requirement of processing the data in a reasonable time period.
Processing Methods
The implementation of the data processing flow is generally a stepwise process that incrementally applies higher order processing in a logical and efficient manner. The steps identified are:
The following is a discussion of the methods and parameters that have been recommended for each step and a review of the current status and plans for implementation of the processing algorithms.
Radiometric Calibration
Radiometric calibration of the AVHRR visible and near-infrared channels (channels 1 and 2) is difficult because there is not always reliable preflight calibration, no onboard calibration, and difficulty with inflight calibration. Preflight calibration coefficients can change while the instrument is in storage or after launch because of the space environment. Degradation of AVHRR sensors after launch is well documented (e.g., Rao, 1987; Price, 1987; Holben et al., 1990). These studies have used a variety of approaches such as ground-based measurements from stable sites such as homogeneous desert targets to monitor the degradation of the sensors. Although the results from the different approaches often are not in close agreement, there is an increasing consensus on suitable coefficients. The calibration method that is recommended for this data set accounts for sensor degradation by using coefficients developed by Teillet and Holben (1994). Their calculation uses a desert calibration to develop time-dependent calibration coefficients for the AVHRR sensor on NOAA-11. Use of calibration coefficients involves extrapolation of the most recent calibration results for processing data on a near real-time basis. Therefore, the time- dependent coefficients are based on a piecewise linear fit of the desert results. A piecewise linear fit is recommended for operational use because, unlike polynomial fits, they will not change retroactively when new data are added to the end of the time series.
The equations for radiometric calibration to radiance and reflectance for the visible and near infrared channels are described by Teillet and Holben (in press). The calibration coefficients for AVHRR thermal channels 3, 4, and 5 are derived onboard the satellite using a view of a stable blackbody and deep space as a reference. The calibration process converts raw digital counts to radiance as described by Kidwell (1991). The radiance values for all channels are stored with 16-bit precision.
Atmospheric Correction
The impact of atmospheric effects on the AVHRR channel 1 and 2 data and NDVI can be significant. Four principle atmospheric factors, water vapor, aEDCols, ozone, and Rayleigh scattering, are considered to have the most impact.
The corrections for ozone and Rayleigh scattering are straightforward (Teillet, 1991). Appropriate Rayleigh scattering correction must include an adjustment for topography. Recommended reference values for Rayleigh optical depths for standard pressure and temperature conditions are available (Teillet, 1990). The local elevation adjustment can be derived using ETOPO5 which is the best available global digital terrain data at this time and for this task has sufficient accuracy.
The correction for ozone should be based on actual measurements derived from the Total Ozone Mapping Spectrometer (TOMS) or other appropriate sensors. However, access and utilization of these data can be difficult. Using the concentration values from standard climatic tables with latitudinal and seasonal dependence is acceptable and is the approach that has been implemented for this data set.
Proper use of the radiative transfer code is important. Bandpass calculations of 0.005-micrometer spacing or better are recommended. No specific radiative transfer code is recommended but it should be noted that results from different code using large optical depths for aEDCols and off-nadir viewing angles greater than 60 degrees can vary significantly.
Several approaches for correction of water vapor exist but there is no community agreement on a feasible method. The basic problem is determination of the spatial and temporal variability of water vapor concentrations. The same circumstances affect aEDCol corrections. Therefore no water vapor or aEDCol correction will be applied.
The input to the atmospheric correction process is radiance values from the calibrated visible and near- infrared channels. The output of the atmospheric correction process is surface reflectance (in percent) of the visible and near-infrared, albeit without corrections for water vapor and aEDCols.
Computation of NDVI
The NDVI is the difference of near-infrared (NIR, channel 2) and visible (VIS, channel 1) reflectance values normalized over the sum of channels 1 and 2 ((NIR-VIS)/(NIR+VIS)). The NDVI equation produces values in the range of -1.0 to 1.0, where increasing positive values indicate increasing green vegetation and negative values indicate nonvegetated surface features such as water, barren, ice, and snow or clouds. To obtain the most precision, the NDVI is derived from calibrated channels 1 and 2 data in 16- bit precision, prior to geometric registration and resampling.
It is recommended to scale the computed NDVI results to 8-bit range to minimize the volume and optimize analysis and display. This is acceptable because there is very little quantitative understanding of the relationship between excessive mathematical precision in the NDVI and physical measurements of the vegetation condition. If greater precision is required, the NDVI can be computed from the calibrated channel 1 and 2 values that are stored with 16-bit precision.
The scaling method is chosen to emphasize different types of vegetation or vegetation condition. The EDC routinely uses a formula that scales the NDVI from -1.0 to 1.0 as 0 to 200, where each value represents 1.0 percent of the total possible range. However, the total possible range is never used. The Canada Centre for Remote Sensing and others use a scaling from -0.05 to 0.67 NDVI values for the 8-bit range (0-255) or values 0 to 200 representing NDVI values from -0.05 to 0.95, where each value represents 0.5 percent, values less than -0.05 are scaled to 0 and values greater than 0.95 are scaled to 200. The NDVI scaling method for this project has not yet been finalized. Further research is needed to determine a scaling method that adequately characterizes vegetation on a global scale.
Geometric Registration
Geometric registration involves precise transformation of the image from the sensor-based projection to an earth surface-based projection. This process includes calculating a satellite model, matching ground and image-based control points, and transformation and resampling the data to a map projection coordinate system. The satellite model is also used to compute satellite zenith, solar zenith, and relative azimuth viewing angles for each pixel.
Based on a comprehensive evaluation of map projections for global data sets (Steinwand et al., in press), the Interrupted Goode Homolosine is recommended for the data set. The Interrupted Goode Homolosine (Goode, 1925; Steinwand) has two important features. First, it is an equal area projection that facilitates spatial analysis. Second, it essentially divides the world into 12 regions that can be mosaicked into a global map. The regionalization of the global map has advantages for data handling.
The most important factor to consider in geometric registration is the positional accuracy of each pixel as it is moved from the sensor-based projection to the surface-based projection. The recommendation for this project is a positional accuracy of 1000 meters or less. A systematic correction using the satellite model alone cannot achieve this positional accuracy. Thus, it is necessary to perform control point matching and correction for terrain elevation. If terrain is not accounted for in the satellite's geometric model, registration errors of up to 12 kilometers can occur for extremely off nadir pixels in areas of high relief. Incorporating heights above the reference ellipsoid derived from digital elevation models into the satellite's geometric model can greatly reduce these registration errors. The digital elevation data used for this process is ETOPO5.
Most methods for control point matching use automatic correlation of image segments with ground control points and then integrate the adjustments derived from the correlation with the orbital model (Cracknell and Paithoonwattanakij 1989; Kelly and Hood 1991; Brunel and Marsouin 1987). The problem with this approach is the identification of enough control points on a global basis to ensure the required registration accuracy. The problem is complicated by the routine presence of clouds that prevents a successful correlation process. The problem is minimized by having an extremely dense and evenly distributed set of control points.
Control points for AVHRR data are commonly identifiable features along coast lines, lakes, and rivers. Some prominent physiographic features are used in arid regions where hydrologic features are sparse. One approach is to use satellite imagery to develop the control points. Control points are chosen from on near- nadir, cloud-free base image of 1-km AVHRR data. Another possible approach is to select the control points from Landsat Multispectral Scanner (MSS) data and then degrade the image to AVHRR resolution. Although the use of MSS data improved the precision of the control point marking process it is time consuming and expensive to develop.
Yet another approach is to use hydrologic features from vector data sets such as the Digital Chart of the World (DCW) and the World Vector Shoreline (WVS). The vector data representing coastlines, shorelines, and other major hydrologic features are rasterized to AVHRR resolution and warped into satellite projection. A comparable edge feature is derived from the AVHRR data at approximately the same location using edge extraction filtering methods. The correlation is performed on several control points and adjustments are made. The advantage of this approach is that DCW and WVS are available worldwide. The accuracy of the DCW, however, is variable, depending on the geographic location and original map source. The WVS is considered more accurate for shorelines than DCW, but DCW is the best source of inland water bodies and rivers.
The prototype procedure that was developed using the DCW had an accuracy of 0.8 rmse to 1.3 rmse pixels based on image to map verification and, therefore, meets the recommended standards in most cases. The accuracy of the registration was improved slightly by constraining the correlation process. However, the accuracy was still slightly greater than 1.0 rmse pixels. The most limiting factor is still the overall accuracy of the DCW and WVS. However, the important point was to achieve multi-temporal (image to image) registration rather than absolute positional accuracy in order to prevent "blurring" in the compositing process. Verification of image to image registration had an accuracy of less than 1 pixel rms.
Compositing
The first factor to consider in the compositing process is the length of the compositing period. Compositing periods of 7, 10, and 14 days have been used most commonly. The choice of the period is usually based on the length of time necessary to obtain a composite with minimal cloud contamination and or the amount of time necessary to observe meaningful changes in surface characteristics. The compositing period that is recommended for the prototype products is approximately 10 days created by month. Thus, January has three composites of 10, 10, and 11 days; February has 10, 10, and 9 or 8 depending on whether it is a leap year, and so on. This procedure has the advantage of creating calendar month composites, which is a common reporting period for agronomic and biophysical characteristics.
The recommended method is maximum NDVI compositing. The NDVI is examined pixel by pixel for each observation during the compositing period to determine the maximum value. The retention of the highest NDVI value reduces the number of cloud-contaminated pixels and selects the pixels nearest to nadir (Holben, 1986).
Other compositing techniques, including those using multiple criteria such as maximum NDVI and maximum apparent temperature, are currently being investigated. The results of these studies may provide a basis for a different compositing technique.
Processing Flow
The order of the processing steps can vary, depending on the desired types of products. The processing flow for the global land 1-km AVHRR 10-day composites was designed to take advantage of computer resources and provide optimum product flexibility.
The first step was to read in the portion of an orbital segment corresponding to the region of the Interrupted Goode Homolosine projection that was being processed. The satellite model was used to compute the three solar/satellite viewing angles. The angles are eventually used in the atmospheric correction process.
The next step was to perform the control point matching and terrain correction process and develop the adjustments to the orbital model. The control points and the terrain data were warped into the satellite projection for the correlation process. A transformation grid was computed but not applied in this step.
The calibration of the five AVHRR channels to radiance was completed next. The date of the observation was determined and the appropriate time adjusted coefficients for channel 1 and 2.
Next, the geometric registration was completed for the five AVHRR spectral channels and the three viewing angles. The transformation grid developed in the control point correlation process was applied. The grid postings were at 10-km intervals to be consistent with the ETOPO5 digital elevation data. Intervening pixel locations were interpolated.
The next step was to compute the NDVI and perform the maximum NDVI compositing. The output of the compositing process was defined as a 10-band image that includes the NDVI value for each pixel selected by the maximum value compositing for the 10-day period, the radiometrically calibrated channel 1-5 values, the satellite viewing geometry data, and a date index value. The date index band is provided to allow a user to identify the specific date and scene identification of the observation for each pixel in the composite period. Since there are multiple observations of the same ground location each day, especially in higher latitude regions, an ancillary table of the specific date and scene information was provided. The viewing geometry information was provided for several reasons including the identification of pixels that would be considered off nadir or for use in user defined atmospheric correction algorithms. The steps described so far were repeated on all orbital segments for the 10-day period.
The final step was atmospheric correction. Typically, atmospheric correction is performed in conjunction with calibration. However, in a recent study (Cihlar and Huang, personal communication) it was shown that using atmospherically corrected data in the maximum NDVI compositing process increased the probability of selecting pixels with higher satellite zenith angles, and preferably in the backscatter direction. Based on these findings and one other important factor, the atmospheric correction will be applied following the compositing process.
The rationale for applying atmospheric correction after compositing is feasible because the maximum NDVI value selection process yields the same results whether the NDVI is computed from raw digital numbers or calibrated radiance. The same pixel is selected from the array of choices. In this process, the NDVI used in the compositing process were computed from the calibrated radiance values that were subsequently be used in the atmospheric correction process. After the compositing process was completed, the channel 1 and 2 radiance values used in the computation of the NDVI that were selected as the maximum value were atmospherically corrected and converted to surface reflectance. The reflectance values were used to recompute the NDVI for the composite.
Perhaps the most significant advantage of this approach is that it provides an opportunity to apply improved atmospheric corrections to the data in the future. For example, there is currently no consensus method for water vapor or aEDCol corrections on a global basis. However, ongoing research and validation is likely to provide a suitable method at some point in the future. When this occurs, it will be possible to reprocess the archived radiance values and recompute a better atmospherically corrected NDVI. The computation will be completed without having to use extensive resources to calibrate, geometrically register, and composite the raw data.
The basic assumption is that the original maximum value compositing based on the uncorrected NDVI is acceptable. If the improved atmospheric correction procedures would eliminate the initial concern regarding selection of off nadir views or a new compositing process is developed, the argument for saving processing resources is mute. However, this process still offers the most flexibility given the current state of the art.
The First Prototype Products
The challenge of processing a global land 1-km AVHRR data set is great. One of the areas of concentration is the development of the 10-day composite product. A benchmark 10-day composite period was identified to test development of the processing system. The benchmark period chosen was June 21- 30, 1992. This is the summer/winter solstice and represents a period with extreme sun angles and spatial variation in global vegetation condition.
The first test was to create a single day 1-km global image. The date, June 24, 1992, was selected from the benchmark period. The objective of the test was to determine the feasibility of registering the half-orbit passes into the Interrupted Goode Homolosine projection using a systematic satellite model correction (i.e., without ground control points).
Using the Goode's Interrupted Homolosine at 1-km resolution, the global image is 17,347 lines by 40,031 samples (694 megabytes per band) for an 8-bit image. The uncalibrated VIS and NIR data of the 14 orbits were processed. The individual orbits were overlapped left to right to create the final global image. The band combination 2, 1, 1 as red, green, blue was used in creating a false color image.
The second test was to create a global 1-km 10-day composite for the benchmark period of June 21-30, 1992. The processing for the composite includes calibration, geometric registration, computation of NDVI, and compositing. No atmospheric correction was performed on the data. The calibration of channels 1 and 2 adheres to the recommended procedure that accounts for sensor degradation over time. The thermal channels were corrected using the standard procedure as described by NOAA (Kidwell, 1991). The NDVI values were computed from calibrated reflectance maintained in 16-bit precision. Afterwards, all five channels and the NDVI were converted to byte data to reduce the data volume using standard EDC procedures (Eidenshink, 1992). The geometric registration was performed using DCW for ground control points, and the data were transformed to the Interrupted Goode Homolosine projection. The data were composited using maximum NDVI compositing and a 10-band image, as described earlier, was created.
The input to the processing was the stitched orbital segments. There were 386 scenes stitched into 130 orbital segments for the 10-day period. The raw data volume for the 386 scenes was approximately 21 gigabytes. The problems associated with processing extremely large data volumes were minimized by splitting the world into the 12 regions defined by the Interrupted Goode Homolosine projection. The portion of each orbital segment that intersects a region was extracted and processed. As multiple orbital segments were processed for a region, they were composited. After the processing for all the regions was complete they were mosaicked to form the global image, and a land/sea mask was applied to the data. The 10-day composite is represented as a greenness image.
Band Description of Composite Images
| Band | Description | Band | Description |
|---|---|---|---|
| 1 | AVHRR channel 1 | 6 | NDVI |
| 2 | AVHRR channel 2 | 7 | Satellite zenith |
| 3 | AVHRR channel 3 | 8 | Solar zenith |
| 4 | AVHRR channel 4 | 9 | Relative azimuth |
| 5 | AVHRR channel 5 | 10 | Date Index |
| AccessibilityFOIAPrivacyPolicies and Notices | |
  |
|
