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  M A L C O L M HAIG GRAY B.Sc, University of British Columbia, 1978 A THESIS SUBMITTED IN PARTIAL F U L F I L L M E N T OF T H E REQUIREMENTS F O R T H E D E G R E E OF M A S T E R 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  In presenting  this  thesis i n partial  f u l f i l m e n t of the  r e q u i r e m e n t s f o r an a d v a n c e d d e g r e e a t t h e of B r i t i s h Columbia, I agree that it  freely  the L i b r a r y s h a l l  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 .  agree t h a t permission f o r extensive for  University  s c h o l a r l y p u r p o s e s may  for  financial  of  The U n i v e r s i t y o f B r i t i s h 1956 Main M a l l V a n c o u v e r , Canada V6T 1Y3  Columbia  Date  ii  my  It is thesis  s h a l l n o t be a l l o w e d w i t h o u t my  permission.  Department  thesis  be g r a n t e d by t h e h e a d o f  copying or p u b l i c a t i o n of t h i s  gain  further  copying of t h i s  d e p a r t m e n t 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 . understood that  I  make  written  ABSTRACT The radiometry of satellite imagery is influenced by ground cover, local topography, and atmosphere. In order to increase the accuracy 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 illumination 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 analysis. Results of this analysis indicate that solar illumination angle has the largest effect on target pixel irradiance followed by atmosphere 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 imagery for topographic and atmospheric effects. Visual assessment of the corrected imagery indicates that many but not all of the topographic effects have been reduced. Comparisons between computer classified imagery and the forest cover map show an improvement in correctly classified pixels from 54% for the original image to 72% for the corrected image.  iii  CONTENTS  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 2.2 Research Objective 2.3 Methodology 2.4 Previous Work  4 4 5 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 4.2 Data Set 4.2.1 Landsat Imagery 4.2.2 Digital Elevation Model 4.2.3 Digital Forest Cover Maps 4.3 Model Parameter Estimation Method 4.3.1 Target Selection 4.4 Regression Analysis 4.5 Correction Parameters  17 18 18 18 19 19 19 20 23  CHAPTER 5  25  5.1 Assessment of Correction Accuracy 5.1.1 Atmospheric Parameters 5.1.2 Visual Comparison 5.1.3 Classification Results  25 25 26 27  CHAPTER 6  29  6.1 Summary 6.2 Conclusion  29 30  References  32  Appendix  35  Tables  37  Figures  42  v  LIST OF  TABLES  1  Solar Irradiance at the Top of the Atmosphere  37  2  Regression Results for Band 5  38  3  Confusion Matrices for Classification Results  39  4  Point Source Coordinates for the Asymmetric Sky Model  40  5  Point Source Coordinates for the Uniform Tesselation of a Hemisphere  41  VI  LIST OF  FIGURES  1  Variation of Optical Depth and Upward and Downward Transmission with Elevation  42  2  Components of Illumination  43  3  Definition of Incident Angle, Emergent Angle and Phase Angle  44  4  Reflectance Maps  5  Target Radiance vs Elevation  46  6  Original Image  47  7  Corrected Image  48  8  Classification Results  49  9  Asymmetric Sky Illumination Model  50  44,45  10 Uniform Sky Illumination Model  vii  51  ACKNOWLEDGEMENTS 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 programming 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 appreciation is extended to fellow graduate students at the Laboratory for Computational Vision for their software tools made freely available 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  INTRODUCTION  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: better 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 orientation 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 reflectance 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 characteristics 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 distribution 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 effects. 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 imageformation 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  THESIS  STATEMENT  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 problematic. 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  RESEARCH OBJECTIVE  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 atmospheric parameters for the time of satellite overpass.  4  2.3  METHODOLOGY  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 situation 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 considered 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 ( D E M ) . 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 radiometry as a function of direct solar illumination angle, depth of atmosphere and sky illumination. The equation is linearized and multiple linear regression analysis is used estimate model parameters. The coefficient of determination (R ), 2  standard error, and residuals  were examined to judge the suitability of the various models investigated. Physical plausibility 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  PREVIOUS  WORK  Lambertian (diffuse) reflectance models have been widely tested. Smith et al. [1980]) investigated the Lambertian assumption for Landsat MSS data for Ponderosa pine (Pinus 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 targets 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. Moderate 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 transmission 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 relief. 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 acceptable 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  IMAGE FORMATION  MODEL  The image formation process must be understood before unwanted factors can be corrected 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 properties 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  in the image plane where x = u and y = v. Therefore coordinates (x,y)  (u,v)  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  ATMOSPHERE  For the model developed here the atmosphere is assumed to be an optically thin, horizontally 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 transmission 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)  =  r e- * 0  z/H  where r = r(0) is the optical thickness at sea level and parameter H is the scale height. r  0  For a horizontally homogeneous atmosphere the upward and downward transmissions become T {z)  =  u  T {z) d  =  e~ W r  e-'W' ^ 03  Figure 1 illustrates the variation of optical thickness and upward and downward transmissions 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 Lambertian 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  ILLUMINATION  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  E  0  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 is relatively constant 0  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 4  th  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 atmosphere 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 hemispherical 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 explicitly is difficult because of the complex interaction between surface shape, surface reflectance factors and illumination characteristics. Components of illumination are illustrated 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 E  0  measured perpen-  dicular to the beam of light arriving with incident angle * , one obtains for points not in shadow L  r  = — Eo cos i 7T  where L  r  is radiance.  When illuminated by a hemispherical uniform source (eg. radiance L  0  over the visible hemisphere one obtains L  r  = Lp 0  (1 + cos e) -  12  uniform skylight) with  where e is slope (see H o r n & Sjoberg [1979] for derivation of these equations). Here the dependence on e arises because surface elements see differing amounts of sky depending on surface slope. T h i s f o r m assumes a level horizon and does not consider the amount of the sky obscured b y surrounding terrain. F o r the study area topographic position as well as slope controls the amount of sky seen. These combined correspond to the simple model of a diffuse surface i l l u m i n a t e d by a collimated (sun) source a n d a u n i f o r m hemispherical source (uniform skylight). Pi-,  L = — EQ R  •  COS I  r  + p L  0  (1 +  -  c  o  s  e  )  2  7T  R a d i a n c e is controlled by the incident angle of the sun at the target and the slope of the target. T h i s is a simplification i n m a n y respects; the surface cover is not strictly L a m b e r t i a n , skylight is not u n i f o r m l y distributed a n d surrounding topography also determines the amount of sky seen by a target.  A l s o this model does not yet include atmospheric  transmission factors, p a t h radiance or the dependence of skylight on elevation. I n choosing and testing a m o d e l of the image formation process it is i m p o r t a n t to keep i n m i n d the p h y s i c a l reality that corresponds to y o u r m o d e l . P h y s i c a l models correctly characterize simple worlds and can be elaborated as the need is demonstrated. Before m o d i f y i n g a model, one should ensure that it is sufficiently complete to allow confidence i n rejecting assumptions. O t h e r B R D F ' s , i n particular M i n n a e r t surfaces other t h a n L a m b e r t i a n were not tested by the models used i n this research. T h e focus for this work was to thoroughly test the L a m b e r t i a n assumption by explicitly modeling sun i l l u m i n a t i o n , skylight i l l u m i n a t i o n a n d p a t h radiance.  3.1.5  TOPOGRAPHY  T h e topography of the study area is described by a digital elevation model ( D E M ) which assigns elevation values on a regular horizontal grid. If the surface is i n the f o r m of z = f (x, y) , then a surface n o r m a l vector is -df(x, dx  y)  -df(x,y) '  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 reflectance 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 (p , qo) then the vector [—po , —qo , 1] 0  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 +  COS I  cos e = cos (*)  =  qg  0  VI + P + q 2  2  1 yjl+pl  + q  2  0  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+c  2  14  e  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 empirical studies. Synthetic images based on these reflectance models are illustrated in figure 4.  3.2  IMAGE FORMATION  EQUATION  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 orthographic projection can be applied. • A t 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 distribution (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  T L d  .  top  cos i + L  3ky  , 1 + cos e.  {  J  path  is the image irradiance  L  m  p is the reflectance factor  T Td Ui  L  are the upward and downward transmission of the atmosphere  is the solar irradiance at the top of the atmosphere and is constant for each band  top  i is the angle between the surface normal and the sun direction  Lk s  y  is the sky radiance  e is the angle between the surface normal and the vertical  Ep h is the path radiance at  16  Chapter 4  4.1  S T U D Y SITE  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 1 5 % 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 4.2.1  DATA SET LANDSAT IMAGERY  A m u l t i s p e c t r a l image consisting of four spectral bands f r o m L a n d s a t 2 was obtained for the area. T h e image, (frame i d 22428 17522) was acquired about 9:52 A M P S T on September 15, 1981. Snow cover is a m i n i m u m for late summer, early fall images. D u r i n g the over flight the s u n was at a n elevation of 38.0° w i t h a n a z i m u t h of 145.0°. T h i s image was registered t o another L a n d s a t 2 image, w h i c h h a d itself been rectified to the digital elevation m o d e l using a technique developed b y L i t t l e [1982]. T h e image was resampled to a 1 0 0 m x 100m pixel size. T h i s was done t o m a t c h the pixel size of the i n i t i a l version of the d i g i t a l forest cover m a p .  4.2.2  DIGITAL ELEVATION  MODEL  T h e d i g i t a l elevation m o d e l was produced b y m a n u a l l y d i g i t i z i n g elevation d a t a along ridges a n d channels f r o m the 1:50,000 C a n a d a N a t i o n a l T o p o g r a p h i c S y s t e m ( N T S ) mapsheet 82 F / 9 (St. M a r y Lake) ( L i t t l e [1982]).  T h e ridge a n d channel structure  was represented i n i t i a l l y using a T r i a n g u l a t e d Irregular N e t w o r k ( T I N ) (Peucker et a l . [1978]). A 1 0 0 m g r i d representation was produced f r o m the T I N . G r i d coordinates i n the D E M correspond to the U n i v e r s a l Transverse M e r c a t o r ( U T M ) m a p projection used i n C a n a d i a n N T S maps. T h e gradient (p, q) at each point i n the D E M is estimated b y =  /(x+l,y)  -f(x-l,y) Ax  P  =  f{x,y  9  +  l)-f{x,y-l)  A y  where x a n d y are the D E M grid spacings expressed i n the same units as f(x,y).  The  gradient is used to determine values for cos i a n d cos e. T h u s , these angles are determined f r o m average elevations spaced 2 0 0 m horizontally apart. T h i s is adequate for describing the surface orientation of m u c h of the terrain. Some areas are inaccurately described at this coarse scale. T h e 9 3 0 m to 2730m M S L elevation range is quantized into 256 levels y i e l d i n g a v e r t i c a l increment of 7.03 m/level. T h i s was done so that the D E M could be displayed a n d analysed as a n image.  18  4.2.3  DIGITAL FOREST COVER  MAPS  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  MODEL PARAMETER ESTIMATION  METHOD  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 topographic and atmospheric attributes computed from the D E M (independent variables). 4.3.1  TARGET  SELECTION  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 M S L • 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  REGRESSION ANALYSIS  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 D E M (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 .  TdL (cos top  (1 + cos e)\  ..  t) + L  I + E th  aky  pa  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. T  u  , Tj , L  sky  and Epath,  were approximated by a linear function of elevation plus a  constant. T  u  = a + b(z)  T = c + d(z) d  L E  eky  p a t h  =  f-g{z)  = h-  j(z)  This results in an equation with seven terms involving . . . . , . , cos i , [z)cos i , [zfcos  . (1 + cose) (1 + cose) 2 (1 + cose) 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 significant 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 correlation which existed between many of the independent variables. Other combinations of cos i , (z)cos i , (  1+c  2  °* ) e  }  (z) °  1+c s  , z , z were tried, before selecting four vari-  e  2  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)  COS l ,  2  ,  Z, Z  This corresponds to a physical model where target radiance depends upon solar illumination 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 , z skylight as the independent variables.  2  and calculated nonuniform  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  foe a c h  total  possible  skylight  illumtnatton  r  • j /.  \  amount of a uniform skylight distrib 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 illustrated 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  CORRECTION  PARAMETERS  The final image formation model chosen for correction purposes is based on an equation of four variables namely, cos i , z , z , and the calculated amount of a uniform "skylight" 2  illumination model visible to each pixel taking into account obscuring terrain (i.e. the ratio  «^^^^wJ^/^J"^i^S^^) <  r  <  0 <  hereafter called topo. A l l 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.  BAND 4 L  R = 0.46 SE = 0.014 2  =  m  0.541 + 0.0781 cos i - 0.000379 z + 0.000000110 z  BAND 5 L  R = 0.60 SE = 0.012 2  =  m  0.350 + 0.0903 cos i - 0.000276 z + 0.0000000862 z  BAND 6 L  - 0.0193 topo  2  BAND 7  M  2  R = 0.71 SE = 0.024  = 0.856 + 0.239 cos i - 0.000779 z + 0.000000218 z  m  L  - 0.0185 topo  2  2  - 0.0456 topo  R = 0.74 SE = 0.065 2  = 2.25 + 0.693 cos i - 0.00210 z + 0.000000571 z  2  - 0.139 topo  Where L is expressed in ( ffi[« ) > in metres MSL and cos i and topo are unitless. Topo is the ratio m  actual uniform  ct  r  z  .kyUaht illumination  toted possible uniform  skylight  R  2  •  t h  c o e f f i c i e n t  illumination  23  o f  determination and SE is the standard  Descriptive measures of these variables for the target area are as follows.  Descriptive Statistics for 329 cases STD D E V VARIABLE  MINIMUM  MAXIMUM  MEAN  STD D E V  AFTER CORRECTION  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  L r cud = L- 0.541 CO  re  0.0781 cos i + 0.000379 z - 0.000000110 z + 0.0185 topo 2  This calculation was done in floating point mode then converted to integer values. A n 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. Generally 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 atmospheric transmission. The general decrease in target radiance with elevation suggests that the effect of path radiance and sky irradiance dominates over the effect of atmospheric 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  Visual 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 appearance 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. Visually 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 accurate. 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 decreases 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  SUMMARY  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 classification 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 coefficient 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 distribution 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 correct 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  CONCLUSION  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  REFERENCES • 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) :201231. • 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 T M 81988 Godard Space Flight Centre, Greenbelt , Maryland. • Kimes, D.S., J.A.Smith & K.J.Ranson (1980), "Vegetation Reflectance Measurments 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 Distribution, U.S. Geological Survey, Arlington, VA 22202 (1979). • Lee, T.K. (1983), "Modeling and Estimating Atmosphere Effects in Landsat Imagery", MSc Thesis, Dept. of Computer Science, U.B.C., Vancouver, B.C.. • Little, J.J. (1982), "Automatic Registration of Landsat Images using Features Derived from Digital Terrain Models", Proc. Workshop on Computer Vision: Representation 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 Correction 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 Spacecraft Data", Proc. 2nd International Symposium on Remote Sensing of Environment, pp.895-934. Valley, S.L. (ed.) (1965), Handbook of Geophysics and Space Environments, McGrawHill. 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 Correction 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 L  1 + cos e  _ r  = L  0  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 byfigure10. 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  <f»* * W ^ j * ™ ™ * ™ for  but rather the ratio ground  seen  total  possible  skylight  a Lambertian surface, taking into account shadowing by adjacent terrain.  36  illumination  SOLAR IRRADIANCE A T T H E T O P OF T H E ATMOSPHERE 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  band 2 .52 - .60 10.3 13.9  band 3 .63 - .69 6.6 8.9  band 4 .76 - .90 10.8 14.6  band 5 1.55 - 1.75 3.3 4.5  Landsat MSS: Wavelength Interval (pm) Percentage of Solar Constant Solar Irradiance mW/cm 2  Landsat T M : Wavelength Interval (/xm) Percentage of Solar Constant Solar Irradiance (mW/crri ) 2  band 1 .45 - .52 10.2 13.9  based on solar constant of 135.3 mW/crri .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 m W / c m at aphelion to 139.9 mW/cm at perihelion. 2  2  3  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 integrating the standard solar spectral irradiance curve over the specified wavelength interval.  37  band 7 2.08 - 2.35 1.6 2.1  M U L T I P L E REGRESSION RESULTS F O R VARIOUS COMBINATIONS OF INDEPENDENT VARIABLES FOR BAND 5  60 60 59 58 57 57 57 57 56 56 54  SE  cos i  z  1.48 1.49 1.48 1.52 1.52 1.53 1.53 1.57 1.54 1.54 1.58  99.9 99.9 99.9  99.9 99.9 99.9 99.9 99.9 99.9 99.9 99.9 16.6 37.6 99.9  99.9 99.9  topo sky ±±£f££ z c o 8 i Confidence Level Attained (%) 99.9 99.9 99.9 99.9 99.9 99.9 99.9 92.7 99.9 99.9 99.9 99.9 99.9 99.9 99.9 99.9 99.5 99.9 98.9 99.1 84.7 99.9 z'  2  z topo  99.2 99.9  z sky  99.9  Where R 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/"« .'°''. 2  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  ORIGINAL  IMAGE  FOREST 17291 3638 13 63 21005  Image Classification Pixel Count C L E A R C U T WATER 11492 4421 7087 1534 145 259 0 943 18724 7157  FOREST 26823 5548 47 951 33369  IMAGE CORRECTED FORSUN ILLUMINATION A N G L E Image Classification Truth Map Pixel Count Pixel Count C L E A R C U T WATER ALPINE 6167 280 432 33702 FOREST 6608 103 0 12259 CLEARCUT 113 260 2 422 WATER 0 136 2112 3199 ALPINE 12888 779 5546 49582 TOTALS  ALPINE 498 0 5 2193 2696  33702 12259 422 3199 49582  Truth Map Pixel Count FOREST CLEARCUT WATER ALPINE TOTALS  IMAGE CORRECTED FORSUN ILLUMINATION ANGLE A N D ELEVATION 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 CLEARCUT 102 50 267 3 422 WATER 947 0 91 2161 3199 ALPINE 33647 12941 414 2578 49582 TOTALS I M A G E CORRECTED F O RS U N ILLUMINATION A N G L E , ELEV. A N D TOPO 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 CLEARCUT 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 IRRADIANCE M O D E L for Sun Position : azimuth 145.0° elevation 38.0° number 1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47  AZIMUTH 140 124 111 091 080 064 043 023 001 338 352 053 036 083 100 121 105 103 130 127 117 138 140 131  ELEVATION 10 12 10 11 24 14 26 20 18 17 53 51 75 49 73 63 44 22 56 42 30 42 30 21  number 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48  AZIMUTH 150 166 179 199 210 226 247 267 289 312 298 237 254 207 190 169 185 187 160 163 173 152 150 159  ELEVATION 10 12 10 11 24 14 26 20 18 17 53 51 75 49 73 63 44 22 56 42 30 42 30 21  Table 4. Point source coordinates for asymmetric sky irradiance model.  40  UNIFORM TESSELATION OF A HEMISPHERE 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36  AZIM 0 60 120 180 240 300 0 60 120 180 240 300 00 60 120 180 240 300 00 60 120 180 240 300 00 60 120 180 240 300 30 90 150 210 270 330  ELEV 10 10 10 10 10 10 28 28 28 28 28 28 45 45 45 45 45 45 60 60 60 60 60 60 75 75 75 75 75 75 64 64 64 64 64 64  37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72  AZIM 30 90 150 210 270 330 30 90 150 210 270 330 19 41 79 101 139 161 199 221 259 281 319 341 19 41 79 101 139 161 199 221 259 281 319 341  ELEV 37 37 37 37 37 37 06 06 06 06 06 06 50 50 50 50 50 50 50 50 50 50 50 50 04 04 04 04 04 04 04 04 04 04 04 04  073 074 075 076 077 078 079 080 081 082 083 084 085 086 087 088 089 090 091 092 093 094 095 096 097 098 099 100 101 102 103 104 105 106 107 108 109  AZIM 14 46 74 106 134 166 194 226 254 286 314 346 11 49 71 109 131 169 191 229 251 289 311 349 23 37 83 97 143 157 203 217 263 277 323 337 00  ELEV 35 35 35 35 35 35 35 35 35 35 35 35 18 18 18 18 18 18 18 18 18 18 18 18 22 22 22 22 22 22 22 22 22 22 22 22 90  Table 5. Position coordinates for the uniform tesselation of a hemisphere.  41  ELEVATION (m)  Figure 1. Variation of optical depth and upward and downward 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 conditions 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 conditions 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  1  < W  • • • • • • • • • • • • • • • • • • •  • • • •  ELEVATION RANGE  ^  •  ,  1  m  • • • • —  1  • • • • • * •  OF TARGET PIXELS  • • • • ••* «  • • ••  \  X. \  BAND 5  "  BAND 4  + ••  _•  BAND 6  0.1  • • • • • • •  •  • •  _  • • • • • • • • • •  • • •  0.2  -  \BAND 7  •  900  1200  y  •  1600  2000  2400  ELEVATION (m MSL)  Figure 5. Estimated variation of target radiance with elevation after correction for sun illumination angle and topo variable. Radiance values normalized to zero at 1200m M S L for all bands. 46  2800  Figure 6. Original image (false colour composite, band 4 displayed 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 maximum 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 distribution for sun position; elevation 38°, azimuth 145" based on Stewart (1984] and Steven [1977]. Contours are relative irradiance. Point source positions simulate distribution.  50  Figure 10. Uniform tesselation of a hemisphere. Centres of hexogons correspond to point source positions. Radial angle transformed 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