- Email: [email protected]
Topographic effects in AVHRR NDVI data
ELSEVIER
Topographic Effects in AVHRR NDVI Data D. W. Burgess,* P. Lewis,* and J.-P. A. L. Mulled A n investigation is made of the errors induced in the normalized difference vegetation index (NDVI) by topographic variation. The Advanced RAdiometric RAy Tracer (ARARAT) is used to simulate AVHRR imagery for various viewing and illumination conditions at 1.1-km and 50-m resolution. Topographic error is calculated as the difference between the variable terrain simulation and a flat plane simulation. Correlation between this error and various topographic factors such as shadowing and adjacement hill illumination are examined. The importance of sky occlusion as a topographic factor is highlighted. While topographic errors of 13.5% are found in the 50-m resolution simulations, these errors are reduced markedly in the 1.1-km dataset to approximately 3%.
INTRODUCTION Over recent years, the normalized difference vegetation index (NDVI), derived f~om Advanced Very High Resolution Radiometer (AVHRR) data, has become the de facto standard for monitoring vegetation growth on a global seale. It is now recognized in the literature that a number of factors lead to variations in the NDVI. These include atmospheric effects (Cihlar et al., 1994), soil brightness and canopy density (Huete et al., 1985; Borel and Gerstl, 1994), and sensor degradation (Gutman, 1991). It has also been recognized that the NDVI is not invariant under different viewing and illumination geometries and a number of empirical bulk corrections have been proposed to compensate for the surface anisotropy (Gutman, 1991; Goward et al., 1991; Taylor et al:, 1985).
*Manaaki-Whenua-Landcare Research, Wellington,New Zealand tDepartment of Geography,UniversityCollege London, United Kingdom *Department of Photogrammetryand Surveying,UniversityCollege London, United Kingdom Address correspondence to D. W. Burgess, Manaaki-WhenuaLandcare Research, Box 38-491, Wellington,New Zealand. Received 10 August 199.4;revised 21 June 1995. REMOTE SENS. ENVIRON. 54:223-232 (1995) ©Elsevier Science Inc., 1995 655 Avenueof the Amerieas, New York, NY 10010
A more complex effect on the NDVI arises due to topographic variation. Encompassed in this issue are the effects of shadow, adjacent hill illumination, sky occlusion, and slope orientation with respect to the bidirectional reflectance distribution function (BRDF) of the cover type (Schaaf et al., 1994). This article presents results of a series of simulation studies using the Advanced RAdiometric RAy Tracer (ARARAT) (Lewis and Muller, 1992), which simulates the radiation regime, for a given set of illumination and viewing conditions, using 3-D models of individual or fields of plants with associated lower geometric boundary conditions described by a digital elevation model (DEM) which is used to represent soil roughness or macroscopic terrain effects. Simulations over AVHRR wavebands and viewing and illumination geometries have been carried out both at the nominal 1.1-km resolution (nadir viewing) and at a higher resolution of 50 m (nadir viewing), the latter allowing topographic effects to be more easily identified.
BACKGROUND In general, successful NDVI studies are carried out in large, relatively fiat parts of the world such as the extensive cropping areas of the midwest of the United States, or areas of Africa. This type of terrain is ideal, given that the resolution of the data at nadir is relatively low at 1.1 km, and a large fiat test site with contiguous ground cover allows statistically significant measurements to be made. However, the consequence of this is that the issue of the effect of topography on the NDVI measure has largely been ignored in the literature. If there is a dependence on topography, then there are many terrestrial areas where it is not reasonable to ignore the effect.
Radlometric Topographic Effects The topographic effect is defined as the variation in radiance from inclined surfaces compared with radiance from a horizontal surface, as a function of the orientation 0034-4257 / 95 / $9.50 SSDI 0034-4257(95)00155-T
224 Burgess et al.
of the surface relative to the light source and sensor position (Holben and Justice, 1980). In other words, radiometric topographic effects can all ultimately be described as resulting from the "sun-sensor--target" geometry. However, to fully address the impact of topography, it is necessary to consider neighborhood effects alSO. If the earth's surface were an isotropic reflector, giving the same surface-leaving radiance regardless of viewing angle, then variations in the radiance received at the satellite would be simply a combination of atmospheric effects on the upwelling radiance and illumination angle effects. However, it has been shown that this is not the case (Smith et al., 1980; Gutman, 1991). A BRDF is needed to characterise the behavior of direct illumination on, and reflectance of, a given cover type, dependent on the incident and exitant angles of the radiance, as well as the radiometric and geometric properties of the land cover. In practice, this means that the solar zenith angle, view zenith angle, view azimuth, and terrain slope all affect the radiance measured for a particular ground cover. Therefore, the terrain slope at a given point is responsible for the point topographic effect. A less obvious effect of topography has been investigated by Teillet and Staenz (1992) with reference to AVHRR and simulated Moderate Resolution Imaging Spectroradiometer (MODIS) imagery. They studied the percentage change in NDVI measurements due to the different atmospheric path lengths resulting from terrain elevation, and showed that, with 23-km visibility, the change in AVHRR NDVI was 0.79% per 100-m elevation, for standard vegetation. They also showed that the change in NDVI was greater for less vegetated targets. Other topographic effects arise due to the context of the point within the surrounding terrain. These could be termed neighborhood topographic effects. The most obvious example is the result of direct illumination by the sun, causing shadows. The magnitude of this effect in raw satellite data is dependent on the solar elevation, decreasing as the Sun approaches zenith. Harder to quantify and model are the effects of diffuse illumination on sloping surfaces. The proportion of diffuse illumination is dependent on the turbidity of the atmosphere including the presence of cloud. Proy et al. (1989) identified two major interactions of diffuse illumination and topography. First, in steep terrain, pixels are often hidden from part of the sky hemisphere by adjacent terrain. This reduces the total diffuse irradiance reaching the point and hence affects its radiance. Second, terrain can be illuminated by reflected radiation from adjacent slopes, which has the effect of falsely increasing the reflected radiation from that point. Proy has shown that, while these diffuse irradiance effects may be negligible for directly illuminated pixels, they
may represent up to 33 % of the illumination contribution of shadowed pixels. Previous Work There have been several investigations made into the effects of topography on high resolution satellite data (e.g., Smith et al., 1980; Hall-K6nyves, 1987; Newton et al., 1991). Studies have also been made in controlled environments using hand-held radiometers (Holben and Justice, 1980; Holben and Justice, 1981; Pinty et al., 1990). In all cases, measurable effects were identified. While high resolution studies are invaluable for providing an insight into the mechanisms governing how topography affects the radiance sensed, for lower resolution data such as AVHRR, practitioners have tended to take a more statistical approach. For a number of years it was believed that the NDVI was invariant under different illumination and viewing conditions, and it has only been recently that researchers have started to consider that this may not be a valid assumption. Gutman (1991) has examined the frequency of individual view angles contributing to 10-day NDVI maximumvalue composite images over 1 month. Contrary to previous theories, he found that the NDVI composite was not biased towards nadir views, with their shorter atmospheric path lengths, but towards forward scattering views. This adds weight to the argument that the effects of surface anisotropy are carried through to the NDVI measures. He goes on to develop a simple empirical correction for AVHRR reflectance data based solely on view zenith angle. This had the effect of reducing the coefficient of variation of the NDVI by a factor of 2 from 0.12 to 0.06, where the mean corrected NDVI was 0.49. Cihlar et al. (1994) examined the BRDF effects in AVHRR data on NDVI compositing. They found them to be significant in four distinct cover types but noted that the anisotropy is reduced when data are atmospherically corrected, which it to be expected since the atmospheric correction process moderates the BRDF effects of the atmospheric phase function. METHOD
Gathering an adequate data set for investigations in this field is inherently very difficult, due to the fact that good sampling of the illumination and viewing geometry is required within the constraints of orbital characteristics, cloud cover, variation in sensor characteristics across the NOAA series, and atmospheric effects. This study has therefore taken advantage of the unique capabilities provided by a software simulator called the Advanced Radiometer Ray Tracer (ARARAT) for simulating AVHRR data. AVHRR Band 1 and Band 2 data have been simulated, with the resulting images yielding "at ground"
Topographic Effects in A VHRR NDVI Data 225
radiance measures. These radiances have then been converted to NDVI using the usual formula: (b2- bl)/ (b2 + bl). Although this step is more usually carried out on data that has been processed to reflectances, it has been shown that the relationship between radiancebased and reflectance-based NDVI measures is a simple one. It results from the multiplicative effect of the near infrared/visible ratio of' solar spectral irradiance, and can be approximated to a linear relationship over a NDVI range of 0.0-0.8 (Goward et at., 1991). While this does lead to a modest decrease in the range of the NDVI, 1 this was not considered significant for the generalized comparisons made in this study. The first series of simulations have been designed to look at the magnitude of variation possible due solely to varying Sun-sensor-target geometry. For this purpose, a flat plane DEM has been used, and the solar and view zenith angles have been varied. In this case the Ahmad-Deering BRDF model for prairie grass in June has been used (Ahmad and Deering, 1992). The second stimulation focuses on higher order effects of topography. A 12.5-m resolution DEM has been used with the simulation performed at 50-m resolution. Although considerably finer than actual AVHRR resolution, this detail has been necessary to identify the higher-order mechanisms of topographic effect. Because of the "bias" introduced by a BRDF, these simulations have first been run using a Lambertian reflectance function; then another simulation has been performed over the same area using the previously mentioned BRDF. Finally, a simulation has been performed at the AVHRR nominal resolution of 1.1 km over a DEM. An investigation of the correlation between this result and actual AVHRR data has been made. The variation in apparent topographic effect between this and the previous simulation has been examined. An Overview of ARARAT The first use of Monte-Carlo ray-tracing for satellite sensor data [Landstat Thematic Mapper (TM)] simulation was described by Muller and Dalton (1988). ARARAT is a second generation version which more correctly models the physics and can simulate a wider range of geometric primitives. ARARAT is described in detail in Lewis and Muller (1992). The aim here is only to give an overview of the functionality of the system. The basic principle is that of a Monte Carlo ray tracer (Cook et al., 1984; Kajiya, 1986; Glassner, 1989). Geometric models are rendered for a given set of illumination and viewing conditions. The sensor is characterized by a location and a pointing vector, for a given instanta-
l A radiance NDVI range of [0.0-0.8] corresponds to a reflectance NDVI range of [0.18-0.88J; thus the extent of the range is reduced from 0.8 NDVI units to 0.7 NDVI units.
Figure 1. ARARATsimulation of a millet field.
neous field of view. A modulation transfer function for the optics and finite aperture effects can be incorporated, but is not used in these simulations. The irradiance field is modeled by a gridded angular "map" of spectral sky radiance distribution and a (direct) solar spectral irradiance component. The original design focus of ARARAT was the simulation of ground based radiometer measurements. For this purpose, 3-D models of plants can be rendered, often with a DEM to represent the underlying soil structure. Figure 1 is an example of the kind of image that ARARAT produces at this scale. Because of the detailed physical basis to the simulation process, the software can be used at any spatial scale. Recent work by some of the authors, on the NASA First ISLSCP Field Experiment (FIFE) test site, employed ARARAT to simulate MODIS and AVHRR imagery. For satellite image simulations, account must be taken of the effects of atmospheric attenuation and scattering on both the directional irradiance function and the path from the ground to the sensor. Only the former are modeled in ARARAT at the present time; therefore, ARARAT simulates at-ground radiances. In order to simulate the data from different sensors, over
226 Burgess et al.
Imaging G e o m e t r y at Wellington, N Z ~"
~
'
~
0
Relative Azimuth Angle
AVHRR NOAA-9 L a L - 41 S L o n g . - 175
15 D e c e m b e r 2 January
270 s i ' 255 \
~
105
135
225
~
r
Solar Zenith A n g l e I
50.00 - 55.00
I 1
4s.00- 50.00 4o.00- 45.oo
1
3s.00- 4o.oo 30.00 - 35.00 25.00 - 30,00
J
Figure 2. X-SATVIEW (Barnsley et al., 1994) plot of the
summer imaging geometry of the NOAA-9 AVHRR satellite for Wellington, New Zealand.
the various spectral bands, ARARAT performs Monte Carlo sampling of a broadband sensor spectral response function (Lewis and Muller, 1992) as part of the rendering process. The simulator incorporates the Zibordi and Voss (1989) sky irradiance model which takes such inputs as ozone concentration, water vapor concentration, and aerosol scattering albedo. It also incorporates a number of bidirectional reflectance models including the Ahmad and Deering model (Ahmad and Deering, 1992) and the Pinty-Verstraete-Dickinson (PVD) model (Pinty et al., 1990). ARARAT uses reverse ray tracing, which is more computationally efficient than forward ray tracing for simulating satellite imagery (Muller and Dalton, 1988). The history of a ray is traced from the sensor to the surface, then out to the illumination source. Rays are allowed to interact with surrounding surfaces as they are reflected, until they finally escape the model to the sky. Direct illumination sampling is achieved by firing a ray directly at the sun and diffuse illumination sampling is done by firing sample rays over a hemisphere around the local surface normal according to the local surface reflectance function. Diffuse contributions from adjacent slopes are implicitly included in this algorithm. Not surprisingly, this process is computationally very expensive.
The Study Domain: Geometry A study was made of the viewing and illumination geometry of AVHRR data captured over New Zealand. Currently, afternoon passes are received daily from the NOAA-11 satellite, for the purpose of NDVI monitoring. Figure 2, which was generated using the package X-SATVIEW (Barnsley et al., 1994), shows the nominal geometry for an afternoon NOAA satellite, in midsum-
mer, for the data captured over central New Zealand. Because of the orbital characteristics of the satellite, this geometry does not vary greatly up and down the country. The relative azimuth angle is the angle between the .view and solar azimuth. From this figure, it is clear that the majority of images are very close to the principal plane, with relative azimuth angles close to 0 ° or 180 °. The solar zenith angle falls in the range of 25-55 ° , with the highest solar elevations (smallest solar zenith angles) occurring in the forward scatter direction. In midwinter, the data move out of the principal plane and are oriented along the 135 ° line. The solar zenith angle drops to a range of 65-85 ° . It should be noted that, while the viewing geometry varies significantly across a full AVHRR swath in a given scene, the variation over a small region of interest is small when compared to inter-scene variations. One restriction of the sky model used by ARARAT is that it models the atmosphere as parallel planes. This means that the simulator cannot be used for solar zenith angles beyond 70 ° , for which the model would not accurately reflect the true optical depth. For this reason the solar zenith angles investigated have been restricted to 25-55 ° - the lower limit being the summer limit for New Zealand and the upper limit chosen to be well within the acceptable simulation limits and to provide useful illumination levels. As indicated previously, the simulations are very time-consuming to run, so that further restrictions have been placed on the domain being investigated. Holben and Justice (1980) found that the greatest errors occurred in the principal plane. Since the summer imagery is oriented near the principal plane, and covers the solar zenith angle range already selected, only the principal plane situation is studied.
The Study Domain: Topography A site was chosen near the east coast of the North Island of New Zealand. It is hilly terrain, which has gradients ranging from 0 ° to 49 ° with an average of 20 °. The predominant ground cover is rye grass and clover, which is used for pastoral farming. In the 5 km x 5 km square site used for the high resolution simulation, the altitude ranges from 40 m to 500 m above sea level with a mean altitude of 195 m. Figure 3 shows the DEM used for this test area, which has a resolution of 12.5 m. The 22 km x 22 km area of the low resolution simulation (Fig. 9c) ranges in altitude from 21 m to 620 m above sea level with a mean elevation of 189 m. The average gradient is 16 ° within a range of 4-64 °. The digital elevation model used in this case was at 100 m resolution.
Simulating Imagery The first step in simulating the imagery was creating a sky model. In this case, a sky model appropriate to New
Topographic Effects in AVHRR NDVI Data 227
80.0
Band 2
¸%¸¸:¸ ~--~ 6 0 . 0
8 c
40.0
0--. • "'(3
@. /
-\'@. A. -.~. -.~--0 ~.~.Z~.~x_ -
.~ ~<
~..\-~
~.
.
sza=25 ®G---®-:E3 sza-~o ~" , ~- --~ i -4> sza=35 / u, / YJ- - - ~ ~ - ~ sza=40 / /~ ~ ~ sza--45 . / 7"< ~I~ ~ sza=SO
/~x- . & .
.~-A- - --A-- - -~
p<-
/"
~ .
.
.,~-~
/
, [,
]:>
./
®
¢n
Band 1 20.0
o'%.0
~;.o
o.o
+o.o
-100.0
View Zenith Angle (degrees)
Forwardscatter - Backscatter
Figure 3. 5 km x 5 km DEM of test area for simulation series 2 (12.5-m resolution).
Zealand atmospheric conditions was required. The local precipitable water vapor concentration was taken from tephigram data received at the meteorological station at Paraparaumu Airport near Wellington. The ozone concentration coefficient over New Zealand was obtained by interrogating the online database of Total Ozone Mapping Spectrometer (TOMS) data at the Geophysical Data Facility at Rutherford Appleton Laboratory, Oxford, United Kingdom. Appropriate coefficients of the/~ngstr6m formula were obtained from ~ngstr6m (1961) and others were set to average maritime values. The behavior o f the simulator can be customized in a number of ways. Because it contains an element of random sampling, results can be different from one simulation to the next. Therefore, it is vital to set the simulation parameters so that the results are convergent. The principal way to achieve this is by increasing the number of rays fired per pixel. After trials, it was found that 100,000 rays/pixel gave satisfactory results, with agreement in the output radiance for successive simulations to four significant figures. Another important parameter is the ray tree depth. This controls the number of intersections of the ray that are followed. As the signal is attenuated at each bounce, its contribution eventually becomes microscopic, if it does not terminate by escaping the terrain to the sky. However, it is imporatant to keep track of small contributions since their cumulative effect can be substantial. Experience has shown that, in general, a ray tree depth of 10 is adequate. This was verified during the simulations performed by monitoring the illumination contributions at various depths of the ray tree to ensure that the percent contributions at level 10 were small. These contributions ranged from 0.0005% down to < 0.000001%. Another customizable feature of ARARAT is the range of BRDF models that can be utilized. For the
Figure 4. Flat plane simulated radiances for June prairie grass for a range of solar zenith angles.
purposes of this study, the Ahmad-Deering model (Ahmad and Deering, 1992) was used with a cover type of prairie grass. This model was considered the best available to date and the inclusion of prairie grass parameters made it undoubtedly the best option for simulating a pasture type of cover. A Lambertian surface with a "green vegetation" reflectance was also used in a number of simulations. RESULTS Simulation Series 1--Flat Plane with a BRDF The simulations of AVHRR sensor response over a flat plane were carried out for solar zenith angles of 2555 ° in 5 ° steps, over a range of view zenith angles from 60 ° backscatter to 60 ° forwardscatter in steps of 15 °. The Ahmad-Deering (1992) June prairie grass BRDF model was used. In each case, the hotspot was also stimulated. The results are illustrated in Figure 4. It can be seen from this graph that both viewing and illumination angle have a large impact on "at ground radiance." Variations of up to 50% are apparent across different solar zenith angles. Figure 5 shows how this variation is propagated through to the NDVI measure, where the NDVI ranges from 0.47 to 0.68. Here we find that the backscatter hotspot becomes a minimum and the maximum NDVI values occur in the forwardscatter direction. This is consistent with the findings of Gutman (1991) in his study of the viewing geometry of pixels contributing to 10-day NDVI maximum value composites over the U.S. Great Plains. The variation in NDVI across the range of view zenith angles is approximately 34%, and the greatest variation across the solar zenith angles, at nadir, is approximately 11%. These figures have significant
228 Burgess et al.
0.65
0.60
I z
0.55
sza=26
0.50
sza=30 sza=35 sza=40 sza--45 sza=50 sza=55
Figure 6. NDVI image of simulation 2 with a
0.45
solar zenith angle of 25 ° . Area corresponds to 1 km x 5 km at right edge of DEM. 0"4000.0
i 50.0
0.0
-50.0
-100.0
View Zenith Angle (degrees) Forwardscatter- Backscatter
N D V I % _ error(x,y) Figure 5. Flat plane NDVI variation for June prairie grass for a range of solar zenith angles.
implications when considering the role of the NDVI measure in land cover monitoring. Even restricting the viewing geometry to + 30 °, NDVI variation is of the order of 24%. Simulation Series 2 - - A c t u a l Terrain with a
Lambertian Surface The second set of simulations was performed over actual terrain using the 5 km x 5 km DEM described previously, with a Lambertian surface reflectance model. The processing overhead with this series was far greater than for the first for two reasons. First, a 100 x 20 pixel image strip down the DEM was simulated rather than a single point, as in the flat plane case. Second, the processing of a single point is far more complex in varying terrain, since many interactions occur during ray tracing with adjacent hill facets. Because of this greater overhead, it was only possible to simulate three solar zenith angles (25 °, 40 °, and 55 °) for a nadir viewing position. "AVHRR" pixels (50 m) were generated over the 12.5-m DEM to produce a 100 × 20 image corresponding to a strip down the right-hand side of the test area shown in Figure 3. Corresponding flat plane simulations were performed for each terrain simulation to provide an NDVI measure with no terrain effect. These were essentially the same as the first series except that they used the Lambertian surface rather than a BBDF. The NDVI was calculated across each terrain simulation and compared to the flat plane NDVI. Figure 6 shows the resulting NDVI image for the case of a solar zenith angle of 25 °. Errors were expressed as a percentage of the flat plane values as follows:
(INDVIt~,,ai,(x,y) - NDVI~t(x,y)I / x 100, =\ NDVI~a,(x,y) /
(1)
for each pixel (x,y) in the scene. The mean scene N D V I % _ e r r o r was then calculated for each solar zenith angle. Table 1 shows that these errors were consistently in the region of 13.5% regardless of solar zenith angle. It must be appreciated that these figures are not a representation of total NDVI error due to viewing and illumination geometry and terrain effects but solely the latter. The former effect is eliminated by normalizing the error to the fiat plane generated for the same viewing conditions. These results were compared to a number of terrain factors in an attempt to establish the cause of this error. As part of its output, ARARAT generates a list of ray tree depth contributions. This means that, for each level in the tree of ray interactions, information is gathered on the diffuse and direct illumination contributions incident on the terrain surface at that point. From the first level of this tree we get the percent direct and percent diffuse illumination directly coming to the pixel. Therefore, the percent direct illumination is a measure of how much the pixel is in shadow. Table 1 shows that for Band 2 there was a strong negative correlation between percent direct illumination and percent NDVI error, which implies that shadowed pixels gave greater NDVI error than directly illuminated ones. From all the levels below level 1 in the ray tree we get the illumination information for the adjacent hill points that reflect light onto the target pixel. By adding up all these contributions, it is possible to obtain a measure of the total percent adjacent hill illumination. The extent of adjacent hill illumination is largely due to the complexity of the terrain. Steep sided valleys will cause light to be reflected from one side to the other. Table 1 shows that percent adjacent hill immuniation is highly correlated with NDVI error.
Topographic Effects in A VHRR NDVI Data 229
Table i. Results from Second Set of Simulations Correlation of % NDVI Error with: % Illumination Solar Zenith Angle
Mean % NDVI Error
% Sky Visible
25
13.5
- 0.99
0.99
0.99
- 0.62
- 0.98
40
13.4
- 0.98
0.98
0.98
- 0.48
- 0.92
55
13.5
- 0.91
0.95
0.93
- 0.61
- 0.82
Of greatest interest in determining significant terrain factors was the correlation of the percent NDVI error with a percent sky visible image. This image was created also using ARARAT but with no sun, a uniform sky hemisphere and a perfectly white diffuse reflectance function. Under these conditions, the resultant image is equivalent to the percent of sky visible at each point. Figure 7 shows this image. The minimum percent sky visible value is 19.8%; the maximum is 99.6% and the mean is 68.2%. In Table 1, it can be seen that this image was consistently highly correlated with the percent NDVI error image. Obviously this factor is related to diffuse illumination, however, the beauty of it, as a potential tool for normalizing NDVI data for terrain effects, is that it is viewing and illumination geometry independent. Neither the Band 2 direct illumination measure nor the percent adjacent hill illumination measure have this advantage. Simulation 3 - A c t u a l Terrain Using a BRDF Having identified, first, the variation in NDVI caused by viewing and illumination geometry over a fiat plane using a BRDF, and second, the effect of actual terrain on the NDVI when a Lambertian surface is assumed, the natural progression is to look at what happens when a BRDF is used in an actual terrain simulation. Utilizing a BRDF leads to further processing overheads, which Figure 7. Percent sky visible image over area corresponding to the DEM in Figure 3.
Adjacent Hill Band 1 Band 2
Direct Band 1 Band 2
unfortunately made multiple simulations in this series impractical. The one completed simulation of the 100 x 20 pixel strip took approximately i week to complete, running in parallel on 18-24 Sun and Silicon Graphics Workstations. This simulation used the same BRDF as in the first fiat plane series and was done at a solar zenith angle of 55 ° and a view zenith angle of - 5 0 °, that is, near the hot spot. The resulting mean percent NDVI error was 16.1%, which was greater than in the second simulation (13.5%) but not by as much as might be expected in light of the results of the first fiat plane results. However, closer examination of the nature of BRDFs in general, and in particular the one utilized in this simulation, reveals why the difference in error is not great. As the slope orientation in hilly terrain varies
Figure 8. Prairie grass BRDF at 662 nm showing terrain variation line.
Praire Grass BRDF X= 6 6 2 p r n
0.065 0.06 0.055 0.05 ~ o114 3~ 0.0.003~
5
0.025~O "40 5
•
: sza
- v z a = 5°
230
B u r g e s s e t al.
(a)
Simulation 4--Actual Terrain at l.l-km Resolution Having established that various topographic effects appear to have an impact on the measured NDVI in 50 m resolution data, the next step is to see if this variation is still present at the scale of AVHRR data. An area of 20 x 20 pixels was simulated (22 km x 22 km) over the test site as indicated previously, for a solar zenith angle of 55 ° with a nadir viewing position and a Lambertian green surface. Band 2 of the 1.1-km resolution resulting image is shown in Figure 9. A coarser resolution DEM of 100 m was used in this simulation, which was considered appropriate given the desired simulation output resolution of 1.1 km and the processing overheads of using the finer, 12.5-m resolution DEM over such a large
(b)
area.
(c) Comparison of actual and simulated data: a) Simulation 4, Band 2; b) AVHRR data, Band 2; c) DEM of area for a) and b). F i g u r e 9.
over a small geographical area, both the effective viewing and illumination angle change with respect to the surface normal, but the angle between the sun and viewing direction does not. Figure 8 shows the AhmadDeering model for Prairie grass at 662 nm (Ahmad and Deering, 1992). The line of triangles indicates the case where the view zenith angle plus the solar zenith angle is equal to 5 ° . This is the case in the third simulation with a solar zenith angle of 55 ° and a view zenith angle of - 5 0 °. It can be seen that as the terrain varies, the BRDF is sampled along a line: vza = - sza + a,
(2)
where v z a is the view zenith angle, s z a is the solar zenith angle, and a is a constant representing the angle between the v z a and s z a . In general, the variation in the BRDF along this line is less than the variation over view zenith angles for a given solar zenith angle; therefore, variations in overall viewing geometry across an AVHRR scene and between distinct scenes will have a greater effect on the measured radiance than variation caused by slope orientation.
Again, flat plane values were also simulated and the resulting NDVI values compared. In this case the average percent difference between the terrain and flat plane NDVI values was approximately 3%. This compares with the 13.5% variation of simulation 2 and the 16.1% variation of simulation 3. In light of the lower errors in the coarser resolution simulation, it was speculated that perhaps this was due to the coarser resolution DEM, and not a direct result of larger pixel size. A smaller 10 km x 10 km subscene was resimulated using the high resolution 12.5-m DEM. An increase of approximately 1% in the mean percent NDVI error was found. This suggests that there is only a slight sensitivity to the DEM resolution used. Therefore, it would appear that, at the coarser resolution, the effects of topography on the resulting radiance are moderated since, within one pixel area of 1.1 km 2, a larger sample of slopes and illumination conditions are represented. An attempt was made to correlate the simulated data with actual data received over the same area. Figure 9b shows actual AVHRR data corresponding to the simulation in Figure 9a. Figure 9c shows the DEM of the same area. Similarities can be seen in the flatter area in the bottom left quandrant; however, the overriding impact of the dark area of native forest on the southeastern face of the ridge in the top left quadrant makes it hard to discern any correspondence in this area. Efforts to correlate these datasets were therefore largely unsuccessful. It is clear that the nominal contiguous pasture cover of the test site actually contained a significant amount of variation in the form of shelter belts and small stands of native bush. Where these areas were clearly indicated on a 1:50,000 scale map, for example, an area of bush mentioned previously, they were masked out, along with a few pixels of ocean in the lower right-hand corner of the scene. However, correlation over the unmasked areas was still not possible. This was attributed to small scale variation in the cover characteristics, which were not modeled by the simulator and difficulties in accurate subpixel registration of the datasets.
Topographic Effects in A VHRR NDVI Data 231
CONCLUSIONS Valuable insights have been gained into the mechanisms of topographic effects on the NDVI measure. The fiat plane series of simulations highlighted the significance of the variation in measurement due solely to viewing and illumination geometry. It also showed that the variation in viewing geometry had a greater effect than the sun position. The second series of simulations yielded a topographic error of 13.5 % in the NDVI values irrespective of illumination angle for 50-m resolution pixels. Shadowing, adjacent hill illumination and skT occlusion were well correlated to this error, suggesting that these mechanisms contribute to the variation. The impact of the BRDF in complex terrain was briefly investigated. It was found that, while surface anisotropy does add to the variability in the measured radiance, the constraints on the sampling of the BRDF imposed by the geometry of the system means that the variation is generally not as great as that caused by view angle variations for a given illumination angle. The fourth simulation showed that the previously identified NDVI errors were significantly reduced in 1.1-kin resolution data (3%). When this error is compared to the overall variation in NDVI due to Sunsensor geometry demonstrated in the first series of simulations (34% across all view zenith angles and a maximum of 11% across solar zenith angles), it is clear that the "bulk effect" of viewing and illumination geometry, irrespective of terrain orientation, is the dominant source of error in 1.1-kin resolution data. We suggest that, f~r the operational correction of AVHRR NDVI data, it is reasonable to ignore high-order topographic effects such as sky occlusion and adjacent hill illumination; however, bulk corrections should be made for viewing and illumination geometry variations across the full AVHRR swath in order to accurately compare NDVI measurements in different parts of the scene and in other scenes. Funds for this study were provided by The New Zealand Foundation for Research, Science and Technology under Contract C09208. D. Burgess wishes to acknowledge the support of the British Council, with whose help a significant portion of this work was undertaken in the United Kingdom, and Professor Harley and the Department of Photogrammetry and Surveying at University College London for support during the 6-month sabbatical. REFERENCES
Ahmad, S. P., and Deering, D. W. (1992), A simple analytical function for bidirectional reflectance, J. Geophys. Res. 97(D17):18,867-18,886. AngstrOm, A. (1961), Techniques of determining the turbidity of the atmosphere, Tellus 13:214-223.
Barnsley, M. J., Strahler, A. H., Morris, K. P., and Muller J.-P. (1994), Sampling the surface bidirectional reflectance distribution function (BRDF): 1. Evaluation of current and future satellite sensors, Remote Sens. Rev. 8:271-311. Borel, C. C., and Gerstl., S. A. W. (1994), Nonlinear spectral mixing models for vegetative and soil surfaces, Remote Sens. Environ. 47:403-416. Cihlar, J., Manak, D., and Voisin, N. (1994), AVHRR bidirectional reflectance effects and compositing, Remote Sens. Environ. 48:77-88. Cook, R., Porter, T., and Carpenter, L. (1984), Distributed ray tracing, Comput. Graphics (Proc. SIGGRAPH) 18(3): 137-145. Glassner, A. (1989), An Introduction to Ray Tracing, Academic New York, 1989. Goward, S. N., Markham, B., Dye, D. G., Dulaney, W., and Yang, J. (1991), Normalized difference vegetation index measurements from the Advanced Very High Resolution Radiometer, Remote Sens. Environ. 35:257-277. Gutman, G. G. (1991), Vegetation indices from AVHRR: an update and future prospects, Remote Sens. Environ. 35: 121-136. Hall-K6nyves, K. (1987), The topographic effect on LANDSAT data in gently undulating terrain in Southern Sweden, Int. J. Remote Sens. 8(2):157-168. Holben, B. N., and Justice, C. O. (1980), The topographic effect on spectral response from nadir-pointing sensors, Photogramm. Eng. Remote Sens. 46(9):1191-1200. Holben, B., and Justice, C. (1981), An examination of spectral band ratioing to reduce the topographic effect on remotely sensed data, Int. J. Remote Sens. 2(2):115-133. Huete, A. R., Jackson, R. D., and Post, D. F. (1985), Spectral response of a plant canopy with different soil backgrounds, Remote Sens. Environ. 17:37-53. Kajiya, J. T. (1986), The rendering equation, Comput. Graphics 20(4):143-150. Lewis, P., and Muller, J.-P. (1992), The Advanced Radiometric Ray Tracer: ARARAT for plant canopy reflectance simulation, in Proceedings of the XVI~ Congress of the ISPRS, Washington, DC, Commission 7, Vol. XXIX, pp. 26-34. Muller, J.-P., and Dalton, M. (1988), Application of ray-tracing to satellite image understanding, in Proceedings oflGARSS "88 Symposium, Edinburgh, Scotland, 13-16 September, pp. 1187-1190. Newton, A., Muller, J.-P., and Pearson, J. (1991), SPOT DEM shading for Landsat-TM topographic correction, in Proceedings of IGARSS "91, Espoo, Finland, 3-6, June Vol. 2, pp. 655-659. Pinty, B., Verstraete, M. M., and Dickinson, R. E. (1990), A physical model of the bidirectional reflectance of vegetation canopies. 2. Inversion and validation, J. Geophys. Res. 95(D8):11,767-11,775. Proy, C., Tanrd, D., and Deschamps, P. Y. (1989), Evaluation of topographic effects in remotely sensed data, Remote Sens. Environ. 30:21-32. Schaaf, C. B., Li, X., and Strahler, A. H. (1994), Topographic effects on bidirectional and hemispherical reflectance calculated with a geometric-optical model, IEEE Trans. Geosci. Remote Sens. 32(6):1186-1193. Smith, J. A., Lin, T. L., and Ranson, K. J. (1980), The Lam-
232 Burgess et al.
bertian assumption and Landsat data, Photogramm. Eng. Remote Sens. 46(9):1183-1189. Taylor, B. F., Dini, P. W., and Kidson, J. W. (1985), Determination of seasonal and interannual variation in New Zealand pasture growth from NOAA-7 data, Remote Sens. Environ. 18:177-192. Teillet, P. M., and Staenz, K. (1992), Atmospheric effects due
to topography on MODIS vegetation index data simulated from AVIRIS imagery over mountainous terrain, Can. J. Remote Sens. 18(4):283-291. Zibordi, G., and Voss, K. J. (1989), Geometrical and spectral distribution of sky radiance: comparison between simulations and field measurements, Remote Sens. Environ. 27: 343-358.
Topographic Effects in AVHRR NDVI Data D. W. Burgess,* P. Lewis,* and J.-P. A. L. Mulled A n investigation is made of the errors induced in the normalized difference vegetation index (NDVI) by topographic variation. The Advanced RAdiometric RAy Tracer (ARARAT) is used to simulate AVHRR imagery for various viewing and illumination conditions at 1.1-km and 50-m resolution. Topographic error is calculated as the difference between the variable terrain simulation and a flat plane simulation. Correlation between this error and various topographic factors such as shadowing and adjacement hill illumination are examined. The importance of sky occlusion as a topographic factor is highlighted. While topographic errors of 13.5% are found in the 50-m resolution simulations, these errors are reduced markedly in the 1.1-km dataset to approximately 3%.
INTRODUCTION Over recent years, the normalized difference vegetation index (NDVI), derived f~om Advanced Very High Resolution Radiometer (AVHRR) data, has become the de facto standard for monitoring vegetation growth on a global seale. It is now recognized in the literature that a number of factors lead to variations in the NDVI. These include atmospheric effects (Cihlar et al., 1994), soil brightness and canopy density (Huete et al., 1985; Borel and Gerstl, 1994), and sensor degradation (Gutman, 1991). It has also been recognized that the NDVI is not invariant under different viewing and illumination geometries and a number of empirical bulk corrections have been proposed to compensate for the surface anisotropy (Gutman, 1991; Goward et al., 1991; Taylor et al:, 1985).
*Manaaki-Whenua-Landcare Research, Wellington,New Zealand tDepartment of Geography,UniversityCollege London, United Kingdom *Department of Photogrammetryand Surveying,UniversityCollege London, United Kingdom Address correspondence to D. W. Burgess, Manaaki-WhenuaLandcare Research, Box 38-491, Wellington,New Zealand. Received 10 August 199.4;revised 21 June 1995. REMOTE SENS. ENVIRON. 54:223-232 (1995) ©Elsevier Science Inc., 1995 655 Avenueof the Amerieas, New York, NY 10010
A more complex effect on the NDVI arises due to topographic variation. Encompassed in this issue are the effects of shadow, adjacent hill illumination, sky occlusion, and slope orientation with respect to the bidirectional reflectance distribution function (BRDF) of the cover type (Schaaf et al., 1994). This article presents results of a series of simulation studies using the Advanced RAdiometric RAy Tracer (ARARAT) (Lewis and Muller, 1992), which simulates the radiation regime, for a given set of illumination and viewing conditions, using 3-D models of individual or fields of plants with associated lower geometric boundary conditions described by a digital elevation model (DEM) which is used to represent soil roughness or macroscopic terrain effects. Simulations over AVHRR wavebands and viewing and illumination geometries have been carried out both at the nominal 1.1-km resolution (nadir viewing) and at a higher resolution of 50 m (nadir viewing), the latter allowing topographic effects to be more easily identified.
BACKGROUND In general, successful NDVI studies are carried out in large, relatively fiat parts of the world such as the extensive cropping areas of the midwest of the United States, or areas of Africa. This type of terrain is ideal, given that the resolution of the data at nadir is relatively low at 1.1 km, and a large fiat test site with contiguous ground cover allows statistically significant measurements to be made. However, the consequence of this is that the issue of the effect of topography on the NDVI measure has largely been ignored in the literature. If there is a dependence on topography, then there are many terrestrial areas where it is not reasonable to ignore the effect.
Radlometric Topographic Effects The topographic effect is defined as the variation in radiance from inclined surfaces compared with radiance from a horizontal surface, as a function of the orientation 0034-4257 / 95 / $9.50 SSDI 0034-4257(95)00155-T
224 Burgess et al.
of the surface relative to the light source and sensor position (Holben and Justice, 1980). In other words, radiometric topographic effects can all ultimately be described as resulting from the "sun-sensor--target" geometry. However, to fully address the impact of topography, it is necessary to consider neighborhood effects alSO. If the earth's surface were an isotropic reflector, giving the same surface-leaving radiance regardless of viewing angle, then variations in the radiance received at the satellite would be simply a combination of atmospheric effects on the upwelling radiance and illumination angle effects. However, it has been shown that this is not the case (Smith et al., 1980; Gutman, 1991). A BRDF is needed to characterise the behavior of direct illumination on, and reflectance of, a given cover type, dependent on the incident and exitant angles of the radiance, as well as the radiometric and geometric properties of the land cover. In practice, this means that the solar zenith angle, view zenith angle, view azimuth, and terrain slope all affect the radiance measured for a particular ground cover. Therefore, the terrain slope at a given point is responsible for the point topographic effect. A less obvious effect of topography has been investigated by Teillet and Staenz (1992) with reference to AVHRR and simulated Moderate Resolution Imaging Spectroradiometer (MODIS) imagery. They studied the percentage change in NDVI measurements due to the different atmospheric path lengths resulting from terrain elevation, and showed that, with 23-km visibility, the change in AVHRR NDVI was 0.79% per 100-m elevation, for standard vegetation. They also showed that the change in NDVI was greater for less vegetated targets. Other topographic effects arise due to the context of the point within the surrounding terrain. These could be termed neighborhood topographic effects. The most obvious example is the result of direct illumination by the sun, causing shadows. The magnitude of this effect in raw satellite data is dependent on the solar elevation, decreasing as the Sun approaches zenith. Harder to quantify and model are the effects of diffuse illumination on sloping surfaces. The proportion of diffuse illumination is dependent on the turbidity of the atmosphere including the presence of cloud. Proy et al. (1989) identified two major interactions of diffuse illumination and topography. First, in steep terrain, pixels are often hidden from part of the sky hemisphere by adjacent terrain. This reduces the total diffuse irradiance reaching the point and hence affects its radiance. Second, terrain can be illuminated by reflected radiation from adjacent slopes, which has the effect of falsely increasing the reflected radiation from that point. Proy has shown that, while these diffuse irradiance effects may be negligible for directly illuminated pixels, they
may represent up to 33 % of the illumination contribution of shadowed pixels. Previous Work There have been several investigations made into the effects of topography on high resolution satellite data (e.g., Smith et al., 1980; Hall-K6nyves, 1987; Newton et al., 1991). Studies have also been made in controlled environments using hand-held radiometers (Holben and Justice, 1980; Holben and Justice, 1981; Pinty et al., 1990). In all cases, measurable effects were identified. While high resolution studies are invaluable for providing an insight into the mechanisms governing how topography affects the radiance sensed, for lower resolution data such as AVHRR, practitioners have tended to take a more statistical approach. For a number of years it was believed that the NDVI was invariant under different illumination and viewing conditions, and it has only been recently that researchers have started to consider that this may not be a valid assumption. Gutman (1991) has examined the frequency of individual view angles contributing to 10-day NDVI maximumvalue composite images over 1 month. Contrary to previous theories, he found that the NDVI composite was not biased towards nadir views, with their shorter atmospheric path lengths, but towards forward scattering views. This adds weight to the argument that the effects of surface anisotropy are carried through to the NDVI measures. He goes on to develop a simple empirical correction for AVHRR reflectance data based solely on view zenith angle. This had the effect of reducing the coefficient of variation of the NDVI by a factor of 2 from 0.12 to 0.06, where the mean corrected NDVI was 0.49. Cihlar et al. (1994) examined the BRDF effects in AVHRR data on NDVI compositing. They found them to be significant in four distinct cover types but noted that the anisotropy is reduced when data are atmospherically corrected, which it to be expected since the atmospheric correction process moderates the BRDF effects of the atmospheric phase function. METHOD
Gathering an adequate data set for investigations in this field is inherently very difficult, due to the fact that good sampling of the illumination and viewing geometry is required within the constraints of orbital characteristics, cloud cover, variation in sensor characteristics across the NOAA series, and atmospheric effects. This study has therefore taken advantage of the unique capabilities provided by a software simulator called the Advanced Radiometer Ray Tracer (ARARAT) for simulating AVHRR data. AVHRR Band 1 and Band 2 data have been simulated, with the resulting images yielding "at ground"
Topographic Effects in A VHRR NDVI Data 225
radiance measures. These radiances have then been converted to NDVI using the usual formula: (b2- bl)/ (b2 + bl). Although this step is more usually carried out on data that has been processed to reflectances, it has been shown that the relationship between radiancebased and reflectance-based NDVI measures is a simple one. It results from the multiplicative effect of the near infrared/visible ratio of' solar spectral irradiance, and can be approximated to a linear relationship over a NDVI range of 0.0-0.8 (Goward et at., 1991). While this does lead to a modest decrease in the range of the NDVI, 1 this was not considered significant for the generalized comparisons made in this study. The first series of simulations have been designed to look at the magnitude of variation possible due solely to varying Sun-sensor-target geometry. For this purpose, a flat plane DEM has been used, and the solar and view zenith angles have been varied. In this case the Ahmad-Deering BRDF model for prairie grass in June has been used (Ahmad and Deering, 1992). The second stimulation focuses on higher order effects of topography. A 12.5-m resolution DEM has been used with the simulation performed at 50-m resolution. Although considerably finer than actual AVHRR resolution, this detail has been necessary to identify the higher-order mechanisms of topographic effect. Because of the "bias" introduced by a BRDF, these simulations have first been run using a Lambertian reflectance function; then another simulation has been performed over the same area using the previously mentioned BRDF. Finally, a simulation has been performed at the AVHRR nominal resolution of 1.1 km over a DEM. An investigation of the correlation between this result and actual AVHRR data has been made. The variation in apparent topographic effect between this and the previous simulation has been examined. An Overview of ARARAT The first use of Monte-Carlo ray-tracing for satellite sensor data [Landstat Thematic Mapper (TM)] simulation was described by Muller and Dalton (1988). ARARAT is a second generation version which more correctly models the physics and can simulate a wider range of geometric primitives. ARARAT is described in detail in Lewis and Muller (1992). The aim here is only to give an overview of the functionality of the system. The basic principle is that of a Monte Carlo ray tracer (Cook et al., 1984; Kajiya, 1986; Glassner, 1989). Geometric models are rendered for a given set of illumination and viewing conditions. The sensor is characterized by a location and a pointing vector, for a given instanta-
l A radiance NDVI range of [0.0-0.8] corresponds to a reflectance NDVI range of [0.18-0.88J; thus the extent of the range is reduced from 0.8 NDVI units to 0.7 NDVI units.
Figure 1. ARARATsimulation of a millet field.
neous field of view. A modulation transfer function for the optics and finite aperture effects can be incorporated, but is not used in these simulations. The irradiance field is modeled by a gridded angular "map" of spectral sky radiance distribution and a (direct) solar spectral irradiance component. The original design focus of ARARAT was the simulation of ground based radiometer measurements. For this purpose, 3-D models of plants can be rendered, often with a DEM to represent the underlying soil structure. Figure 1 is an example of the kind of image that ARARAT produces at this scale. Because of the detailed physical basis to the simulation process, the software can be used at any spatial scale. Recent work by some of the authors, on the NASA First ISLSCP Field Experiment (FIFE) test site, employed ARARAT to simulate MODIS and AVHRR imagery. For satellite image simulations, account must be taken of the effects of atmospheric attenuation and scattering on both the directional irradiance function and the path from the ground to the sensor. Only the former are modeled in ARARAT at the present time; therefore, ARARAT simulates at-ground radiances. In order to simulate the data from different sensors, over
226 Burgess et al.
Imaging G e o m e t r y at Wellington, N Z ~"
~
'
~
0
Relative Azimuth Angle
AVHRR NOAA-9 L a L - 41 S L o n g . - 175
15 D e c e m b e r 2 January
270 s i ' 255 \
~
105
135
225
~
r
Solar Zenith A n g l e I
50.00 - 55.00
I 1
4s.00- 50.00 4o.00- 45.oo
1
3s.00- 4o.oo 30.00 - 35.00 25.00 - 30,00
J
Figure 2. X-SATVIEW (Barnsley et al., 1994) plot of the
summer imaging geometry of the NOAA-9 AVHRR satellite for Wellington, New Zealand.
the various spectral bands, ARARAT performs Monte Carlo sampling of a broadband sensor spectral response function (Lewis and Muller, 1992) as part of the rendering process. The simulator incorporates the Zibordi and Voss (1989) sky irradiance model which takes such inputs as ozone concentration, water vapor concentration, and aerosol scattering albedo. It also incorporates a number of bidirectional reflectance models including the Ahmad and Deering model (Ahmad and Deering, 1992) and the Pinty-Verstraete-Dickinson (PVD) model (Pinty et al., 1990). ARARAT uses reverse ray tracing, which is more computationally efficient than forward ray tracing for simulating satellite imagery (Muller and Dalton, 1988). The history of a ray is traced from the sensor to the surface, then out to the illumination source. Rays are allowed to interact with surrounding surfaces as they are reflected, until they finally escape the model to the sky. Direct illumination sampling is achieved by firing a ray directly at the sun and diffuse illumination sampling is done by firing sample rays over a hemisphere around the local surface normal according to the local surface reflectance function. Diffuse contributions from adjacent slopes are implicitly included in this algorithm. Not surprisingly, this process is computationally very expensive.
The Study Domain: Geometry A study was made of the viewing and illumination geometry of AVHRR data captured over New Zealand. Currently, afternoon passes are received daily from the NOAA-11 satellite, for the purpose of NDVI monitoring. Figure 2, which was generated using the package X-SATVIEW (Barnsley et al., 1994), shows the nominal geometry for an afternoon NOAA satellite, in midsum-
mer, for the data captured over central New Zealand. Because of the orbital characteristics of the satellite, this geometry does not vary greatly up and down the country. The relative azimuth angle is the angle between the .view and solar azimuth. From this figure, it is clear that the majority of images are very close to the principal plane, with relative azimuth angles close to 0 ° or 180 °. The solar zenith angle falls in the range of 25-55 ° , with the highest solar elevations (smallest solar zenith angles) occurring in the forward scatter direction. In midwinter, the data move out of the principal plane and are oriented along the 135 ° line. The solar zenith angle drops to a range of 65-85 ° . It should be noted that, while the viewing geometry varies significantly across a full AVHRR swath in a given scene, the variation over a small region of interest is small when compared to inter-scene variations. One restriction of the sky model used by ARARAT is that it models the atmosphere as parallel planes. This means that the simulator cannot be used for solar zenith angles beyond 70 ° , for which the model would not accurately reflect the true optical depth. For this reason the solar zenith angles investigated have been restricted to 25-55 ° - the lower limit being the summer limit for New Zealand and the upper limit chosen to be well within the acceptable simulation limits and to provide useful illumination levels. As indicated previously, the simulations are very time-consuming to run, so that further restrictions have been placed on the domain being investigated. Holben and Justice (1980) found that the greatest errors occurred in the principal plane. Since the summer imagery is oriented near the principal plane, and covers the solar zenith angle range already selected, only the principal plane situation is studied.
The Study Domain: Topography A site was chosen near the east coast of the North Island of New Zealand. It is hilly terrain, which has gradients ranging from 0 ° to 49 ° with an average of 20 °. The predominant ground cover is rye grass and clover, which is used for pastoral farming. In the 5 km x 5 km square site used for the high resolution simulation, the altitude ranges from 40 m to 500 m above sea level with a mean altitude of 195 m. Figure 3 shows the DEM used for this test area, which has a resolution of 12.5 m. The 22 km x 22 km area of the low resolution simulation (Fig. 9c) ranges in altitude from 21 m to 620 m above sea level with a mean elevation of 189 m. The average gradient is 16 ° within a range of 4-64 °. The digital elevation model used in this case was at 100 m resolution.
Simulating Imagery The first step in simulating the imagery was creating a sky model. In this case, a sky model appropriate to New
Topographic Effects in AVHRR NDVI Data 227
80.0
Band 2
¸%¸¸:¸ ~--~ 6 0 . 0
8 c
40.0
0--. • "'(3
@. /
-\'@. A. -.~. -.~--0 ~.~.Z~.~x_ -
.~ ~<
~..\-~
~.
.
sza=25 ®G---®-:E3 sza-~o ~" , ~- --~ i -4> sza=35 / u, / YJ- - - ~ ~ - ~ sza=40 / /~ ~ ~ sza--45 . / 7"< ~I~ ~ sza=SO
/~x- . & .
.~-A- - --A-- - -~
p<-
/"
~ .
.
.,~-~
/
, [,
]:>
./
®
¢n
Band 1 20.0
o'%.0
~;.o
o.o
+o.o
-100.0
View Zenith Angle (degrees)
Forwardscatter - Backscatter
Figure 3. 5 km x 5 km DEM of test area for simulation series 2 (12.5-m resolution).
Zealand atmospheric conditions was required. The local precipitable water vapor concentration was taken from tephigram data received at the meteorological station at Paraparaumu Airport near Wellington. The ozone concentration coefficient over New Zealand was obtained by interrogating the online database of Total Ozone Mapping Spectrometer (TOMS) data at the Geophysical Data Facility at Rutherford Appleton Laboratory, Oxford, United Kingdom. Appropriate coefficients of the/~ngstr6m formula were obtained from ~ngstr6m (1961) and others were set to average maritime values. The behavior o f the simulator can be customized in a number of ways. Because it contains an element of random sampling, results can be different from one simulation to the next. Therefore, it is vital to set the simulation parameters so that the results are convergent. The principal way to achieve this is by increasing the number of rays fired per pixel. After trials, it was found that 100,000 rays/pixel gave satisfactory results, with agreement in the output radiance for successive simulations to four significant figures. Another important parameter is the ray tree depth. This controls the number of intersections of the ray that are followed. As the signal is attenuated at each bounce, its contribution eventually becomes microscopic, if it does not terminate by escaping the terrain to the sky. However, it is imporatant to keep track of small contributions since their cumulative effect can be substantial. Experience has shown that, in general, a ray tree depth of 10 is adequate. This was verified during the simulations performed by monitoring the illumination contributions at various depths of the ray tree to ensure that the percent contributions at level 10 were small. These contributions ranged from 0.0005% down to < 0.000001%. Another customizable feature of ARARAT is the range of BRDF models that can be utilized. For the
Figure 4. Flat plane simulated radiances for June prairie grass for a range of solar zenith angles.
purposes of this study, the Ahmad-Deering model (Ahmad and Deering, 1992) was used with a cover type of prairie grass. This model was considered the best available to date and the inclusion of prairie grass parameters made it undoubtedly the best option for simulating a pasture type of cover. A Lambertian surface with a "green vegetation" reflectance was also used in a number of simulations. RESULTS Simulation Series 1--Flat Plane with a BRDF The simulations of AVHRR sensor response over a flat plane were carried out for solar zenith angles of 2555 ° in 5 ° steps, over a range of view zenith angles from 60 ° backscatter to 60 ° forwardscatter in steps of 15 °. The Ahmad-Deering (1992) June prairie grass BRDF model was used. In each case, the hotspot was also stimulated. The results are illustrated in Figure 4. It can be seen from this graph that both viewing and illumination angle have a large impact on "at ground radiance." Variations of up to 50% are apparent across different solar zenith angles. Figure 5 shows how this variation is propagated through to the NDVI measure, where the NDVI ranges from 0.47 to 0.68. Here we find that the backscatter hotspot becomes a minimum and the maximum NDVI values occur in the forwardscatter direction. This is consistent with the findings of Gutman (1991) in his study of the viewing geometry of pixels contributing to 10-day NDVI maximum value composites over the U.S. Great Plains. The variation in NDVI across the range of view zenith angles is approximately 34%, and the greatest variation across the solar zenith angles, at nadir, is approximately 11%. These figures have significant
228 Burgess et al.
0.65
0.60
I z
0.55
sza=26
0.50
sza=30 sza=35 sza=40 sza--45 sza=50 sza=55
Figure 6. NDVI image of simulation 2 with a
0.45
solar zenith angle of 25 ° . Area corresponds to 1 km x 5 km at right edge of DEM. 0"4000.0
i 50.0
0.0
-50.0
-100.0
View Zenith Angle (degrees) Forwardscatter- Backscatter
N D V I % _ error(x,y) Figure 5. Flat plane NDVI variation for June prairie grass for a range of solar zenith angles.
implications when considering the role of the NDVI measure in land cover monitoring. Even restricting the viewing geometry to + 30 °, NDVI variation is of the order of 24%. Simulation Series 2 - - A c t u a l Terrain with a
Lambertian Surface The second set of simulations was performed over actual terrain using the 5 km x 5 km DEM described previously, with a Lambertian surface reflectance model. The processing overhead with this series was far greater than for the first for two reasons. First, a 100 x 20 pixel image strip down the DEM was simulated rather than a single point, as in the flat plane case. Second, the processing of a single point is far more complex in varying terrain, since many interactions occur during ray tracing with adjacent hill facets. Because of this greater overhead, it was only possible to simulate three solar zenith angles (25 °, 40 °, and 55 °) for a nadir viewing position. "AVHRR" pixels (50 m) were generated over the 12.5-m DEM to produce a 100 × 20 image corresponding to a strip down the right-hand side of the test area shown in Figure 3. Corresponding flat plane simulations were performed for each terrain simulation to provide an NDVI measure with no terrain effect. These were essentially the same as the first series except that they used the Lambertian surface rather than a BBDF. The NDVI was calculated across each terrain simulation and compared to the flat plane NDVI. Figure 6 shows the resulting NDVI image for the case of a solar zenith angle of 25 °. Errors were expressed as a percentage of the flat plane values as follows:
(INDVIt~,,ai,(x,y) - NDVI~t(x,y)I / x 100, =\ NDVI~a,(x,y) /
(1)
for each pixel (x,y) in the scene. The mean scene N D V I % _ e r r o r was then calculated for each solar zenith angle. Table 1 shows that these errors were consistently in the region of 13.5% regardless of solar zenith angle. It must be appreciated that these figures are not a representation of total NDVI error due to viewing and illumination geometry and terrain effects but solely the latter. The former effect is eliminated by normalizing the error to the fiat plane generated for the same viewing conditions. These results were compared to a number of terrain factors in an attempt to establish the cause of this error. As part of its output, ARARAT generates a list of ray tree depth contributions. This means that, for each level in the tree of ray interactions, information is gathered on the diffuse and direct illumination contributions incident on the terrain surface at that point. From the first level of this tree we get the percent direct and percent diffuse illumination directly coming to the pixel. Therefore, the percent direct illumination is a measure of how much the pixel is in shadow. Table 1 shows that for Band 2 there was a strong negative correlation between percent direct illumination and percent NDVI error, which implies that shadowed pixels gave greater NDVI error than directly illuminated ones. From all the levels below level 1 in the ray tree we get the illumination information for the adjacent hill points that reflect light onto the target pixel. By adding up all these contributions, it is possible to obtain a measure of the total percent adjacent hill illumination. The extent of adjacent hill illumination is largely due to the complexity of the terrain. Steep sided valleys will cause light to be reflected from one side to the other. Table 1 shows that percent adjacent hill immuniation is highly correlated with NDVI error.
Topographic Effects in A VHRR NDVI Data 229
Table i. Results from Second Set of Simulations Correlation of % NDVI Error with: % Illumination Solar Zenith Angle
Mean % NDVI Error
% Sky Visible
25
13.5
- 0.99
0.99
0.99
- 0.62
- 0.98
40
13.4
- 0.98
0.98
0.98
- 0.48
- 0.92
55
13.5
- 0.91
0.95
0.93
- 0.61
- 0.82
Of greatest interest in determining significant terrain factors was the correlation of the percent NDVI error with a percent sky visible image. This image was created also using ARARAT but with no sun, a uniform sky hemisphere and a perfectly white diffuse reflectance function. Under these conditions, the resultant image is equivalent to the percent of sky visible at each point. Figure 7 shows this image. The minimum percent sky visible value is 19.8%; the maximum is 99.6% and the mean is 68.2%. In Table 1, it can be seen that this image was consistently highly correlated with the percent NDVI error image. Obviously this factor is related to diffuse illumination, however, the beauty of it, as a potential tool for normalizing NDVI data for terrain effects, is that it is viewing and illumination geometry independent. Neither the Band 2 direct illumination measure nor the percent adjacent hill illumination measure have this advantage. Simulation 3 - A c t u a l Terrain Using a BRDF Having identified, first, the variation in NDVI caused by viewing and illumination geometry over a fiat plane using a BRDF, and second, the effect of actual terrain on the NDVI when a Lambertian surface is assumed, the natural progression is to look at what happens when a BRDF is used in an actual terrain simulation. Utilizing a BRDF leads to further processing overheads, which Figure 7. Percent sky visible image over area corresponding to the DEM in Figure 3.
Adjacent Hill Band 1 Band 2
Direct Band 1 Band 2
unfortunately made multiple simulations in this series impractical. The one completed simulation of the 100 x 20 pixel strip took approximately i week to complete, running in parallel on 18-24 Sun and Silicon Graphics Workstations. This simulation used the same BRDF as in the first fiat plane series and was done at a solar zenith angle of 55 ° and a view zenith angle of - 5 0 °, that is, near the hot spot. The resulting mean percent NDVI error was 16.1%, which was greater than in the second simulation (13.5%) but not by as much as might be expected in light of the results of the first fiat plane results. However, closer examination of the nature of BRDFs in general, and in particular the one utilized in this simulation, reveals why the difference in error is not great. As the slope orientation in hilly terrain varies
Figure 8. Prairie grass BRDF at 662 nm showing terrain variation line.
Praire Grass BRDF X= 6 6 2 p r n
0.065 0.06 0.055 0.05 ~ o114 3~ 0.0.003~
5
0.025~O "40 5
•
: sza
- v z a = 5°
230
B u r g e s s e t al.
(a)
Simulation 4--Actual Terrain at l.l-km Resolution Having established that various topographic effects appear to have an impact on the measured NDVI in 50 m resolution data, the next step is to see if this variation is still present at the scale of AVHRR data. An area of 20 x 20 pixels was simulated (22 km x 22 km) over the test site as indicated previously, for a solar zenith angle of 55 ° with a nadir viewing position and a Lambertian green surface. Band 2 of the 1.1-km resolution resulting image is shown in Figure 9. A coarser resolution DEM of 100 m was used in this simulation, which was considered appropriate given the desired simulation output resolution of 1.1 km and the processing overheads of using the finer, 12.5-m resolution DEM over such a large
(b)
area.
(c) Comparison of actual and simulated data: a) Simulation 4, Band 2; b) AVHRR data, Band 2; c) DEM of area for a) and b). F i g u r e 9.
over a small geographical area, both the effective viewing and illumination angle change with respect to the surface normal, but the angle between the sun and viewing direction does not. Figure 8 shows the AhmadDeering model for Prairie grass at 662 nm (Ahmad and Deering, 1992). The line of triangles indicates the case where the view zenith angle plus the solar zenith angle is equal to 5 ° . This is the case in the third simulation with a solar zenith angle of 55 ° and a view zenith angle of - 5 0 °. It can be seen that as the terrain varies, the BRDF is sampled along a line: vza = - sza + a,
(2)
where v z a is the view zenith angle, s z a is the solar zenith angle, and a is a constant representing the angle between the v z a and s z a . In general, the variation in the BRDF along this line is less than the variation over view zenith angles for a given solar zenith angle; therefore, variations in overall viewing geometry across an AVHRR scene and between distinct scenes will have a greater effect on the measured radiance than variation caused by slope orientation.
Again, flat plane values were also simulated and the resulting NDVI values compared. In this case the average percent difference between the terrain and flat plane NDVI values was approximately 3%. This compares with the 13.5% variation of simulation 2 and the 16.1% variation of simulation 3. In light of the lower errors in the coarser resolution simulation, it was speculated that perhaps this was due to the coarser resolution DEM, and not a direct result of larger pixel size. A smaller 10 km x 10 km subscene was resimulated using the high resolution 12.5-m DEM. An increase of approximately 1% in the mean percent NDVI error was found. This suggests that there is only a slight sensitivity to the DEM resolution used. Therefore, it would appear that, at the coarser resolution, the effects of topography on the resulting radiance are moderated since, within one pixel area of 1.1 km 2, a larger sample of slopes and illumination conditions are represented. An attempt was made to correlate the simulated data with actual data received over the same area. Figure 9b shows actual AVHRR data corresponding to the simulation in Figure 9a. Figure 9c shows the DEM of the same area. Similarities can be seen in the flatter area in the bottom left quandrant; however, the overriding impact of the dark area of native forest on the southeastern face of the ridge in the top left quadrant makes it hard to discern any correspondence in this area. Efforts to correlate these datasets were therefore largely unsuccessful. It is clear that the nominal contiguous pasture cover of the test site actually contained a significant amount of variation in the form of shelter belts and small stands of native bush. Where these areas were clearly indicated on a 1:50,000 scale map, for example, an area of bush mentioned previously, they were masked out, along with a few pixels of ocean in the lower right-hand corner of the scene. However, correlation over the unmasked areas was still not possible. This was attributed to small scale variation in the cover characteristics, which were not modeled by the simulator and difficulties in accurate subpixel registration of the datasets.
Topographic Effects in A VHRR NDVI Data 231
CONCLUSIONS Valuable insights have been gained into the mechanisms of topographic effects on the NDVI measure. The fiat plane series of simulations highlighted the significance of the variation in measurement due solely to viewing and illumination geometry. It also showed that the variation in viewing geometry had a greater effect than the sun position. The second series of simulations yielded a topographic error of 13.5 % in the NDVI values irrespective of illumination angle for 50-m resolution pixels. Shadowing, adjacent hill illumination and skT occlusion were well correlated to this error, suggesting that these mechanisms contribute to the variation. The impact of the BRDF in complex terrain was briefly investigated. It was found that, while surface anisotropy does add to the variability in the measured radiance, the constraints on the sampling of the BRDF imposed by the geometry of the system means that the variation is generally not as great as that caused by view angle variations for a given illumination angle. The fourth simulation showed that the previously identified NDVI errors were significantly reduced in 1.1-kin resolution data (3%). When this error is compared to the overall variation in NDVI due to Sunsensor geometry demonstrated in the first series of simulations (34% across all view zenith angles and a maximum of 11% across solar zenith angles), it is clear that the "bulk effect" of viewing and illumination geometry, irrespective of terrain orientation, is the dominant source of error in 1.1-kin resolution data. We suggest that, f~r the operational correction of AVHRR NDVI data, it is reasonable to ignore high-order topographic effects such as sky occlusion and adjacent hill illumination; however, bulk corrections should be made for viewing and illumination geometry variations across the full AVHRR swath in order to accurately compare NDVI measurements in different parts of the scene and in other scenes. Funds for this study were provided by The New Zealand Foundation for Research, Science and Technology under Contract C09208. D. Burgess wishes to acknowledge the support of the British Council, with whose help a significant portion of this work was undertaken in the United Kingdom, and Professor Harley and the Department of Photogrammetry and Surveying at University College London for support during the 6-month sabbatical. REFERENCES
Ahmad, S. P., and Deering, D. W. (1992), A simple analytical function for bidirectional reflectance, J. Geophys. Res. 97(D17):18,867-18,886. AngstrOm, A. (1961), Techniques of determining the turbidity of the atmosphere, Tellus 13:214-223.
Barnsley, M. J., Strahler, A. H., Morris, K. P., and Muller J.-P. (1994), Sampling the surface bidirectional reflectance distribution function (BRDF): 1. Evaluation of current and future satellite sensors, Remote Sens. Rev. 8:271-311. Borel, C. C., and Gerstl., S. A. W. (1994), Nonlinear spectral mixing models for vegetative and soil surfaces, Remote Sens. Environ. 47:403-416. Cihlar, J., Manak, D., and Voisin, N. (1994), AVHRR bidirectional reflectance effects and compositing, Remote Sens. Environ. 48:77-88. Cook, R., Porter, T., and Carpenter, L. (1984), Distributed ray tracing, Comput. Graphics (Proc. SIGGRAPH) 18(3): 137-145. Glassner, A. (1989), An Introduction to Ray Tracing, Academic New York, 1989. Goward, S. N., Markham, B., Dye, D. G., Dulaney, W., and Yang, J. (1991), Normalized difference vegetation index measurements from the Advanced Very High Resolution Radiometer, Remote Sens. Environ. 35:257-277. Gutman, G. G. (1991), Vegetation indices from AVHRR: an update and future prospects, Remote Sens. Environ. 35: 121-136. Hall-K6nyves, K. (1987), The topographic effect on LANDSAT data in gently undulating terrain in Southern Sweden, Int. J. Remote Sens. 8(2):157-168. Holben, B. N., and Justice, C. O. (1980), The topographic effect on spectral response from nadir-pointing sensors, Photogramm. Eng. Remote Sens. 46(9):1191-1200. Holben, B., and Justice, C. (1981), An examination of spectral band ratioing to reduce the topographic effect on remotely sensed data, Int. J. Remote Sens. 2(2):115-133. Huete, A. R., Jackson, R. D., and Post, D. F. (1985), Spectral response of a plant canopy with different soil backgrounds, Remote Sens. Environ. 17:37-53. Kajiya, J. T. (1986), The rendering equation, Comput. Graphics 20(4):143-150. Lewis, P., and Muller, J.-P. (1992), The Advanced Radiometric Ray Tracer: ARARAT for plant canopy reflectance simulation, in Proceedings of the XVI~ Congress of the ISPRS, Washington, DC, Commission 7, Vol. XXIX, pp. 26-34. Muller, J.-P., and Dalton, M. (1988), Application of ray-tracing to satellite image understanding, in Proceedings oflGARSS "88 Symposium, Edinburgh, Scotland, 13-16 September, pp. 1187-1190. Newton, A., Muller, J.-P., and Pearson, J. (1991), SPOT DEM shading for Landsat-TM topographic correction, in Proceedings of IGARSS "91, Espoo, Finland, 3-6, June Vol. 2, pp. 655-659. Pinty, B., Verstraete, M. M., and Dickinson, R. E. (1990), A physical model of the bidirectional reflectance of vegetation canopies. 2. Inversion and validation, J. Geophys. Res. 95(D8):11,767-11,775. Proy, C., Tanrd, D., and Deschamps, P. Y. (1989), Evaluation of topographic effects in remotely sensed data, Remote Sens. Environ. 30:21-32. Schaaf, C. B., Li, X., and Strahler, A. H. (1994), Topographic effects on bidirectional and hemispherical reflectance calculated with a geometric-optical model, IEEE Trans. Geosci. Remote Sens. 32(6):1186-1193. Smith, J. A., Lin, T. L., and Ranson, K. J. (1980), The Lam-
232 Burgess et al.
bertian assumption and Landsat data, Photogramm. Eng. Remote Sens. 46(9):1183-1189. Taylor, B. F., Dini, P. W., and Kidson, J. W. (1985), Determination of seasonal and interannual variation in New Zealand pasture growth from NOAA-7 data, Remote Sens. Environ. 18:177-192. Teillet, P. M., and Staenz, K. (1992), Atmospheric effects due
to topography on MODIS vegetation index data simulated from AVIRIS imagery over mountainous terrain, Can. J. Remote Sens. 18(4):283-291. Zibordi, G., and Voss, K. J. (1989), Geometrical and spectral distribution of sky radiance: comparison between simulations and field measurements, Remote Sens. Environ. 27: 343-358.