Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Radiometric correction of satellite imagery for topographic and atmospheric effects Gray, Malcolm Haig 1986

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


831-UBC_1986_A6_7 G72.pdf [ 7.19MB ]
JSON: 831-1.0096705.json
JSON-LD: 831-1.0096705-ld.json
RDF/XML (Pretty): 831-1.0096705-rdf.xml
RDF/JSON: 831-1.0096705-rdf.json
Turtle: 831-1.0096705-turtle.txt
N-Triples: 831-1.0096705-rdf-ntriples.txt
Original Record: 831-1.0096705-source.json
Full Text

Full Text

RADIOMETRIC CORRECTION OF SATELLITE IMAGERY FOR TOPOGRAPHIC AND ATMOSPHERIC EFFECTS by MALCOLM HAIG GRAY B.Sc, University of British Columbia, 1978 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF T H E REQUIREMENTS FOR T H E D E G R E E OF MASTER OF SCIENCE in T H E FACULTY OF G R A D U A T E STUDIES (Faculty of Forestry /Remote Sensing) We accept this thesis as conforming to the reauired standard T H E UNIVERSITY OF BRITISH COLUMBIA August 18, 1986 ©Malcolm Haig Gray, 1986 I n p r e s e n t i n g t h i s t h e s i s i n p a r t i a l f u l f i l m e n t o f t h e r e q u i r e m e n t s f o r an advanced degree a t t h e U n i v e r s i t y o f B r i t i s h C o l u m b i a , I agree t h a t t h e L i b r a r y s h a l l make i t f r e e l y a v a i l a b l e f o r r e f e r e n c e and s t u d y . I f u r t h e r a gree t h a t p e r m i s s i o n f o r e x t e n s i v e c o p y i n g o f t h i s t h e s i s f o r s c h o l a r l y p u r p o s e s may be g r a n t e d by t h e head o f my department o r by h i s o r h e r r e p r e s e n t a t i v e s . I t i s u n d e r s t o o d t h a t c o p y i n g o r p u b l i c a t i o n o f t h i s t h e s i s f o r f i n a n c i a l g a i n s h a l l n o t be a l l o w e d w i t h o u t my w r i t t e n p e r m i s s i o n . Department o f The U n i v e r s i t y o f B r i t i s h C o l u m b i a 1956 Main Mall V a n c o u v e r , Canada V6T 1Y3 Date ii A B S T R A C T The radiometry of satellite imagery is influenced by ground cover, local topography, and atmosphere. In order to increase the accu-racy of ground cover identification from satellite imagery, effects due to topography and atmosphere must be removed. These effects can be estimated by modeling the image-formation process. For this thesis an image-formation model is developed and tested on Landsat MSS data over a mountainous region. Solar illumina-tion angle, atmosphere depth, and sky illumination are calculated with the help of a digital elevation model. A digital forest cover map is used to select a target forest type for which model parameters are estimated using regression analy-sis. Results of this analysis indicate that solar illumination angle has the largest effect on target pixel irradiance followed by at-mosphere depth. Sky illumination as calculated, was significantly correlated with target pixel irradiance but in a negative sense. This correlation suggests that inter reflection (also called mutual illumination) from adjacent terrain may be a small but significant source of illumination. The estimated model parameters are used to correct the im-agery for topographic and atmospheric effects. Visual assessment of the corrected imagery indicates that many but not all of the to-pographic effects have been reduced. Comparisons between com-puter classified imagery and the forest cover map show an improve-ment in correctly classified pixels from 54% for the original image to 72% for the corrected image. i i i C O N T E N T S Abstract iii List of Tables vi List of Figures vii Acknowledgements viii CHAPTER 1 1 1.1 Introduction 1 CHAPTER 2 4 2.1 Thesis Statement 4 2.2 Research Objective 4 2.3 Methodology 5 2.4 Previous Work 6 CHAPTER 3 8 3.1 Image Formation Model 8 3.1.1 Geometry 9 3.1.2 Atmosphere 10 3.1.3 Illumination 11 3.1.4 Surface Photometric Properties 12 3.1.5 Topography 13 3.2 Image Formation Equation 15 iv CHAPTER 4 17 4.1 Study Site 17 4.2 Data Set 18 4.2.1 Landsat Imagery 18 4.2.2 Digital Elevation Model 18 4.2.3 Digital Forest Cover Maps 19 4.3 Model Parameter Estimation Method 19 4.3.1 Target Selection 19 4.4 Regression Analysis 20 4.5 Correction Parameters 23 CHAPTER 5 25 5.1 Assessment of Correction Accuracy 25 5.1.1 Atmospheric Parameters 25 5.1.2 Visual Comparison 26 5.1.3 Classification Results 27 CHAPTER 6 29 6.1 Summary 29 6.2 Conclusion 30 References 32 Appendix 35 Tables 37 Figures 42 v L I S T O F T A B L E S 1 Solar Irradiance at the Top 37 of the Atmosphere 2 Regression Results for Band 5 38 3 Confusion Matrices for 39 Classification Results 4 Point Source Coordinates for the 40 Asymmetric Sky Model 5 Point Source Coordinates for the 41 Uniform Tesselation of a Hemisphere V I L I S T O F F I G U R E S 1 Variation of Optical Depth and Upward and 42 Downward Transmission with Elevation 2 Components of Illumination 43 3 Definition of Incident Angle, Emergent 44 Angle and Phase Angle 4 Reflectance Maps 44,45 5 Target Radiance vs Elevation 46 6 Original Image 47 7 Corrected Image 48 8 Classification Results 49 9 Asymmetric Sky Illumination Model 50 10 Uniform Sky Illumination Model 51 vii A C K N O W L E D G E M E N T S I wish to express my sincere gratitude to my supervisor, Dr. Robert J. Woodham for his support and insight which have guided my work and for his insistance on clear thinking and concise writing. Grateful acknowledgement goes to Marc Majka whose program-ming assistance was invaluable to me. I also want to thank the people of the B.C. Ministry of Forests, Inventory and Planning Branch for providing the digital forest cover maps. My apprecia-tion is extended to fellow graduate students at the Laboratory for Computational Vision for their software tools made freely avail-able and for their helpful criticism. Finally I would like to thank Nedenia Krajci for drafting the figures and my wife Margot for assistance typing the manuscript. The funding for this thesis was provided by the Donald S. McPhee Fellowship Fund (1983-84), and by the MacMillan Bloedal Fellowship in Forest Mensuration (1984-85). viii Chapter 1 1.1 I N T R O D U C T I O N Earth observing satellites have given man a new view of his world. Images gathered from space are being used to map geographic features, inventory resources, and monitor change. The use of this information source is expanding for many reasons including: bet-ter images, cheaper display devices, and a growing body of knowledge on the subject. As man's capacity to exploit the earth has grown so must his knowledge and understanding. Satellite derived information can help provide that knowledge. Given its global scale, comprehensive nature, and frequent update, prehaps it can also promote understanding of the earth. Much of the world's forest resources are found in mountainous areas. Monitoring these resources using satellite remote sensing is difficult. Features due to topography and atmosphere frequently dominate the imagery, making identification of ground cover problematic. The motivation for this thesis is to increase the usefulness of satellite imagery for forest monitoring by removing topographic and atmospheric effects. The goal of remote sensing is to recover information from the imaged scene. For the research undertaken for this thesis the desired information is a surface reflectance factor, an intrinsic property of the surface independent of the imaging situation. Reliable recovery of reflectance factors demands an accurate model of the imaging process. The mathematical formulation of the model is embodied in the image-formation equation. It is necessary to invert the equation in order to express the reflectance factors as a function of image irradiance and other parameters and thereby compute the desired surface description. 1 Major aspects of any image-formation model are; the reflectance properties and orien-tation of the surface, the effects due to the atmosphere, and the geometric and spectral characteristics of the sensor and source(s) of illumination. Modeling the reflectance properties of the surface and atmospheric effects is the heart of the problem in satellite imaging. A diffuse or Lambertian reflectance model is a good first approximation of the re-flectance characteristics of the solid part of the earth's surface (Little [1980], Horn &; Bachman[1978]). For this reason and also because of its simple mathematical form, this model has been widely used. A more general relationship is given by the photometric function proposed by Minneart, which has been used to determine the surface charac-teristics of planetary bodies (Minneart [1941]). Other relatively simple, semi-empirical models have been investigated in attempts to characterize the reflectance properties of complex surfaces such as forest canopies (Kimes et al. [1980]). Atmospheric effects are extremely variable. Under clear sky conditions, a simple atmospheric model is used to predict the attenuation factors, path radiance, and distri-bution of the scattered light responsible for skylight. A difficulty arises in the fact that the atmosphere and ground cover interact, so that skylight is not only a function of the sun and atmosphere but also depends upon the reflectance properties and topography of the surrounding terrain. For this thesis an image-formation model is developed and tested on Landsat MSS data over a mountainous region. Topographic and atmospheric attributes of interest are computed from a digital elevation model (DEM). A digital forest cover map is used to select target areas for which model parameters are estimated using regression analysis. These parameters are used to correct the imagery for topographic and atmospheric ef-fects. An assessment is made of the correction accuracy both visually and by comparisons between computer classified imagery and the forest cover map. Chapter 2 describes the methodology employed in constructing and testing the image-formation model. Previous work in estimating topographic and atmospheric effects in satellite imagery is reviewed. Chapter 3 examines the image-formation equation in detail. Chapter 4 describes the data set, and develops the method for estimating model 2 parameters. Analysis of the results suggest modification of the model. The model is changed and the process repeated. Parameters estimated for the final chosen model are used to correct the image. Chapter 5 assesses the accuracy of the radiometric correction. The original image is visually compared to the corrected image. Ground cover classification results of both images are compared. Chapter 6 summarizes the results and discusses their relevance to current remote sensing issues. 3 Chapter 2 2.1 T H E S I S S T A T E M E N T Much of the world's forest resources are found in mountainous areas. Monitoring these resources using satellite remote sensing is often difficult. Features due to topography and atmosphere frequently dominate the imagery making ground cover identification prob-lematic. These effects can be quantitatively estimated by modeling the image-formation process. Satellite imagery can be corrected for these unwanted factors in order to increase its usefulness for forest monitoring purposes. 2.2 R E S E A R C H O B J E C T I V E Radiometric correction of Landsat data has been primarily motivated by the need to directly compare data obtained from different sites and at different times. In areas of rugged terrain, radiometric correction is required even if the measurements in one data set are to be directly compared only to those from other locations in the same data set. The goal of the research undertaken for this thesis is to improve simple (e.g. forest versus nonforest) ground cover classification from satellite imagery over mountainous terrain by correcting for topographic and atmospheric effects. It is preferable for wide application that the correction scheme not require independent measurements of atmo-spheric parameters for the time of satellite overpass. 4 2.3 M E T H O D O L O G Y Before constructing and testing a model of the image-formation process one must choose a conceptual framework for that model. For this research the primary choice was between a physical model or a statistical model. This choice is important to the approach used in setting up the model and for the interpretation of the results of fitting the model to the data. Many of the procedures used in testing the model will be the same regardless of the choice made. Physical models are based on a theoretical understanding of the physics of the situ-ation augmented by results from empirical studies. All factors believed to be important are explicitly included. Generally, physical plausibility is used as an important criterion in judging the results of fitting the model to the data. Statistical models differ from physical models in the selection of factors to be consid-ered and the criteria used to judge the results of fitting the model to the data. Usually such models employ additional variables compared to the corresponding physical models. These additional variables are often transformations and combinations of the variables used for the physical model. Statistical measures are the primary means of judging the correctness of the model with physical plausibility being secondary. Considering the work of previous researchers in this field, and considering the data set available for analysis, the choice was made to base the model on physical principles. The solar illumination angle, depth of atmosphere, and amount of sky visible was calculated for each pixel from the digital elevation model (DEM). This was used to explicitly model direct solar illumination, atmospheric transmission, path radiance, and sky illumination. Actual measurements of atmospheric parameters at the time of satellite overpass are unavailable. For this reason as well as the difficulty in obtaining these measurements routinely, these parameters were estimated from the model rather than being supplied to the model as auxiliary input. The forest cover map was used to choose and edit target areas. As the motivation for this work was to increase forest classification accuracy it was felt that using a forest cover type to estimate model parameters would be consistent with this goal. When selecting target forest types consideration was given to satisfying model assumptions. In particular 5 forest stands with closed canopies are felt to better approximate a Lambertian surface than stands with open canopies. Preferably for model verification, a target forest type should be widely distributed over a range of elevations, orientations and topographic positions. The image-formation model is stated as an equation expressing target pixel radiome-try as a function of direct solar illumination angle, depth of atmosphere and sky illumi-nation. The equation is linearized and multiple linear regression analysis is used estimate model parameters. The coefficient of determination (R2), standard error, and residuals were examined to judge the suitability of the various models investigated. Physical plau-sibility was a primary requirement in choosing the model to be used for the correction process. The parameters estimated for the chosen model and target area are used to correct the image by inverting the image formation equation. Visual and automatic comparisons between the original image, the corrected image, and the forest cover map are used to gauge the success or failure of the correction process. 2.4 P R E V I O U S W O R K Lambertian (diffuse) reflectance models have been widely tested. Smith et al. [1980]) investigated the Lambertian assumption for Landsat MSS data for Ponderosa pine (Pi-nus ponderosa Laws.) forest in Colorado. Their model did not explicitly include sky irradiance and path radiance. They concluded that their model was most valid for tar-gets with slopes of less than 25 degrees and solar illumination angles of less than 45 degrees. This is not surprising because sky irradiance becomes an important source of illumination (relative to direct solar illumination) at high illumination angles. Justice et al. [1980] compared four techniques for reducing the topographic effect in Landsat MSS data of the forested mountains of Pennsylvania. Their models did not explicitly include sky irradiance or path radiance. Results showed that their application of the Minneart function substantially reduced the topographic effect whereas the Lambertian model slightly increased the topographic effect. Teillet et al. [1982] formulated Lambertian and non Lambertian corrections taking 6 into account atmospheric effects as well as topographic effects, for Landsat MSS data and airborne MSS data of two areas in the Coast Mountains of British Columbia. Mod-erate but statistically significant correlations were found between MSS data and effective illumination angle. Similar correlations were found with angle of reflection to the sensor. These correlations were stronger for MSS bands 6 and 7 than for band 4 and 5. Ground cover type classification accuracies for corrected imagery were not significantly improved except for one specific model that included path radiance. Ahern et al. [1977] have shown it is possible to extract both path radiance and trans-mission information from Landsat data using clear water bodies as dark targets of known low albedo. This method was applied to areas of the Canadian shield exhibiting low re-lief. To determine path radiance and transmission factors as a function of elevation using this method would require clear water lakes at many elevations. Sjoberg and Horn[l983] used published estimates of path radiance, optical depth and sky irradiance to correct a Landsat MSS image of mountainous terrain. The resulting albedo map was an effective aid in adjusting the model parameters to obtain an ac-ceptable representation. Lee [1983] used the same model as Sjoberg but attempted to estimate the atmospheric parameters by examining the image intensities as a function of elevation and by examining image intensities across shadow boundaries. The work undertaken for this thesis is based on data from the same study area used by Woodham and Lee [1985]. The present work makes use of additional information unavailable to Woodham and Lee, in the form of a digital forest cover map. This allows much flexibility in choosing and analysing target areas. Also the forest cover information is used to validate the corrected image in addition to subjective visual inspection. The skylight models used here are elaborated to include the effects of obscuring terrain and a nonuniform skylight distribution. 7 Chapter 3 3.1 I M A G E F O R M A T I O N M O D E L The image formation process must be understood before unwanted factors can be cor-rected for. By explicitly modeling the process we can test and improve our understanding. There are five major factors to consider in an imaging situation: 1. the geometry of the source, viewer, and scene 2. the medium through which energy is transmitted, the atmosphere 3. the spatial and spectral distribution of illumination 4. the intrinsic radiometric properties of the viewed surface 5. the topography of the viewed surface, it's shape In remote sensing using Landsat imagery these factors are determined as follows: 1. The relative geometry of the earth, sun, and satellite is determined by the date and time of image acquisition. Landsat views the world from approximately overhead everywhere at the same local solar time, with the position of the sun in the sky varying according to the season. 2. The atmosphere scatters and absorbs energy passing from the sun to the surface and from the surface to the satellite. The scattering of energy by the atmosphere is responsible for sky illumination and path radiance. Path radiance is energy scattered towards the sensor by the column of air viewed by the sensor. 8 3. Direct solar illumination is the primary source of illumination and also contributes diffuse sky illumination by it's interaction with the atmosphere. Reflection from surrounding terrain may also be a significant source of indirect illumination and it also contributes to skylight and path radiance. 4. The ground cover and it's condition and structure determine the photometric prop-erties of the surface. 5. The topography is the shape of the surface and can be described in various ways such as a contour map or for this work by a digital elevation model that assigns elevation values on a regular horizontal grid. 8 . 1 . 1 G E O M E T R Y Landsats 1, 2 and 3 have a nominal altitude of 900 km with a maximum scan angle (deviation from the vertical) of 5.78°m (Landsat Data Users Handbook [1979]). Thus relief displacement is small, requiring an elevation change of 507 (at maximum scan angle) for a 1 pixel horizontal displacement (Wong [1984]). Because of the small viewing angle, Landsats 1,2 and 3 scan a relatively small portion of the earth's surface, allowing us to ignore the curvature of the earth. The imaging geometry can be approximated by orthographic projection. With an orthographic projection the mapping between object space and image space becomes trivial. A point (x,y,z) in the object space can be considered to map into point (u,v) in the image plane where x = u and y = v. Therefore coordinates (x,y) and (u,v) can be used interchangeably. As a result of their sun-synchronous orbit, Landsat satellites always pass over a given part of the earth at approximately the same local solar time (9:52 AM for the image used here). Consequently the major geometric variable is the position of the sun. Given a date, time, and position on the earth the elevation and azimuth coordinates of the sun can be determined (Horn [1978]). 9 3.1.2 A T M O S P H E R E For the model developed here the atmosphere is assumed to be an optically thin, horizon-tally homogeneous layer. Optically thin considers only primary scattering in describing the scattering mechanism. This is often a reasonable assumption but should be tested by an inspection of the imagery. Localized haze and thin linear clouds have been detected on seemingly clear imagery over the study area for a different date [Sept. 25 1974] than that used for this research. Optical thickness is the attribute which mostly determines the atmospheric trans-mission factors, path radiance and sky irradiance (skylight). The optical thickness r of an atmospheric mass measures the total extinction experienced by a light beam passing through it. Optical thickness depends upon the density and size distribution of particles in the atmosphere and on their scattering and absorption properties. For the visible and near infrared portion of the spectrum, three classes of particles are relevant: molecules, aerosols and ozone. Molecules of air are small with respect to the wavelength of light and contribute to Rayleigh scattering. Areosols are particles large with respect to the wavelength of light and contribute to Mie scattering. Ozone absorbs radiation. Based on published standard (i.e. clear sky) atmospheres (Valley[1965]), the Rayleigh, the Mie, and the total optical thickness varies almost exponentially with elevation. Let r(z) = r0e-z/H* where r 0 = r(0) is the optical thickness at sea level and parameter Hr is the scale height. For a horizontally homogeneous atmosphere the upward and downward transmissions become Tu{z) = e~rW Td{z) = e-'W'03^ Figure 1 illustrates the variation of optical thickness and upward and downward trans-missions with elevation , for green light with wavelengths between 500 to 600 nanometres. For the 1000 m to 2000 m elevation range the upward and downward transmissions are essentially linear with elevation. 10 There have been a number of theoretical studies of radiative transfer mechanisms within the atmosphere that allow one to estimate sky irradiance and path radiance. These studies assume a horizontally homogeneous atmosphere over a horizontal Lam-bertian surface (ie. a diffusely reflecting surface). Results predict sky irradiance and path radiance as a function of optical thickness and average background albedo. Since optical thickness depends on elevation so do sky irradiance and path radiance. For the simple case of an optically thin, horizontally homogeneous atmosphere over a horizontal Lambertian surface of a single albedo one can show that both sky irradiance and path radiance vary almost exponentially with elevation (Turner & Spencer [1972]). Sky irradiance and path radiance may also be influenced by local topography and albedo. The study site has rugged terrain with elevations ranging from 900 to 2700 m above sea-level (MSL). The variety of ground cover present on the site results in marked variation in local albedo. A simple model of sky irradiance and path radiance which does not consider local attributes may be inadequate for application to actual complex scenes. 3.1.3 I L L U M I N A T I O N Direct solar illumination is the most important in Landsat imagery. A target not in shadow sees an attenuated solar beam with irradiance EQTJCOS i where E0 is the solar irradiance at the top of the atmosphere, Ti is the downward transmission through the atmosphere and t is the angle of incidence at the target. Values of E 0 applicable to the Landsat bands are given in table 1 (Woodham and Lee [1985]). E 0 is relatively constant for a quiet sun with a yearly variation of approximately 3.5 %. Sky irradiance or skylight is poorly understood and documented. Generally for clear sky radiance most studies show the majority of the skylight to come from near the sun position (e.g. Stewart [1980] finds 50 % within 40 degrees of the sun) and with some increase towards the horizon . Skylight is stronger for shorter wavelengths because Rayleigh scattering is proportional to the inverse 4th power of wavelength. For example, a clear sky is bright in the visible part of the spectrum but dark in the infrared part of the spectrum. Both the spectral and spatial distribution of skylight require more study. Because of 11 the dependence upon optical thickness skylight generally decreases as elevation increases. Skylight includes radiation reflected from adjacent terrain that is scattered by the atmo-sphere back to the target. Thus skylight may be locally controlled by the topography and albedo of surrounding terrain. Three skylight illumination models are developed. Two are based on a uniform hemi-spherical source, one with slope alone controlling the amount of illumination and the other with slope and obscuring terrain controlling the amount of illumination. The third model utilizies a spatially nonuniform source with slope and obscuring terrain controlling the amount of illumination. Radiation arising from adjacent terrain, sometimes called mutual illumination, or inter-reflection can be a significant source of illumination. Modeling inter-reflection ex-plicitly is difficult because of the complex interaction between surface shape, surface reflectance factors and illumination characteristics. Components of illumination are il-lustrated in figure 2. 3.1.4 SURFACE PHOTOMETRIC PROPERTIES The bidirectional reflectance distribution function (BRDF) was introduced by Nicodemus et al. [1977] as a unified notation for the specification of reflectance in terms of both the incident and reflected beam geometry. The B R D F is an intrinsic property of the surface material and determines how bright a surface will appear when viewed from a given direction and illuminated from another. A Lambertian or diffuse surface appears equally bright regardless of view angle. When illuminated by a collimated source (eg. the sun) with irradiance E0 measured perpen-dicular to the beam of light arriving with incident angle * , one obtains for points not in shadow Lr = — Eo cos i 7T where Lr is radiance. When illuminated by a hemispherical uniform source (eg. uniform skylight) with radiance L0 over the visible hemisphere one obtains (1 + cos e) Lr = L0p -12 where e is slope (see Ho rn & Sjoberg [1979] for der ivat ion of these equations). Here the dependence on e arises because surface elements see differing amounts of sky depending on surface slope. Th is form assumes a level horizon and does not consider the amount of the sky obscured by surrounding terrain. For the study area topographic posit ion as wel l as slope controls the amount of sky seen. These combined correspond to the simple model of a diffuse surface i l luminated by a col l imated (sun) source and a uni form hemispherical source (uniform skyl ight). Pi-, • r (1 + c o s e ) LR = — EQ COS I + p L0 -7T 2 Radiance is control led by the incident angle of the sun at the target and the slope of the target. Th is is a s impl i f icat ion in many respects; the surface cover is not str ict ly Lambert ian, skylight is not uniformly distr ibuted and surrounding topography also determines the amount of sky seen by a target. A lso this model does not yet include atmospheric transmission factors, path radiance or the dependence of skyl ight on elevation. In choosing and testing a model of the image formation process it is important to keep in m ind the physical real ity that corresponds to your model. Phys ica l models correctly characterize simple worlds and can be elaborated as the need is demonstrated. Before modi fy ing a model, one should ensure that it is sufficiently complete to al low confidence in rejecting assumptions. Other B R D F ' s , i n part icular Minnaert surfaces other than Lamber t ian were not tested by the models used in this research. The focus for this work was to thoroughly test the Lamber t i an assumption by expl ic i t ly model ing sun i l luminat ion, skyl ight i l luminat ion and path radiance. 3.1.5 T O P O G R A P H Y The topography of the study area is described by a digita l elevation model ( D E M ) which assigns elevation values on a regular horizontal grid. If the surface is in the form of z = f (x, y) , then a surface normal vector is -df(x, y) -df(x,y) dx ' dy ' 13 Letting = df(x,y) dx _ df{x, y) dy The normal surface vector can be rewritten as [—p, —q, 1]. (p, q) is called the gradient of the function f(x,y). The two dimensional space of (p, q) is used to represent surface orientation. The reflectance of many surfaces depends only on surface orientation. Surface re-flectance can be written as a function of incident, emergent and phase angles (i , e , g) respectively (see figure 3). For distant illumination and orthographic projection the phase angle g is constant. Letting the light source direction have gradient (p0, qo) then the vector [—po , —qo , 1] points in the direction of the light source. For nadir looking sensors the vector [0,0, 1] points in the direction of the viewer. Expressing the cosine of the angle between two vectors as a normalized dot product of the vectors one obtains l + PPO + qg0 COS I cos e = cos (*) = V I + P 2 + q2 1 yjl+pl + q20 The reflectance map R(p,q), introduced by Horn [1977], determines scene radiance as a function of the surface gradient. A reflectance map is a useful representation because is compiles the relevant information about surface material, light source distribution and viewer position into a single function. Four reflectance maps are computed for the analysis undertaken here, all assuming a Lambertian surface viewed from directly overhead. These correspond to the following scene illumination situations; 1. illuminated by the sun alone, reflectance map proportional to (cos i). 2. illuminated by hemispherically uniform skylight with target slope alone controlling the amount of illumination, reflectance map proportional to ((1+c2°* e>). 14 3. illuminated by hemispherically uniform skylight with the amount determined by target slope and sky obscuring surrounding terrain. 4. the same as 3 , but utilizing a spatially asymmetric skylight model based on em-pirical studies. Synthetic images based on these reflectance models are illustrated in figure 4. 3.2 I M A G E F O R M A T I O N E Q U A T I O N Before giving the image formation equation we should state all our initial simplifying assumptions. • The distance of the sun and the satellite are assumed sufficiently distant that or-thographic projection can be applied. • At any time, the Landsat MSS off-nadir viewing angle is small; therefore, we can neglect the dependence of upwards transmission through the atmosphere on viewing angle. • The content and density of the air particles vary within the atmosphere so that the atmosphere appears to have different layers. We will assume that the atmosphere consists of only one layer. The absorbed and reflected properties are laterally homogeneous. • The ground objects are Lambertian reflectors, that is they appear equally bright from any viewing angle. • There is no cloud in the studied imagery. • Polarization effects are not considered. • There are no additional light sources on the surface. • The skylight is initially represented as a hemisphere with a constant intensity dis-tribution (this is modified later with the other skylight models mentioned in sec. 2.1.3). 15 • radiance arising from adjacent terrain is not considered (modified later) Including path radiance and atmospheric transmission the image formation equation of sec. 2.1.4 becomes for targets not in shadow : •L>m — IT Lm is the image irradiance . , 1 + cos e. Td Ltop cos i + L3ky { J path p is the reflectance factor TUiTd are the upward and downward transmission of the atmosphere Ltop is the solar irradiance at the top of the atmosphere and is constant for each band i is the angle between the surface normal and the sun direction Lsky is the sky radiance e is the angle between the surface normal and the vertical Epath is the path radiance 16 Chapter 4 4.1 S T U D Y S I T E The study site is a 25 km by 21 km area surrounding St. Mary Lake in southeastern B . C . ( latitude N 49 36 30, longitude W 116 11 30). It is a mountainous region with elevations ranging from 930 m to 2730 m above mean sea level (MSL). Forest is the predominant ground cover accounting for 70% of the total area. Major species are, in order of abundance; • Lodgepole pine ( Pinus contorta Dougl.) • Balsam fir ( Abies lasiocarpa [Hook.] Nutt.) • Englemann spruce ( Picea englemannii Parry) • Western larch ( Larix occidentalis Nutt.) • Douglas-fir ( Pseudotsuga menziesii [Mirb.] Franco) One quarter of the forested area is non-commercial alpine forest. Many of the accessible, merchantable stands have been harvested, therefore recent clear cuts cover a significant portion (1/8) of the forested area. Alpine areas (rock, snow, meadows) account for 15% of the total area with rock bluffs not in alpine areas covering an additional 7%. Lakes and rivers cover 5% of the total area. The remaining 3% is occupied by roads, powerline right of ways and open fields. 17 4.2 D A T A S E T 4.2.1 L A N D S A T I M A G E R Y A mult ispectra l image consisting of four spectral bands f rom Landsat 2 was obtained for the area. The image, (frame id 22428 17522) was acquired about 9:52 A M P S T on September 15, 1981. Snow cover is a m in imum for late summer, early fal l images. Dur ing the over flight the sun was at an elevation of 38.0° w i th an az imuth of 145.0°. This image was registered to another Landsat 2 image, which had itself been rectified to the digital elevation model using a technique developed by L i t t l e [1982]. The image was resampled to a 100m x 100m pixel size. Th i s was done to match the pixel size of the in i t ia l version of the dig i ta l forest cover map. 4.2.2 D I G I T A L E L E V A T I O N M O D E L The digita l elevation model was produced by manual ly digit iz ing elevation data along ridges and channels f rom the 1:50,000 Canada Nat iona l Topographic System (NTS) mapsheet 82 F /9 (St. M a r y Lake) (L i t t le [1982]). The ridge and channel structure was represented in i t ia l ly using a Tr iangulated Irregular Network (TIN) (Peucker et al. [1978]). A 100m gr id representation was produced f rom the T I N . G r i d coordinates in the D E M correspond to the Universal Transverse Mercator ( U T M ) map projection used in Canad ian N T S maps. The gradient (p, q) at each point in the D E M is estimated by = / ( x + l , y ) -f(x-l,y)  P Ax = f{x,y + l ) - f { x , y - l )  9 A y where x and y are the D E M grid spacings expressed in the same units as f(x,y). The gradient is used to determine values for cos i and cos e. Thus, these angles are determined f rom average elevations spaced 200m horizontal ly apart. Th is is adequate for describing the surface orientat ion of much of the terrain. Some areas are inaccurately described at this coarse scale. The 930m to 2730m M S L elevation range is quantized into 256 levels y ie ld ing a vert ical increment of 7.03 m/level. Th is was done so that the D E M could be displayed and analysed as an image. 18 4.2.3 D I G I T A L F O R E S T C O V E R M A P S The study site is covered by nine 1:20,000 B .C . Ministry of Forests forest cover maps (82F.058, 059, 060, 068, 069, 070, 078, 079, and 080). These maps have been digitized by the M O F as part of their program to digitize the provincial database of some 7000 maps. The map information is stored in a polygon or vector format. That is, the coordinates of forest type polygon boundaries is the way spatial information is captured. Initially the data was received in a raster (grid) format with 100m pixel size, which determined the pixel size of the image and D E M . Later we received the forest cover maps in a polygon format which we could display at a self selected pixel size. The software to accomplish this was written by Marc Majka at the Laboratory for Computational Vision at U B C . There are potentially 65 attributes attached to each map polygon. Forest type attributes included age, height, % crown closure, species composition, site class, and others. 4.3 M O D E L P A R A M E T E R E S T I M A T I O N M E T H O D The forest cover map is used to select target areas of homogeneous ground cover with a presumed constant reflectance factor. Multiple regression analysis is used for these areas to investigate the correlation between recorded radiance (dependent variable) and topo-graphic and atmospheric attributes computed from the D E M (independent variables). 4.3.1 T A R G E T S E L E C T I O N In selecting target areas one should try to satisfy the model assumptions. In particular some types of ground cover such as snow or bare rock are known to be nonLambertian surfaces. Likewise ground cover with an irregular surface structure such as sparse forest is not desirable because of the complexity added by contained shadows. Eliminating the above still leaves a large selection of cover types to choose from. Because the motivation for this work was to increase forest cover classification accuracy it was felt that using a forest cover type to estimate model parameters would be consistent with this goal. Therefore a forest type of a uniform height and age with a closed canopy would be suggested. Furthermore a forest type consisting of a single species would help to 19 insure spectral homogeneity. Unfortunately most pure single species forest stands are not widely distributed over a range of elevations, orientations and topographic positions. This makes it difficult to determine statistically significant correlations. After much investigation of target radiance statistics and trial regression analysis the target forest type chosen was; • Lodgepole pine (>80% pure) • age : 80—100 yr • height : 10—20 m • crown closure : 60—100 % This particular target had the following topographic attributes ; • elevation range: 1205—2102 m MSL • sun illumination angle range (i) : 21.7°—78.0° • slope range (c) : 7.2°—35.7° To insure homogeneity all mixed pixels (pixels intersected by a polygon boundary) were excluded. Also the target areas were visually edited with the help of 1:12,000 aerial photographs to remove pixels that corresponded to roads, landings, and small bare patches that were not depicted on the 1:20,000 forest cover maps. This resulted in 329 pixels available for regression analysis. 4.4 R E G R E S S I O N A N A L Y S I S Using multiple regression analysis to fit a model to the data has shortcomings. These are of two types, first not meeting the assumptions required by the technique and second shortcomings inherent in the data set. A required assumption is that the independent variables are measured without error. All of the independent variables are derived, at least in part from the D E M . This is an approximation of the surface shape and is not without error. It is difficult to assess 20 the error in the DEM (7 m elevation increments) by comparison to the 1:50,000 NTS topographic map (30 m contour intervals). An estimated vertical error is ± 25 m. The image contains radiometric uncertainty due to sensor performance, quantization level, destriping and resampling. Target areas even after elimination of suspect pixels are still only approximately homogeneous. Furthermore a target with a relatively low albedo such as coniferous forest is not the best choice for estimating transmission factors. Keeping in mind the limitations discussed above, multiple regression analysis is an appropriate technique for testing image formation models and estimating parameters. Using a physically plausible model as a guide the regression analysis is kept as simple as possible. The initial equation tested was . .. (1 + cos e)\ TdLtop(cos t) + Laky I + Epath To investigate the elevation dependence of this equation the data set was stratified into six subsets each representing a elevation range of 150 m. Separate regression analysis were run on each subset. The change in path radiance over these elevations could be estimated from the intercept values of the resulting regression equations, similarly the transmission factors and sky irradiance as a function of elevation could be estimated. Unfortunately there were not enough data points over a sufficiently large range of the independent variables to yield significant or consistent regressions for all the subsets. Another approach was tried which linearized the equation. Tu , Tj , Lsky and Epath, were approximated by a linear function of elevation plus a constant. Tu = a + b(z) Leky = f-g{z) Td = c + d(z) E p a t h = h - j(z) This results in an equation with seven terms involving . . . . , . , . (1 + cose) (1 + cose) 2 (1 + cose) cos i , [z)cos i , [zfcos i , ^  - '- , (z)- - '- , [zyi - '- , (z) Regression analysis showed that not all of these terms were significantly correlated (at the 99% significance level) to target radiance. The coefficients derived for the signif-icant variables were difficult to relate to a physical model. The variables found to 21 be significant varied widely from band to band, probably because of the high corre-lation which existed between many of the independent variables. Other combinations of cos i , (z)cos i , ( 1 + c 2 °* e ) } (z) 1+c°s e , z , z2 were tried, before selecting four vari-ables that meet the criteria of physical plausibility, high coefficient of determination and significance at the 99% confidence level. The four variables were (1 + cos e) 2 COS l , , Z , Z This corresponds to a physical model where target radiance depends upon solar il-lumination angle, target slope, and elevation, with no coupling between these variables. Solar illumination angle is the most important factor for all bands, followed by elevation. Surprisingly there was a negative correlation between target radiance values and uniform skylight illumination controlled by slope. This was puzzling and suggested that the model was unrealistic. Therefore a nonuniform skylight distribution was numerically constructed based on the work of Stewart[1984]) and Steven[1977]). The amount of this nonuniform distribution visible for each pixel was computed (taking into account obscuring terrain). The appendix explains the method used to accomplish this. The regression analysis was repeated with cos i , z , z2 and calculated nonuniform skylight as the independent variables. The results were similar, showing a negative correlation between target radiance and calculated skylight. This result led to the hypothesis that the negative correlation of target radiance with skylight was actually a positive correlation with reflection from surrounding terrain. In other words targets surrounded by sky obscuring terrain were brighter due to reflection from surrounding terrain. To test this idea a third "skylight" model was computed which was in essence the ratio of actual slight illumination fo e a c h • j /. amount of a uniform skylight distribution total possible skylight illumtnatton r \ J 0 visible taking into account obscuring terrain). Details of this calculation are given in the Appendix . Regression analysis again yielded similar results. The coefficients of determination and standard errors obtained for various combinations of independent variables are illus-trated in table 2 for band 5. The coefficients of determination achieved for band 5 were generally higher than band 4 and lower than band 6 or 7. 22 4.5 C O R R E C T I O N P A R A M E T E R S The final image formation model chosen for correction purposes is based on an equation of four variables namely, cos i , z , z2 , and the calculated amount of a uniform "skylight" illumination model visible to each pixel taking into account obscuring terrain (i.e. the ratio «^ ^<^^ <wJr /^^ J"^i^ S^ <0^ ) hereafter called topo. All variables were significantly correlated with target radiance at the 99% level. Examination of regression residuals as a function of the independent variables indicated no discernible trends. The consistent negative correlation between calculated "skylight" for all models and target radiance (after correction for sun illumination angle and elevation) suggested that inter-reflection from adjacent terrain was a small but significant factor. Of the three models, topo is the most direct and unbiased measure of the amount of terrain viewed for each pixel. The topo variable resulted in a slightly higher coefficient of determination and a slightly lower standard error, for all bands than either of the other "skylight" illumination models. The equations derived from the regression analysis for the chosen variables are as follows. B A N D 4 R2 = 0.46 SE = 0.014 Lm = 0.541 + 0.0781 cos i - 0.000379 z + 0.000000110 z2 - 0.0185 topo B A N D 5 R2 = 0.60 SE = 0.012 Lm = 0.350 + 0.0903 cos i - 0.000276 z + 0.0000000862 z2 - 0.0193 topo B A N D 6 R2 = 0.71 SE = 0.024 Lm = 0.856 + 0.239 cos i - 0.000779 z + 0.000000218 z2 - 0.0456 topo B A N D 7 R2 = 0.74 SE = 0.065 LM = 2.25 + 0.693 cos i - 0.00210 z + 0.000000571 z2 - 0.139 topo Where Lmis expressed in ( c tffi[« r) > z in metres MSL and cos i and topo are unitless. Topo is the ratio actual uniform .kyUaht illumination R2 • t h c o e f f i c i e n t o f determination and SE is the standard toted possible uniform skylight illumination 23 Descriptive measures of these variables for the target area are as follows. Descriptive Statistics for 329 cases V A R I A B L E M I N I M U M M A X I M U M M E A N S T D D E V S T D D E V A F T E R C O R R E C T I O N Band 4 0.213 0.306 0.254 0.0190 0.0140 Band 5 0.133 0.220 0.172 0.0182 0.0115 Band 6 0.165 0.432 0.287 0.0445 0.0239 Band 7 0.298 1.020 0.683 0.1276 0.0648 Cos i 0.208 0.929 0.616 0.1600 Elevation 1204 2081 1620 194 Topo 0.176 0.980 0.787 0.1606 The coefficients from the above equations were used to correct the entire image by inverting the equation. Band 4 for example LCOrrecud = L- 0.541 - 0.0781 cos i + 0.000379 z - 0.000000110 z2 + 0.0185 topo This calculation was done in floating point mode then converted to integer values. An constant offset was added if necessary to ensure final image values between 0 and 255. 24 Chapter 5 5.1 Assessment of C o r r e c t i o n A c c u r a c y The motivation for correcting the imagery is to improve the accuracy of ground cover classification. Comparison of computer classified imagery with the forest cover map is the primary method of verifying improvement. Visual comparison of the original and corrected image is also used to judge the validity of the correction process. A third method is comparison of experimentally derived atmospheric parameters with published values. 5.1.1 Atmospheric Parameters The chosen image formation equation used for correction does not directly estimate atmospheric transmission, sky irradiance or path radiance. This makes quantitative comparison with published values difficult. The variation of target radiance with model variables can be qualitatively compared to the expected behavior due to the atmosphere. All three "skylight" illumination models exhibit a significant but negative correlation with target radiance. This suggests that skylight is not important for this particular image and target, or it is closely coupled with sun angle and elevation and can not be adequately separated using regression analysis. Evaluation of the inter reflection or mutual illumination parameter is difficult because of the lack of controlled field studies of this phenomenon. In magnitude it had the least effect on the correction process for all bands. For the target areas, the sun angle (cos i) correction has a greater magnitude than the elevation correction for all bands. The difference in magnitude between the sun angle 25 correction and the elevation correction increases with wavelength. This may result from the greater transparency of the atmosphere to longer wavelength radiation or to stronger reflection from the target for longer wavelength radiation. The estimated dependence of target radiance on elevation is shown in figure 5. Gener-ally target radiance initially decreases with elevation, then increases but does not regain the starting value, except for band 5. This derived relationship between target radiance and elevation combines the elevation effects of path radiance, sky irradiance, and atmo-spheric transmission. The general decrease in target radiance with elevation suggests that the effect of path radiance and sky irradiance dominates over the effect of atmo-spheric attenuation. There is also the possibility of a systematic change in the reflectance of the target pixels with elevation due to ecological factors. Field investigation of this possibility was not undertaken. 5.1.2 V i s u a l Comparison Colour composites of the original image and the corrected image are shown in figures 6 and 7. The colour assignment approximates colour infrared photography with band 4 displayed as blue, band 5 as green, and band 7 as red. The display intensity and contrast were independently adjusted for each image to achieve the best depiction of detail, therefore individual pixel values are not directly comparable between images. The most evident change on the corrected image is the equalization of visual appear-ance for all terrain aspects. This is most easily seen by comparison of the southeast slopes to the northwest slopes of the mountain ridges in the lower left portion of the image. The corrected image does not show alpine areas (bright cyan colour) on the northwest slope of these ridges. This results from the stronger correction for band 7 compared to band 4 or 5 for those areas of high sun angle. This consequence of the correction process is not surprising when one considers the difference in reflectance properties between the forested target areas and alpine areas. Alpine areas are highly reflective compared to forested areas for bands 4 and 5. Shadows have mostly been removed on the corrected image but these areas display little detail. This is not surprising because the original image values are often constant in shadows. 26 Some small topographic features have not been corrected entirely. This is probably due to the difficulty in capturing these abrupt changes in slope with the 200 m spacing used to calculate the gradient [p, q]. This can be seen along ridge tops and ravines in the northeast portion of the image. The elevation correction has the effect of brightening pixels at mid-elevations. Vi-sually it is hard to judge whether this is more correct than the original image. The classification results help to judge this. The most significant visual result is the increased ease in identifying clearcut areas. In particular the logging clearcuts southeast of St. Mary Lake and in the southwest corner are better delineated on the corrected image. 5.1.3 Classification Results The classification is relatively simple with four classes, forest, clearcut, water, and alpine. Although simple, this classification, if accurate, is useful because satellite data is presently being used to map forest depletions (logging clearcuts and fires). The truth map was constructed from the forest cover map. Clearcuts were defined as areas with the following attributes; age < 20 yr or, height < 10m or, crown closure < 20%. Thus this clearcut class does not strictly correspond to logging clearcuts. The rest of the forested areas were aggregated into one forest class. Water and alpine are taken directly from the forest cover map. A null class of powerline right of ways, urban areas, clearings, and nonproductive forest areas (burns, brush etc.) was also defined for the truth map. The colour coded ground truth map is shown in figure 8, as are the classifications of the original and corrected images. The classification was accomplished with a nearest centroid classifier using all four MSS bands. The resulting classification was filtered with a 3x3 maximum frequency operator which replaced the central pixel value with the most frequent value in the 3x3 window. This was done to remove the salt and pepper appearance of the classified image and results in a more map-like representation. Clearcuts above 2260 m MSL were redefined as alpine and alpine zones below 2260 m MSL were redefined as clearcuts, on the truth map and the classified imagery. The classification of the original image identifies much of the southeast aspect slopes 27 as clearcuts and much of the northwest aspect slopes as water. Overall there is too much clearcut and water, alpine areas are approximately correct and there is not enough forest. The classification of the corrected image is a closer match to the truth map. Except for small areas of incompletely corrected shadows the water classification is quite accu-rate. The alpine classification is approximately the same as the uncorrected case ( 69% correctly classified). Forest classification is much better than the uncorrected case ( 80% correctly classified). Clearcuts are better defined with much less forest wrongly classified as clearcuts. There is still a tendency to wrongly identify southeast aspect slopes as clearcut indicating that not all of the topographic effect has been removed. Overall (excluding the null class) 54% of the pixels were correctly classified from the original image versus 72% for the corrected image. The confusion matrices for the original, corrected for sun angle, corrected for sun angle and elevation, and corrected for sun angle, elevation and topo variable are found in table 3. Examination of these matrices leads to the following conclusions: • the sun illumination angle correction achieves the greatest increase in classification accuracy. • the elevation correction in addition to the sun illumination angle correction de-creases the misclassification of forest and clearcuts as water. • the topo variable correction has little effect on classification results. It is difficult to quantitatively judge the effectiveness of the correction method by comparison of classifications of the original and corrected image alone. Ideally one should have an area of flat terrain with similar ground cover to use as a standard for comparison. 28 Chapter 6 6.1 S U M M A R Y The use of satellite remote sensing in forestry is expanding. The primary methods used to extract information from satellite imagery are human interpretation and spectral clas-sification by computer. Imagery of mountainous terrain is often dominated by features due to topography and to a lesser extent the atmosphere. These effects make the correct identification of surface cover difficult. The image formation process must be understood before unwanted factors can be corrected for. By explicitly modeling the process we test and improve our understanding. An image formation model is developed and tested on Landsat MSS data over a mountainous region. The primary assumptions of this model are, the ground objects are diffuse (Lambertian) reflectors, and the atmosphere is horizontally homogeneous. Topographic and atmospheric attributes are computed from a digital elevation model (DEM). These include, sun illumination angle, elevation and sky visibility. A digital forest cover map is used to select target areas of homogeneous forest cover. The forest type selected, after much experimentation, was pure Lodgepole pine, 80 to 100 yrs old, 10 to 20m in height with 60 to 100% crown closure. This forest type was chosen because of it's relatively high crown closure, single species composition, and a wide distribution range of elevation, sun illumination angle, and slope. Multiple regression analysis guided by the criteria of physical plausibility, high coef-ficient of determination and significance at the 99% level resulted in an image formation equation of four variables; sun illumination angle, elevation and elevation squared and calculated skylight. 29 Sun illumination angle had the largest effect on target radiance for all bands. The magnitude of this effect increased with wavelength. There are two probable causes of this increase with wavelength, the increased transparency of the atmosphere for longer wavelength radiation and stronger reflection at longer wavelengths for the target areas. Elevation was the next most important determinant of target radiance. Target radiance initially decreased with elevation then increased, but not back to the original level, except for band 5 (see figure 5). The shape of these curves suggests that path radiance is the dominant factor at low elevations while atmosphere attenuation becomes relatively more important at higher elevations. Three skylight illumination models were tested. Two were based on a uniform dis-tribution of illumination, and one was based on an empirical sun centred asymmetric illumination distribution. All three were negatively correlated with target radiance. This suggests that inter reflection or mutual illumination from adjacent terrain may be a small but significant factor for the target areas. Skylight may not be a significant factor for this particular image and forest target or it's effect may be coupled with sun illumination angle and elevation and therefore not separable using simple regression analysis. Parameters derived from the regression analysis of the target areas were used to cor-rect the image for sun illumination angle, elevation and a topographic factor tentatively identified as inter reflection. Visual comparison of the corrected and original images (false colour composite of band 4, 5 and 7 ) shows that much of the topographic effects have been removed. Logging clearcuts are much easier to delineate. Alpine areas are not correctly portrayed; they are missing from the northwest aspect slopes on the corrected image. Comparisons of computer aided classifications of the original and corrected images with the forest cover map show an increase in pixels correctly classified from 54% to 72%. Most of this increase is due to the sun illumination angle correction. 6.2 C O N C L U S I O N This research has demonstrated that image classification for spectrally distinct classes can be improved by utilizing an image formation model to investigate and correct topographic 30 effects. The results of this analysis need to be validated for the same location on different dates, other locations, and targets of other ground cover classes. The method also needs to be evaluated by attempting more detailed classifications. The data required for this type of analysis is improving in quality and availability. Government agencies and resource industries are digitizing their map databases. This makes the operational integration of this information with digital imagery and DEM's more feasible. Digital elevation models at an effective resolution of approximately 25 m grid spacing are being created by the British Columbia and Alberta governments. The advent of Landsat 5 Thematic Mapper data with 30 m pixel size and additional spectral bands has dramatically increased the information content of publically available satellite imagery. This data is soon to be available in a geocoded form, thereby making integration with map information easier. The SPOT satellite adds the capability of off nadir viewing. This will be useful for investigating atmospheric parameters, particularly if simultaneous views of the same area are acquired by Landsat 5 and SPOT. Combining information on the appearance, shape, and condition of the land results in a more complete and useful description. 31 R E F E R E N C E S • Ahem, F.J., D.G. Goodenough, S.C. Jain, V.R. Rao & G. Rochon (1977), "Use of Clear Lakes as Standard Reflectors for Atmospheric Measurements", Proc. 11th Int. Symp. on Remote Sensing of Environment, pp 731-755. • Horn, B.K.P.(1977), " Understanding Image Intensities" Artificial Intelligence(8) :201-231. • Horn, B.K.P. & B.L.Bachman (1978), "Using Synthetic Images to Register Real Images with Surface Models", Comm. ACM, 21 (11): 914-924. • Horn, B.K.P. (1978), " The Position of the Sun", AI-WP-162, Artificial Intelligence Laboratory, MIT, Cambridge, MA. • Horn, B.K.P. & R.W.Sjoberg (1979), "Calculating the Reflectance Map", Applied Optics, 18 (11) : 1770-1779. • Justice, CO., S.W.Wharton &; B.N.Holben (1980),"Application of Digital Terrain Data to Quantify and Reduce the Topographic Effect on Landsat Data", NASA TM 81988 Godard Space Flight Centre, Greenbelt , Maryland. • Kimes, D.S., J.A.Smith & K.J.Ranson (1980), "Vegetation Reflectance Measur-ments as a Function of Solar Zenith Angle", Photgramm. Engng. Remote Sensing 46 (12) : 1563-1573. • Landsat Data Users Handbook, revised edition. Available from Branch of Distri-bution, U.S. Geological Survey, Arlington, VA 22202 (1979). • Lee, T.K. (1983), "Modeling and Estimating Atmosphere Effects in Landsat Im-agery", MSc Thesis, Dept. of Computer Science, U.B.C., Vancouver, B.C.. • Little, J.J. (1982), "Automatic Registration of Landsat Images using Features De-rived from Digital Terrain Models", Proc. Workshop on Computer Vision: Repre-sentation and Control, ppl78-184, Rindge, New Hampshire • Minneart, M. (1941), "The Reciprocity Principle in Lunar Photometry "Astrophys. J. (93) : 403-410. 32 Nicodemus, F.E., J.C.Richmond, J.J.Hsia, I.W.Ginsburg & T.Limperis (1977), "Geometrical Considerations and Nomenclature for Reflectance", NBS Monogram 160, National Bureau of Standards, Washington, D.C. Peucker, T.K., R.J. Fowler, J.J. Little & D.M. Mark (1978) "The Triangulated Irregular Network", Proc. of the Digital Terrain Models Symp., American Society of Photogrammetry, pp.516-532, St. Louis, May 9-11, 1978. Sjoberg, R.W. & B.K.P.Horn (1983), "Atmospheric Effects in Satellite Imaging of Mountainous Terrain", Applied Optics, (22) : 1702-1716. Smith, J.A., T.L.Lin & K.J.Ranson (1980), "The Lambertain Assumption and Landsat Data", Photogramm. Engng. Remote Sensing, 46 (9):1183-1189. Steven, M.D. (1977), "Standard Distributions of Clear Sky Radiance", Quarterly Journal of the Royal Meteorological Society (103) 457-461. Stewart, K.A. (1984), "The Angular Distribution of Diffuse Solar Radiation Over the Sky Hemisphere", MSc Thesis, Dept. of Geography, U.B.C., Vancouver, B.C.. Teillet, P.M., B. Guidon & D.G. Goodenough (1982), "On the Slope-Aspect Cor-rection of Multispectral Scanner Data", Canadian Journal of Remote Sensing (8) 84-106. Thekaekara, M.P. (1970), "Proposed Standard Values of the Solar Constant and Solar Spectrum", Journal of Environmental Sciences (13) 6-9. Turner, R.E. & M.M. Spencer (1972), "Atmospheric Model for Correction of Space-craft Data", Proc. 2nd International Symposium on Remote Sensing of Environ-ment, pp.895-934. Valley, S.L. (ed.) (1965), Handbook of Geophysics and Space Environments, McGraw-Hill. Wong, F. (1984), "A Unified Approach to the Geometric Rectification of Remotely Sensed Imagery", Computer Science TR-84-6, U.B.C., Vancouver, B.C.. 33 • Woodham, R.J. (1980), "Using Digital Terrain Data to Model Image Formation in Remote Sensing", Proc. SPIE(238)361-369. • Woodham, R.J. & T.K.Lee (1985), "Photometric Method for Radiometric Cor-rection of Multispectral Scanner Data", Canadian Journal of Remote Sensing 11 (2):132-142. 34 Appendix Calculation of Sky Irradiance Models The simplest sky irradiance model is that of a uniformly radiating hemisphere, with the amount visible to a surface element controlled by the surface slope. For an ideal (lossless) Lambertian surface one obtains _ 1 + cos e Lr = L0  A more realistic model would be a sun centred spatial distribution of irradiance, and would take into account the amount of sky obscured by surrounding terrain. A generalized distribution of sky irradiance for temperate latitudes with a sun elevation of 38" was constructed based on the work of Stewart [1984] and Steven [1977]. This was in the form of a contour map of relative irradiance values and is illustrated by figure 9. Forty eight point sources were arranged to simulate this distribution. The coordinates of these point sources are given in table 4. For each point source the Lambertian reflectance map was computed (ie. L=cos i, or L=0 if in shadow for that particular source). All reflectance maps were then summed together with the maximum value set to 255 for display purposes. Figure 4(b) shows the resulting reflectance map for this illumination situation. A third sky irradiance model was constructed. It corresponded to illumination by a uniformly radiating hemisphere taking into account the obscuring effect of surrounding terrain. This model required a tesselation of the sky. This was accomplished with a 35 hexagonal arrangement of 109 point sources as illustrated by figure 10. This pattern was chosen because rotational symmetry reduced the number of azimuth angle computations. The coordinates of these 109 source are given in table 5. The radial coordinates measured from this pattern where transformed by (1 — cos 0)? so that equal areas represent equal solid angles. As above the Lambertian reflectance map was calculated for each source point then summed together. Because of the weighting by cos i illumination sources perpendicular to the target pixel contribute more than sources at a grazing angle. Thus this process does not calculate the ratio but rather the ratio <f»* *W^ j * ™ ™ * ™ for ground seen total possible skylight illumination a Lambertian surface, taking into account shadowing by adjacent terrain. 36 S O L A R I R R A D I A N C E A T T H E T O P O F T H E A T M O S P H E R E Landsat MSS: Wavelength Interval (pm) Percentage of Solar Constant Solar Irradiance mW/cm2 band 4 .5 - .6 13.1 17.7 band 5 .6- .7 11.2 15.1 band 6 .7 - .8 9.1 12.4 band 7 .8 - 1.1 18.4 24.9 Landsat T M : band 1 band 2 band 3 band 4 band 5 band 7 Wavelength Interval (/xm) .45 - .52 .52 - .60 .63 - .69 .76 - .90 1.55 - 1.75 2.08 - 2.35 Percentage of Solar Constant 10.2 10.3 6.6 10.8 3.3 1.6 Solar Irradiance (mW/crri2) 13.9 13.9 8.9 14.6 4.5 2.1 based on solar constant of 135.3 mW/crri2.This value is considered accurate to ±2.1 mW/cm? for a quite sun at the mean sun-to-earth distance. It is estimated that the solar constant varies from 130.9 mW/cm 2 at aphelion to 139.9 mW/cm3 at perihelion. Table 1. Values of solar irradiance at the top of the atmosphere for Landsat satellites (after Woodham & Lee [1984]).Derived from standard A N S I / A S T M E 490 - 1973 "Solar constant and air mass zero solar spectral irradiance tables", proposed in (Thekaekara Jl970]).The value for each Landsat band is obtained by integrat-ing the standard solar spectral irradiance curve over the specified wavelength interval. 37 M U L T I P L E REGRESSION RESULTS F O R VARIOUS COMBINATIONS O F I N D E P E N D E N T V A R I A B L E S F O R B A N D 5 SE cos i z z'2 topo sky ±±£f££ z c o 8 i z topo z sky Confidence Level Attained (%) 60 1.48 99.9 99.9 99.9 99.9 60 1.49 99.9 99.9 99.9 99.9 59 1.48 99.9 99.9 99.9 99.9 58 1.52 99.9 99.9 92.7 99.9 99.2 57 1.52 99.9 99.9 99.9 99.9 57 1.53 99.9 99.9 99.9 99.9 57 1.53 99.9 99.9 99.9 99.9 57 1.57 99.9 99.5 99.9 56 1.54 99.9 16.6 98.9 56 1.54 99.9 37.6 99.1 54 1.58 99.9 84.7 99.9 Where R 2 is expressed as a percentage and SE is expressed in digital image value units. Sky is asymmetric skylight illumination. Topo is the ratio **>'•>y »'»y/"«t.'°''. . * J J «> f total potable v n i / o r m tkyltght i l l u m i n a t i o n Table 2. Regression analyses results for band 5. 38 O R I G I N A L I M A G E Image Classification Truth Map Pixel Count Pixel Count FOREST C L E A R C U T WATER ALPINE 17291 11492 4421 498 33702 FOREST 3638 7087 1534 0 12259 C L E A R C U T 13 145 259 5 422 WATER 63 0 943 2193 3199 ALPINE 21005 18724 7157 2696 49582 TOTALS I M A G E C O R R E C T E D F O R S U N I L L U M I N A T I O N A N G L E Image Classification Truth Map Pixel Count Pixel Count FOREST C L E A R C U T WATER ALPINE 26823 6167 280 432 33702 FOREST 5548 6608 103 0 12259 C L E A R C U T 47 113 260 2 422 WATER 951 0 136 2112 3199 ALPINE 33369 12888 779 5546 49582 TOTALS I M A G E C O R R E C T E D F O R S U N I L L U M I N A T I O N A N G L E A N D E L E V A T I O N Image Classification Truth Map Pixel Count Pixel Count FOREST C L E A R C U T WATER ALPINE 26726 6516 45 414 33702 FOREST 5872 6375 11 0 12259 C L E A R C U T 102 50 267 3 422 WATER 947 0 91 2161 3199 ALPINE 33647 12941 414 2578 49582 TOTALS I M A G E C O R R E C T E D F O R S U N I L L U M I N A T I O N A N G L E , E L E V . A N D T O P O Image Classification Truth Map Pixel Count Pixel Count FOREST C L E A R C U T WATER ALPINE 26966 6249 68 416 33702 FOREST 5940 6301 17 0 12259 C L E A R C U T 104 49 265 4 422 WATER 904 0 100 2195 3199 ALPINE 33914 12599 450 2615 49582 TOTALS Table 3. Confusion matrices from image classification (null class excluded). 39 SUN C E N T R E D A S Y M M E T R I C S K Y I R R A D I A N C E M O D E L for Sun Position : azimuth 145.0° elevation 38.0° number AZIMUTH ELEVATION number AZIMUTH ELEVATION 1 140 10 2 150 10 3 124 12 4 166 12 5 111 10 6 179 10 7 091 11 8 199 11 9 080 24 10 210 24 11 064 14 12 226 14 13 043 26 14 247 26 15 023 20 16 267 20 17 001 18 18 289 18 19 338 17 20 312 17 21 352 53 22 298 53 23 053 51 24 237 51 25 036 75 26 254 75 27 083 49 28 207 49 29 100 73 30 190 73 31 121 63 32 169 63 33 105 44 34 185 44 35 103 22 36 187 22 37 130 56 38 160 56 39 127 42 40 163 42 41 117 30 42 173 30 43 138 42 44 152 42 45 140 30 46 150 30 47 131 21 48 159 21 Table 4. Point source coordinates for asymmetric sky irradiance model. 40 U N I F O R M T E S S E L A T I O N O F A H E M I S P H E R E AZIM E L E V AZIM E L E V AZIM E L E V 1 0 10 37 30 37 073 14 35 2 60 10 38 90 37 074 46 35 3 120 10 39 150 37 075 74 35 4 180 10 40 210 37 076 106 35 5 240 10 41 270 37 077 134 35 6 300 10 42 330 37 078 166 35 7 0 28 43 30 06 079 194 35 8 60 28 44 90 06 080 226 35 9 120 28 45 150 06 081 254 35 10 180 28 46 210 06 082 286 35 11 240 28 47 270 06 083 314 35 12 300 28 48 330 06 084 346 35 13 00 45 49 19 50 085 11 18 14 60 45 50 41 50 086 49 18 15 120 45 51 79 50 087 71 18 16 180 45 52 101 50 088 109 18 17 240 45 53 139 50 089 131 18 18 300 45 54 161 50 090 169 18 19 00 60 55 199 50 091 191 18 20 60 60 56 221 50 092 229 18 21 120 60 57 259 50 093 251 18 22 180 60 58 281 50 094 289 18 23 240 60 59 319 50 095 311 18 24 300 60 60 341 50 096 349 18 25 00 75 61 19 04 097 23 22 26 60 75 62 41 04 098 37 22 27 120 75 63 79 04 099 83 22 28 180 75 64 101 04 100 97 22 29 240 75 65 139 04 101 143 22 30 300 75 66 161 04 102 157 22 31 30 64 67 199 04 103 203 22 32 90 64 68 221 04 104 217 22 33 150 64 69 259 04 105 263 22 34 210 64 70 281 04 106 277 22 35 270 64 71 319 04 107 323 22 36 330 64 72 341 04 108 337 22 109 00 90 Table 5. Position coordinates for the uniform tesselation of a hemisphere. 41 ELEVATION (m) Figure 1. Variation of optical depth and upward and down-ward transmission with elevation (after Sjoberg [1982]). The data is from the U.S. Standard Atmopsphere for molecular and aerosol scattering components (Valley[1965]) for wavelengths between 0.5 and 0.6 nm. 42 Figure 2. Components of illumination. Figure 4. Reflectance maps for the following illumination con-ditions viewed from directly overhead. (a) sun illumination (cos i) (b) asymmetric skylight illumination (sun centred) (c) uniform skylight illumination with amount controlled by slope alone (d) uniform skylight illumination with amount controlled by slope and obscuring terrain 44 Figure 4. Reflectance maps for the following illumination con-ditions viewed from directly overhead. (a) sun illumination (cos i) (b) asymmetric skylight illumination (sun centred) (c) uniform skylight illumination with amount controlled by slope alone (d) uniform skylight illumination with amount controlled by slope and obscuring terrain 45 0.3 0.1 0.2 1 1 1 • • • • ELEVATION < W m • • • ^ RANGE • — • • • OF TARGET • • • • PIXELS • • • • * • • • • • • • • • • • • • • • • ••* « • • , • • • •• • BAND 5 " _ • • • _ \ X . BAND 4 + • • • • • • \ BAND 6 • • • • • • • • • • • • • - \ B A N D 7 y • • • • • 9 0 0 1 2 0 0 1 6 0 0 2 0 0 0 2 4 0 0 2 8 0 0 ELEVATION (m MSL) Figure 5. Estimated variation of target radiance with elevation after correction for sun illumination angle and topo variable. Ra-diance values normalized to zero at 1200m MSL for all bands. 46 Figure 6. Original image (false colour composite, band 4 dis-played as blue, band 5 as green and band 7 as red). 47 Figure 7. Image corrected for sun illumination angle, elevation and topo variable. Same colour assignment as figure 6. 48 B C Figure 8. Classification, results after filtering with a 3x3 maxi-mum frequency filter. Green is forest, red is clearcut class, blue is water, yellow is alpine and black is a null class. (a) original image classification (b) truth map derived from forest cover map (c) corrected image classification 49 Figure 9. Asymmetric sun centred skylight illumination dis-tribution for sun position; elevation 38°, azimuth 145" based on Stewart (1984] and Steven [1977]. Contours are relative irradi-ance. Point source positions simulate distribution. 50 Figure 10. Uniform tesselation of a hemisphere. Centres of hexogons correspond to point source positions. Radial angle trans-formed by (1 — cos 0)» so equal areas represent equal solid angles. 51 


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items