Next Article in Journal
Hyperspectral Classification via Superpixel Kernel Learning-Based Low Rank Representation
Next Article in Special Issue
Deep Memory Connected Neural Network for Optical Remote Sensing Image Restoration
Previous Article in Journal
Vegetation Optical Depth and Soil Moisture Retrieved from L-Band Radiometry over the Growth Cycle of a Winter Wheat
Previous Article in Special Issue
A Variational Model for Sea Image Enhancement
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Ground Reflectance Retrieval on Horizontal and Inclined Terrains Using the Software Package REFLECT

1
Department of Applied Geomatics, Université de Sherbrooke, Sherbrooke, QC J1K 2R1, Canada
2
Department of Geography, Université de Montréal, Montreal, QC H2V 2B8, Canada
3
Agriculture and Agri-Food Canada, Saint-Jean-sur-Richelieu, QC J3B 3E6, Canada
*
Author to whom correspondence should be addressed.
Remote Sens. 2018, 10(10), 1638; https://0-doi-org.brum.beds.ac.uk/10.3390/rs10101638
Submission received: 31 August 2018 / Revised: 7 October 2018 / Accepted: 12 October 2018 / Published: 15 October 2018
(This article belongs to the Special Issue Data Restoration and Denoising of Remote Sensing Data)

Abstract

:
This paper presents the software package REFLECT for the retrieval of ground reflectance from high and very-high resolution multispectral satellite images. The computation of atmospheric parameters is based on the 6S (Second Simulation of the Satellite Signal in the Solar Spectrum) routines. Aerosol optical properties are calculated using the OPAC (Optical Properties of Aerosols and Clouds) model, while aerosol optical depth is estimated using the dark target method. A new approach is proposed for adjacency effect correction. Topographic effects were also taken into account, and a new model was developed for forest canopies. Validation has shown that ground reflectance estimation with REFLECT is performed with an accuracy of approximately ±0.01 in reflectance units (for the visible, near-infrared, and mid-infrared spectral bands), even for surfaces with varying topography. The validation of the software was performed through many tests. These tests involve the correction of the effects that are associated with sensor calibration, irradiance, and viewing conditions, atmospheric conditions (aerosol optical depth AOD and water vapour), adjacency, and topographic conditions.

1. Introduction

Multispectral satellite imagery, especially at high spatial resolution (30 m or finer), represents an invaluable source of information for decision making [1] in various domains in connection with natural resources management, environment preservation, or urban planning and management. Common applications of remotely sensed imagery concern vegetation cover (species, biomass, damages, crop yields, etc.), surface conditions (land use, types of materials, soil physics, and chemistry, etc.) or the inventory of natural resources [2]. This involves the transition from a relative scale (gray scale) of measured physical variables (radiances) to a scale of biophysical variables (biomass, leaf cover, etc.) or a system of classes (e.g., type of land use). This transition is generally accomplished by combining a series of models, algorithms, equations, and/or classification methods. Depending on the spatial, spectral, and radiometric resolution of the sensors, it is possible to identify various classes and conditions of objects. The basic assumption is that the signal recorded by the sensor is directly related to the state and/or nature of the target objects on the ground.
Variation in object reflectance in the solar spectrum is the key parameter for many optical remote sensing applications. However, in addition to object reflectance values, sensor measurements are affected by many factors. Some of them are interrelated, and can be summarised as follows: (1) the atmosphere: composition and optical properties of gases and aerosols; (2) the sun: its position in the sky during data acquisition and its spectral irradiance; (3) the terrain: elevation and topography; (4) ground reflectivity: reflectance anisotropy and reflectivity of the surrounding area; and (5) the sensor: calibration, spectral sensitivity, viewing angle, and altitude. In this article, the REFLECT software package is presented, which considers all of these aspects for the radiometric correction of multispectral satellite images and ground reflectance retrieval.
There are many methods that can be used for ground reflectance retrieval [3]. These methods could be grouped into two main approaches. The first is empirical, and the second is based on the physical modelling of the signal. Empirical methods require a series of recognisable targets on the images of invariable and known reflectance. Therefore, it is therefore by regression to establish the relationship between the signal measured by the sensor and the ground reflectance. Using this relationship, the reflectance of any pixel on the images is then estimated [4,5,6]. The underlying hypothesis here is that the coefficients of the equations, which link the numerical values to the reflectance, are constants for the entire image. However, these coefficients are variable over time, and at best can be used only for one date. The physical modelling approach is by far the most widely used. The atmospheric effects (additive and multiplicative) and ground irradiances (direct and diffuse) are estimated by applying approximate solutions of the radiative transfer equation in the earth–atmosphere system. These solutions are incorporated in atmospheric codes simulating the desired parameters as a function of a series of inputs.
The purpose of this paper is to present the theoretical and practical bases of the REFLECT software package for ground reflectance retrieval based on physical modelling (6S atmospheric code). The first version of this software was developed in the remote sensing laboratory of the Department of Geography at the University of Montreal in order to estimate the reflectance of water in coastal environments based on Landsat-7 ETM+ images [7]. A second version was developed in collaboration with Agriculture and Agri-Food Canada’s Horticulture Research and Development Centre (HRDC), St-Jean-sur-Richelieu, Quebec. The aim of calculating the ground reflectance of crops is based on Landsat-7 ETM+ images [8,9]. The current version that is presented in this paper includes several new features concerning: (a) terrain with variable topography; (b) adjacency effects; (c) atmospheric conditions; and (d) observation conditions and sensor properties. In order to better present the theoretical and practical bases of REFLECT, this article is subdivided into eight sections. The “Material and Methods” are described in Section 2, Section 3 and Section 4, including a brief overview of the Second Simulation of the Satellite Signal in the Solar Spectrum (6S) code, basic elements of REFLECT, and the methodology that was used for validation. The “Results and Discussions” are presented in Section 5, Section 6 and Section 7, where validation tests of the ground reflectance retrieval are shown. These tests involve the correction of effects that are associated with the sensor calibration, irradiance, and viewing conditions, atmospheric conditions (aerosol optical depth AOD and water vapour), adjacency, and topographic conditions. Section 8 concludes the paper.

2. The 6S Atmospheric Code

The version of the 6S code that was used in REFLECT is the one published in 2006 by Vermote, E.F. et al. [10]. Conversely, 6S can estimate the reflectance of a ground object based on spectral radiance measured by a sensor under atmospheric conditions that are also defined by the user. The 6S code assumes a cloud-free atmosphere stratified in parallel layers, in which each one is characterised by its content in gas molecules and aerosol particles. The atmospheric constituents vary in space (water vapour and aerosols). A way of addressing this variability is proposed in REFLECT.
The multi-layer structure of the atmosphere in the 6S code divides the atmospheric column into 34 levels between the ground and the top of the atmosphere, allowing the elevation of the target and the altitude of the sensor to be considered. Thus, the code accounts for the reduction in the number of scattering particles (molecules and aerosols) and the absorbing gases that are above a target at a certain altitude. The atmospheric pressure and temperature profiles are used to adjust the optical depths of the atmospheric constituents, which are more concentrated in the lower layers. Standard models for the vertical distribution of atmospheric gases (US62, MIDSUM, etc.) are used therein.
Gaseous absorption is computed using the Goody model for H2O and the Malkmus model for O2, CO2, O3, N2O, and CH4. Scattering by atmospheric gases is modelled by the Rayleigh theory, while aerosol scattering is described by the Lorenz–Mie theory [11]. The single scattering hypothesis is valid solely if the optical path is very thin (τ << 1) and the scattering albedo is very low (SC << 1). Since these conditions do not always hold [12], multiple scattering is considered using the successive-orders-of-scattering method as an approximation of the radiative transfer equation solution [13].
The 6S code assumes that the atmosphere is limited by a horizontal ground surface. Three cases are then accounted for: (a) the target object is part of a Lambertian surface that extends to infinity; (b) the target object is part of an anisotropic surface (bidirectional properties) that extends to infinity; and (c) the target object occupies the central portion of a regular surface (a circle) of limited extent with Lambertian properties, and is surrounded by a Lambertian surface that extends to infinity with a level of reflectance that is different from that of the target object. Among these three cases, the latter one most closely approaches the conditions that are usually observed in remote sensing. This configuration was selected for REFLECT. However, significant modifications were introduced, as described in Section 3.2.
The code formulates the problem of radiative transfer in terms of reflectance factors, and not in terms of radiance. All of the radiances that are generated in the Earth–atmosphere system are expressed relative to the Lambertian radiance ( L L a m b s a t ) of a target with unit reflectance, which was illuminated and observed under the same conditions as the target, and placed at the top of the atmosphere. This radiance can be expressed as follows (see symbols in Abbreviations):
L L a m b s a t ( λ ¯ ) = 1 π E 0 ( λ ¯ )   d s 2 cos θ s  
Hence, the radiance Lsat, which is measured by the sensor at a spectral band designated by its central wavelength λ ¯ , is described in terms of apparent reflectance factor ρ s a t by:
ρ s a t ( λ ¯ ) = L s a t ( λ ¯ ) / L L a m b s a t ( λ ¯ )  
For the sake of readability, the references to wavelength and the viewing and illumination geometry are sometimes omitted in the equations.
Although atmospheric scattering may be considered independent from gaseous absorption, water vapour presents a problem, since it is often concentrated in the lower atmospheric layers where the aerosols are also concentrated. To take this into account, the code proposes the following formulation (the case of a Lambertian surface):
( ρ s a t ) j = 1 , 2 , 3 = ρ R T g a s O G + ρ A ( T g a s O G + H 2 O ) j = 1 , 2 , 3 + T g a s O G + H 2 O [ T 1 S a l b ρ e n v ( t d i r ρ t a r + t d i f f ρ e n v ) ]  
The atmospheric reflectance is divided into two components that were introduced by the gas molecules (first term on the right) and aerosols (second term on the right). The index OG means “Other Gases” than water vapour. Equation (3) can lead to three different formulations (index j = 1, 2, 3). The first corresponds to an extreme situation where water vapour has a minimal effect on aerosol scattering, since it is concentrated at a layer located below the aerosol layer (j = 1). In this case, we find the formulation usually cited in the literature:
ρ s a t = T g a s ( ρ a t m + T ( θ s ) t d i r ( θ v )   ρ t a r   + T ( θ s ) t d i f f ( θ v )   ρ e n v 1 S a l b ρ e n v )  
The formulation corresponding to j = 3 in Equation (3) also describes an extreme case where the water vapour has a maximal effect, since it is mainly concentrated at a layer above the aerosol layer. Finally, the case j = 2 assumes that half of the water vapour is concentrated at the same layer as the aerosols, and therefore, the water vapour has an effect that is halfway between the two extreme cases.
In Equations (3) and (4), the term T ( θ s ) t d i r ( θ v ) ρ t a r contains the useful signal ρ t a r (Ltar in Figure 1) that must be estimated from satellite measurements by eliminating the multiplicative effects of the atmosphere (total transmittance in the descending path and direct transmittance in the ascending path). The term T ( θ s ) t d i f f ( θ v ) ρ e n v is due to the radiance reflected by the environment surrounding the observed target (adjacency effect) that was introduced into the field of view of the sensor by atmospheric scattering (Lenv in Figure 1). This radiance is also considered as a noise (additive), since it does not contain information about the target viewed by the sensor; it is significant for low-reflectance targets. The term in the denominator is added in order to account for multiple reflections between the ground and the atmosphere, which is characterised by its spherical albedo, Salb. This component can be significant in case of a very hazy atmosphere and a highly reflective environment.

3. REEFLECT: Principles and Main Developments

This section presents the basic elements of the proposed software package: atmospheric properties, surface properties, topographic terrain with variable topography, adjacency effect, observation conditions, and sensor properties.

3.1. Atmospheric Properties

This section shows the modifications made to the atmospheric properties used in radiative transfer calculations of the 6S code for both gaseous absorption and aerosol scattering.

3.1.1. Gaseous Absorption

Except for water vapour and ozone, the content of gases in the atmosphere is fixed and known. It is also possible to better represent actual atmospheric conditions by using meteorological data for the location and date of acquisition [14]. Therefore we have introduced this option in REFLECT by incorporating data on ozone content variation according to latitude and month of the year, as well as the formula of Leckner [14] (Equations (5a) and (5b)) for water vapour content estimation [15].
w = 0.493   H r P s / T a  
where w is the total water vapour content along the length of the atmospheric column in g cm−2, H r is the relative humidity as a fraction of 1, T a is the ambient temperature at ground level in Kelvin, and P s is the partial pressure of water vapour in saturated air given by the equation:
P s = exp ( 26.23 5816 / T a )  

3.1.2. Aerosol Scattering

The knowledge of aerosol properties is essential in remote sensing, since their effect on solar radiation is present throughout the optical spectrum. These properties are calculated using the Mie theory based on the microphysical properties of the particles [16]. However, the chemical composition of aerosols is quite variable, and depends both on the geographic distribution of the sources and on atmospheric dynamics. Determining the proportion of the various types of aerosols at a given location and at a particular time is therefore not easy. In most cases, standard models are used that assume the presence of a mixture of typical aerosol particles (dust, water-soluble, soot, etc.) with known and stable properties. As used in 6S, typical atmospheres (continental, urban, maritime, etc.) are constructed with different mixtures of these particles.
The atmosphere includes a type of aerosol that is called water-soluble, whose particle size increases when they are in contact with water vapour molecules. This contributes to modifying the refractive index and the effective section of the particles, and hence all of the optical properties of this type of aerosol, including the AOD [17,18,19]. Therefore, for the atmosphere types that are rich in water-soluble aerosols, as continental and urban mixtures that are found in the regions most frequently targeted by remote sensing, the AOD is expected to be higher for high values of air relative humidity.
We examined the relationship between the AOD measured by the Microtops II Sun photometer (Solar Light Company, Philadelphia, PA, USA.) at two experimental sites (Saint-Jean-sur-Richelieu and St-Valentin, Quebec, Canada), and the corresponding relative humidity measured at the closest weather station (L’Acadie, Quebec, Canada). Figure 2 shows that there is a proportional relationship (r = 0.62) between the AOD at 550 nm and the relative humidity. We can assume that this relationship is due to the proportion of water-soluble aerosols that are present in the continental atmosphere whose optical depth increases with humidity. Thus, we updated the aerosol model of the 6S code (published in 1983 by the World Meteorological Organization) with a more recent model (OPAC: Optical Properties of Aerosols and Clouds), which was presented by Hess, M. et al. [18], and takes into account the effect of relative humidity on the properties of aerosols and thus on the values of the AOD.
● The OPAC model
The OPAC model [18] describes the optical and microphysical properties of 10 aerosols that are considered the most typical, as well as their mixtures that are usually found in the atmosphere. The optical properties are: the extinction coefficient (EX), the scattering coefficient (SC) (or scattering albedo), the asymmetry parameter (ASY), and the phase function (PH). The data are available for 61 wavelengths between 0.25–40 μm, and eight relative humidity values (for the water-soluble type). OPAC proposes a larger number of atmosphere models (aerosol mixtures) than 6S.
The particle types that are considered in the OPAC model are: (1) water-insoluble (INSO), which is equivalent to dust-like in 6S; (2) several mineral-type aerosols corresponding to the desert dust particles in 6S; (3) water-soluble (WASO), which is composed of various particles such as sulfates, nitrates, and other organic matter, whose size and optical properties are a function of the air relative humidity; the properties of this type of aerosol are considered stable in 6S; and (4) soot, which is identical to that in 6S, which is made up of very absorbent (black) particles composed of carbon. OPAC proposes the following mixtures: continental clean, continental average, continental polluted, urban, desert, maritime clean, maritime polluted, maritime tropical, Arctic, and Antarctic. The properties of a mixture are the averages weighted by the percentages of each aerosol in each mixture.
A comparison of the OPAC dust-type aerosols (INSO and mineral) with the 6S dust aerosol (present in the continental and urban mixtures) has shown that the optical properties (EX, SC, and ASY) of the OPAC particles are more variable than those in the 6S [20]. The WASO aerosols in OPAC are characterised by higher extinction and scattering coefficients than the same type of aerosol that are used in the 6S code (Figure 3). Also, the higher the relative humidity, the higher these coefficients are as well, since the particles are larger (and capture more photons by their larger effective section). The asymmetry parameter of WASO in OPAC is higher for short wavelengths and lower for wavelengths greater than 1.5 μm. The phase functions of INSO in OPAC and of dust in the 6S code are fairly close and very directional [20]. For most wavelengths, the phase function curves of the WASO aerosols in the 6S code are close to those of OPAC for moderate humidity percentages (Figure 4).
● Adaptation of the OPAC aerosol parameters to the 6S code
The aerosol properties (EX, SC, and ASY) of the OPAC model are defined for 61 wavelengths, 17 of which fall between 0.4–3.0 μm (0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.9, 1.0, 1.25, 1.5, 1.75, 2.0, 2.5, and 3.0). Therefore, it has been necessary to calculate by interpolation the properties of the OPAC aerosols in the 10 reference wavelengths of the 6S code (0.4, 0.488, 0.515, 0.55, 0.633, 0.694, 0.86, 1.536, 2.25, and 3.0). The same problem arises for the phase functions: in OPAC, 112 angles between –π and π are considered, while the 6S code divides the circle into 83 angles. Of the 112 OPAC angles, we therefore chose those that are the closest to the 83 angles of the 6S code, and performed interpolations according to wavelengths [20].
● Calculation of relative humidity profiles
Using the WASO aerosol from the OPAC model requires knowing the relative humidity in the layer of the atmosphere where the aerosol-scattering properties are calculated. The relative humidity profile is reconstructed using the formula of Leckner, B. [14] (Equations (5a) and (5b)), which is based on the temperature and pressure profiles according to the atmospheric model used (e.g., MIDSUM). The parameters of the water-soluble aerosols for the relative humidity values calculated at various heights are determined by means of a linear interpolation.
● Methods for estimating the AOD and choice of an approach for REFLECT
Three approaches are commonly for measuring and characterising aerosols: intensive in situ measurement campaigns; ground-based networks; and indirect estimates by remote sensing [21]. The AOD can also be estimated from the visibility measured in weather stations. Measurement campaigns allow to characterise the physical, chemical, and optical properties of aerosols [22,23]. However, since the photoradiometers that are used are very expensive, these measurements are limited to a certain number of stations. The largest ground-based network is AERONET (AErosol RObotic NETwork), which was established by NASA [24]. This network currently includes more than 450 stations that perform spectral radiometric measurements that are capable of determining the most important properties of aerosols (single scattering albedo, size distribution, phase function, asymmetry factor, AOD). Other local networks also exist, such as the Canadian AEROCAN network (which includes about 10 stations and is part of the AERONET network), the EARLINET network (European Aerosol Research Lidar Network), and the AD-Net network (Asian Dust Network). These measurement networks are of insufficient spatial density for specific earth observation applications using high and very high-resolution satellite images.
Global indirect measurements of the AOD are regularly performed from satellite observations above dark targets with low and invariable reflectance values [21]. The satellite approach is greatly appreciated, since it allows monitoring the great spatial and temporal variability of the AOD. For example, the AVHRR (Advanced Very High-Resolution Radiometer) sensor provides AOD estimates in the visible and near infrared bands with a resolution that is close to 1 km [25,26,27,28]. The POLDER (Polarisation and Directionality of the Earth Reflectance) sensor is used to measure the properties of the aerosols above the oceans and continents with a resolution of approximately 6 km [29,30]. Since 2000, the MODIS (Moderate Resolution Imaging Spectroradiometer) sensor has been used to deduce the AOD over the oceans and portions of the ground with daily coverage and for seven multispectral bands between 0.47–2.13 μm, including a blue band for which the ground reflectance is low and aerosol scattering is high [31,32,33,34,35,36]. The AOD estimated from MODIS and MISR (Multi-angle Imaging Spectroradiometer) sensors has a relative error of about 20% compared to AERONET measurements [37,38]. The AOD data obtained from the MODIS sensor are provided at a daily scale for a latitude/longitude grid of 1° [37], with a spatial resolution of 10 km. This resolution is not adequate for the correction of the spatial variability of the AOD in a HR or VHR image (pixels of 30 m or less). Models to estimate AOD at 550 nm from visibility are proposed in the literature, but the results are very approximate [39,40].
We performed a series of in situ measurements of AOD by a Microtops II Sun photometer in the vicinity of Montreal, Quebec, Canada. These measurements were compared to data from the nearest AEROCAN station located in Sherbrooke, Quebec, which was more than 100 km from the measurement sites. In addition, the data from this station were not available for all of the dates of interest. The AOD measured by the Microtops was also compared to estimates performed based on visibility measured by the weather station in St-Hubert, Quebec. For estimating the AOD based on visibility, we used the empirical formula: AOD550,VIS = 3/(VIS − 3), yielding values intermediate between the model of So [39] and Qiu [40]. The comparison showed that the data obtained from the AEROCAN network (Sherbrooke, Quebec, station), when available, are fairly close to the Microtops measurements. The AOD estimated based on visibility (AOD550,VIS) is often far from the measured values, which is probably because the visibility data are not precise measurements (e.g.,: 5 km, 15 km, 23 km…).
Given the state of the art concerning the sources of AOD data, the approach adopted in REFLECT is the dark target method applied to the images to be corrected. This approach is explained in Section 3.4.
● New AOD spectral extrapolation model
For each wavelength of the spectral bands of each sensor, the scattering parameters ( t d i r , t d i f f , t d i r , t d i f f , ρ a t m , or S a l b ) are evaluated by interpolation of their values that were calculated in the 10 reference wavelengths (400 nm, 488 nm, 515 nm, 550 nm, 633 nm, 694 nm, 860 nm, 1536 nm, 2250 nm, and 3750 nm). In the 6S code, the AOD is estimated at these 10 wavelengths by extrapolation based on the AOD at 550 nm (AOD550) according to the power law: A O D ( λ ) = A O D 550 ( λ [ n m ] / 550 ) a . The coefficient a varies between 0 and 4, with the standard value being 2. For each spectral band, these parameters are calculated by integrating the interpolated values for all of the wavelengths covered by the spectral band, taking into account the spectral sensitivity weighting of the sensor and the exoatmospheric solar spectral irradiance [41].
Bouroubi, M.Y. et al. [9] pointed out that the extrapolation law of AOD550 that was used in the 6S code was not adapted to our study area (Montérégie, Quebec), which was possibly because of a dominance of small particles. The study of this issue led us to propose another law based on measurements (at 380 nm, 870 nm, 936 nm, and 1020 nm) that were taken by the Microtops in agricultural (Montérégie, Quebec), urban (Montreal, Quebec), and forest (Laurentians, Quebec) areas. Figure 5 shows that the power law that was used in the 6S code cannot follow the measured values of AOD across wavelengths. For a = 1, the law is suitable for short wavelengths (around 400 nm), but not for long wavelengths (800 nm and more). Conversely, the curve with a = 2 is close to the measurements for long wavelengths, but deviates for short wavelengths. We therefore propose a Gaussian function (Equation (6)) that closely matches the average measurements for the five wavelengths offered by the Microtops.
A O D ( λ ) = K 1   A O D 550 exp ( K 2 ( λ 550 ) 2 )  
The values of the coefficients are K1 = 1.91 and K2 = 0.65. The type of site (rural or urban) does not affect this relationship [20].

3.2. Surface Properties

The reflectance of the environment (adjacent pixels) affects the digital value of the viewed pixel due to scattering within the viewing angle of the sensor. This secondary source of radiance becomes important if the target object is surrounded by a highly reflective object [15]. In this case, the adjacency effect can be similar in magnitude to the reflectance of the target [42]. According to several authors, the reflectance of the environment is the sum of the reflectances of all of the surfaces surrounding the target, and is weighted by a function f(r), which depends on their distance r to the target and the diffusivity of the atmosphere [42,43,44]:
ρ e n v = 0 ρ ( r ) f ( r ) d r  
The environment function f(r) is calculated by taking into account the scattering properties of gases and aerosols. This formulation was included by Tanré, D. et al. [42] in the 5S (Simulation of the Satellite Signal in the Solar Spectrum) code and by Vermote, E.F. et al. [43] in the 6S code. The latter offers an anisotropic model to better account for the scattering properties of the atmospheric constituents in a specific direction.
Calculation of the reflectance of the environment according to the rings method proposed by Vermote, E.F. et al. [45] was introduced by Cavayas, F. et al. [8] in the first version of the REFLECT program. This method assumes that the target occupies the central portion of a circular area of radius r. The tests performed with the rings method led to the conclusion that the phenomenon of adjacency should be more significant than estimated by Vermote’s formulation [45]. An approach with the individual contribution of each pixel of the environment represents more accurately the actual conditions of a scene with a heterogeneous spatial reflectivity. Also, the specular reflection of certain pixels of the environment must be accounted for. These ideas were behind the development of the method used in REFLECT, which was based on: (1) the anisotropic environment function model proposed in the 6S code [45], (2) the independent contribution of surrounding pixels to the adjacency effect when calculating the weighted average of Equation (7), and (3) the presence of both diffuse (Lambertian) and specular components in the reflectance of the environment (Figure 6); hence:
ρ e n v = 0 ρ e ( r ) f ( r )   g ( θ s , θ v , Δ φ )   d r  
where ρ e is reflectance of each surrounding pixel (individual contribution) at a distance r, and f is the anisotropic environment function. The function g , which was inspired from the Phong model [46], takes into account the specularity of each pixel (Figure 7), and it is given by:
g ( γ s p ) = k d + k s ( cos ( γ s p ) ) n  
The coefficients kd and ks are, respectively, the proportions of diffusely and specularly reflected radiances, with kd + ks = 1. The calculation window must respect a compromise between two conditions: (i) accounting for all of the pixels of the environment for which f(r) has a significant value; and (ii) computation time savings.
In the case of a nadir view (Figure 6), the specular direction γ s p is given by:
cos γ s p = cos θ s cos θ e ( i , j ) + sin θ s sin θ e ( i , j ) cos ( φ s p φ e ( i , j ) )  
The angles θs and φ s are respectively solar zenith and azimuth, and φ s p = φ s + 180 ° . For a pixel (i,j):
θ e ( i , j ) = tan 1 ( r ( i , j ) / H s a t )  
φ e ( i , j ) = arccos ( i i 2 + j 2 ) . sign ( j )  
where r(i,j) is the distance between the environment pixel and the central pixel, and Hsat is the height of the satellite. However, since Hsat >> r, it is also possible to consider a single value θ e ( i , j ) = θ v .
● Observation of the specular effect and estimation of the parameters kd, ks, and n
Figure 8 shows a near-infrared (NIR) Landsat-7 ETM+ image of the Hertel Lake on Mount Saint-Hilaire, Montérégie, Quebec (45°32′40″N; 73°09′00″W), which is in the middle of a forest. This is an example of a very dark surface surrounded by very bright surfaces. We can observe that the profiles of the digital numbers of the NIR band follow a certain slope when they are extracted in the direction parallel to the solar azimuth. The profiles of the pixels extracted in the direction orthogonal to the solar azimuth did not show this trend. This phenomenon was observed for other dark surfaces surrounded by bright surfaces and for other sensors (data not shown) suggesting the presence of a specular effect.
The parameters (kd, ks) of Equation (9) can be deduced from these observations. In Figure 8, the apparent reflectances (simplified formula) for the dark pixels 1 and 2 are calculated, considering only the effect of the bright pixels on the dark pixels by:
For point 1: ρ d a r k s a t ( 1 ) = T g a s ρ a t m + T g a s T t d i r   ρ d a r k   + T g a s T t d i f   ρ light ( k d + k s )
For point 2: ρ d a r k s a t ( 2 ) = T g a s ρ a t m + T g a s T t d i r   ρ d a r k   + T g a s T t d i f   ρ light k d
Subtracting the apparent reflectances of pixels 1 and 2 gives: k s = ρ d a r k s a t ( 1 ) ρ d a r k s a t ( 2 ) T g a s T t d i f   ρ light
By calculating the apparent reflectances for the example of Figure 9 and using the atmospheric parameters given by REFLECT, we obtain ks = 0.14 and kd = 0.86. These values are in accordance with those found by the Columbia-Utrecht Reflectance and Texture Database (CUReT) [47]. The exponent “n” is more difficult to estimate; we took the value n = 50, which is intermediate between the cases illustrated in Figure 6 where n = 20 and n = 100.

3.3. New Model for Topographic Correction

The formulation used in REFLECT is a generalisation of Equations (3) and (4) to include terrain with variable topography, as well as non-Lambertian reflectances. We explain here this formulation with expressions involving irradiances and radiances, as a first step, and then translate all of the values into reflectance terms, as used in the 6S code. The various sources of irradiance whose intensity is modulated by an inclined surface are: (1) solar irradiance transmitted directly by the atmosphere; (2) solar irradiance scattered by the atmosphere, which reaches the surface from any direction of the sky (hemispherical source); and (3) total solar irradiance (direct and diffuse), which is reflected by the surrounding area, and which illuminates the surface viewed after multiple reflections between the Earth and the atmosphere (hemispherical source). These three types of irradiance are given by the following equations:
-
Direct solar irradiance: E d i r = E 0 t d i r cos i s , in which is is the angle of incidence of the sun’s rays on an inclined terrain; it is given by: cos i s = cos θ s cos β + sin θ s sin β cos ( φ s α ) . The angles β and α are, respectively, the slope and orientation of the inclined surface.
-
Diffuse sky irradiance: E d i f f , s k y = E d i f f , s k y H o r i z g s k y , where the factor gsky takes into account the losses in intensity of the sky irradiance relative to the irradiance incident to a horizontal terrain. These losses are incurred because parts of the sky are masked by the surface itself, given its inclination.
-
Diffuse irradiance due to multiple Earth–atmosphere reflections: E d i f f , e n v = E ¯ t o t ρ ¯ e n v S a l b 1 ρ ¯ e n v S a l b g e n v ; the factor genv considers the losses because parts of the hemisphere of the sky are masked due to the inclination of the surface itself and of the surrounding area (environment).
These irradiances are reflected in the direction of the sensor by the surface, which has its own factors of bidirectional reflectance and hemispherical–directional reflectance. Hence, the radiances can be written according to their source as follows:
-
Ground radiance due to direct radiation: L d i r = ( 1 / π ) E d i r ρ B
-
Ground radiance due to diffuse sky radiation: L d i f f , s k y = ( 1 / π ) E d i f f , s k y ρ H D
-
Ground radiance due to multiple Earth–atmosphere reflections: L d i f f , e n v = ( 1 / π ) E d i f f , e n v ρ H D
Normalising these radiances in accordance with the 6S code (Equations (1) and (2)) gives:
ρ d i r = t d i r cos i s cos θ s ρ B ; ρ d i f f , s k y = t d i f f g s k y ρ H D ;   and ρ d i f f , e n v = T ¯ ρ ¯ e n v S 1 ρ ¯ e n v S g e n v ρ H D  
ρ B and ρ H D are respectively the bidirectional and hemispherical–directional components of the target reflectance ρ t a r . All of these reflectances are transmitted directly to the sensor with losses due to the atmospheric scattering depending on the transmittance t d i r . In addition, the adjacency effect can be written as a first approximation: ρ a d j = T ¯ t d i f f ρ ¯ e n v .
By homogenising all of these reflectances for the denominator 1 ρ ¯ e n v S and ignoring (in keeping with the 6S code) all of the terms involving the products of two reflectances and the spherical albedo, because their value approaches zero, we will finally get the equation:
ρ s a t = T g a s ( ρ a t m + t d i r cos i s cos θ s ρ B t d i r   + t d i f f g s k y ρ H D t d i r + T ¯ ρ ¯ e n v t d i f f   1 S a l b ρ ¯ e n v )  
Equation (12) is the REFLECT form of Equation (4), which was used in the 6S code. It allows separately modulating the direct and diffuse components of solar radiation. For the direct component, the cosine correction is applied. The definition of the correction factor gsky of the diffuse component requires the use of a diffuse irradiance on the inclined plane model. Loutzenhiser, P.G. et al. [48] presented a review of these models, some of which have already been used in remote sensing studies [49,50,51,52,53]. The one developed by Temps and Coulson [54] is the most elaborate according to Mefti, A. et al. [55]. It is given by Equation (13):
g s k y = ( 1 + cos β 2 ) × ( 1 + sin 3 β 2 ) × ( 1 + cos 2 i s sin 3 θ s )  

3.4. Dark Targets Method for AOD Estimation

The dark targets method developed in REFLECT for AOD estimation is based on clear deep water and very dense forest pixels identification by an automatic (or semi-automatic) histogram thresholding in the NIR and visible spectral bands. For clear deep-water pixels, REFLECT uses the NIR band to distinguish water from land. By keeping among the water pixels those having the lowest values in the available visible bands, turbid or shallow water pixels are avoided. The thresholds in the NIR and in the visible bands (examples: blue, green, and red) are defined, respectively, as the first local minimum (MINNIR) and the first local modes (MODVIS = MODblue, MODgreen, and MODred) of the histograms of the image digital numbers (Figure 10). However, for an image contaminated by clouds, and to avoid selecting cloud-shadowed pixels, a stricter condition can be imposed in the NIR band by defining a lower threshold than MINNIR; for example, the first mode of the histogram MODNIR.
The red–infrared relationship is used to differentiate vegetation from bare soil and water surfaces in the extraction of dense forest pixels (Figure 11). Among these pixels, only those that have an appropriate reflectance in the red band are kept. The condition to select a pixel (with the digital number DN) as dense forest is based on three thresholds: ( D N N I R D N Re d > T 1 ) and ( T 2 D N Re d < T 3 ). To find the T1 threshold, the digital numbers corresponding to reflectances of 0.02 in the red band and 0.15 in the NIR band are simulated for the conditions under which the image was acquired. The atmospheric parameters are calculated using REFLECT atmospheric routines by assuming the presence of a clear atmosphere (AOD = 0.05). The T3 threshold is the first mode of the red band histogram, which was calculated for the pixels that meet the first criterion. The T2 threshold was added to avoid selecting the border pixels of lakes and the edges of cloud-shaded areas. It was set at T2 = T3 − 3.
In order to give the user more options to choose the “acceptance” thresholds of dark targets, we added to the algorithm a manual selection of the “degree of severity” of the darkness conditions. This option enables the user, for water-type dark targets, to vary the thresholds from MINNIR to MODNIR for the NIR and some digital numbers around MODVIS for the visible bands. The same principle is applied for threshold T1 when looking for dense forest pixels.
Another improvement was made to the search process for dark targets due to the following observation: in the case of images with a heterogeneous spatial distribution of the AOD, the histogram thresholding technique selects only dark targets that are located in the parts of the image where the atmosphere is the clearest. Thus, REFLECT allows dividing the scene into NxM sub-images and conducting the search in each of them. This local search enables detecting the spatial variation of the AOD and correcting it. It should be noted that Ouaidrari, H. and Vermote, E.F. [56] also proposed dividing the image into 4 × 4 quadrants to search for the dark targets in an image with a non-homogeneous atmosphere. Ahern, F.J. et al. [57] and Teillet, P.M. et al. [58] also addressed this issue.

3.5. Sensors Integration

We have incorporated in REFLECT most of the HR and VHR sensors that are currently used for Earth observation applications. These sensors are Landsat-5 TM, Landsat-7 ETM+, Landsat-8 OLI, Terra ASTER, SPOT-1 to SPOT-7, Sentinel-2, Ikonos, QuickBird, Pleiades, and WorldView-2 and WorldView-3. New sensors will be incorporated as they come on line. Incorporating a sensor requires: (1) the calibration coefficients “gain” and “offset” for each spectral band (provided with the image as metadata) that are used to convert the digital numbers (given in eight or 16 bits) to apparent radiances in W.m-2.sr-1.μm-1 according to sensor technical documentation; and (2) the spectral sensitivities SSk(λ) used for the averaging of spectral properties ( L k s a t , Tgas, t d i r , t d i f f , t d i r , t d i f f , ρ a t m , or S a l b ) given at the wavelengths λ (with a step of Δ λ = 2.5 nm in the 6S code) within the limits of the spectral band k. Hence, for any parameter p, this integration is performed by the equation:
p ¯ k = ( λ 1 λ 2 p ( λ ) S S ( λ ) E 0 ( λ ) Δ λ ) / ( λ 1 λ 2 S S ( λ ) E 0 ( λ ) Δ λ )  
E0(λ) is the exoatmospheric solar spectral irradiance.
The calibration coefficients are determined during pre-launch operations, but their values drift over the years with the age of the optical system. Certain methods are proposed for making post-launch corrections [59,60], and these coefficients are regularly updated by the satellite operators. Since the satellites cover almost the entire planet, the areas observed often have very different reflectances; consequently, several gains are available on certain sensors to optimise the accuracy of the images (better contrast without saturation). For example, REFLECT allows the user to choose predefined calibration coefficients when it is possible (examples: ETM+ and ASTER) [61]. Otherwise, values from the image metadata file should be used.

3.6. Software Implementation

The REFLECT software was programmed in C++ language. It was provided with a graphic interface and organised in modular form to reduce the computational burden on the main program, which calculates ground reflectance. A fast computation option based on the generation of internal look-up tables of the atmospheric parameters for a combination of AOD and terrain altitude values (which distinguish the various pixels of the image) considerably reduces the computation time. Another tool is provided in the interface for the calculation of atmospheric parameters.

4. Methodology and Data Used for Validation

The accuracy of ground reflectance retrieval by REFLECT was evaluated by means of several tests and using various types of images and related data. These tests involved:
  • Calculation of atmospheric parameters: Temperature and relative humidity data for the image acquisition dates were obtained from the Environment Canada website [62] to estimate the water vapour content in the atmosphere as well as gaseous transmittances, as explained in Section 3.1. The AOD that was used for the calculation of the diffusion parameters was calculated from images using the dark targets method.
  • Correction of sensor gain effect and atmospheric effects: we used five Landsat-5 TM images, two Landsat-7 ETM+ images, one SPOT-1 HRV image, and one SPOT-5 HRG image. Images were acquired over several years between June and August for an area (intersection of the scenes of these images) of approximately 80 km × 65 km around Montreal, Quebec, Canada (Figure 12; Table 1).
    By comparing the fifth and 95th percentiles of the DNs and ground reflectance values calculated by REFLECT, we assessed the correction of sensor gain, atmospheric effects (additive and multiplicative), and conditions of illumination and observation.
  • Comparison of spectral signatures of selected materials: The ground reflectance values of three materials (water, asphalt, and vegetation) identified on the images of the scene in Figure 12 were compared with their simulated values from the USGS ASTER spectral library. The reflectance values of these materials for the spectral bands of the sensors used are simulated using following Equation (similar to Equation (14)):
    R b a n d , m a t e r i a l = ( λ 1 λ 2 ρ m a t e r i a l S S b a n d E 0 ) / ( λ 1 λ 2 S S b a n d E 0 )  
    where ρ m a t e r i a l is the spectral reflectance from the ASTER library for a given material, S S b a n d is the spectral sensitivity of the sensor by band, and E0 is the exoatmospheric spectral solar irradiance.
  • Comparison with ASD measurements: The ground reflectance values calculated for four ETM+ images acquired between 2000–2002 and one WorldView-2 image from 2011 were compared with the spectroradiometer (ASD FieldSpec HH Pro) measurements taken in agricultural fields in the Montérégie region, Quebec, Canada. The ASD reflectance values, which were measured concurrently with image acquisition, were incorporated by the spectral band of the images that was used according to the principle of Equation (16).
  • Correction of adjacency effects: A series of nine Formosat-2 images (Taiwanese experimental sensor operated by SPOT Image, a spatial resolution of 8 m, fast revisit time [63]) acquired between 5 June and 3 July 2005 near Montreal were used to show correction of the adjacency effect (in addition to correction of sensor gain and atmospheric effects) for a vegetated surface surrounded by highly reflective surfaces in the red band.
  • Topographic corrections for flat surfaces: One Ikonos image acquired on 29 August 2002 over the Paulatuk region, Northwest Territories, Canada, was used to test REFLECT’s topographic effects correction model on flat inclined surfaces (roofs of large buildings).
  • Topographic corrections for a forest canopy: Three SPOT-1 HR images acquired during the same period in 1988, 1989, and 1990 over the mountainous region of Tarn, France, as well as a DEM (digital elevation model) of the region enabled us to validate the adaptation of the topographic correction model to a forest canopy.
Validation tests (a) to (c)—where all of the surfaces are considered horizontal—are discussed together in Section 5. Tests (d) and (e) are presented in Section 6, which is devoted to vegetation targets. Tests (f) and (g) deal mainly with topographic corrections, and are presented in Section 7.

5. Reflectance Retrieval for Global Scenes

This section deals with ground reflectance retrieval for a region comprising Montreal, Quebec, Canada, which is common to nine Landsat and SPOT images (Table 1; Figure 11). This operation validates the radiometric corrections that are associated with sensor gain, illumination and observation conditions, and additive and multiplicative atmospheric effects. After an atmospheric parameter calculation step, analysis of the digital numbers and ground reflectance values was carried out as follows: (1) by comparing the fifth and 95th percentiles of the digital numbers and ground reflectance (ρground, which is ρtar of Equation (3), and includes values of all of the available images; and (2) by comparing the spectral signatures from raw (DNs) and corrected images (ρground) of selected materials from the USGS ASTER library (water, asphalt, and vegetation).

5.1. Atmospheric Parameters for TM, ETM+, HRV, and HRG Images

5.1.1. Gaseous Transmittances

Gaseous transmittances are calculated by estimating the water vapour content in the air from the air relative humidity (Hr) and ambient temperature (Ta) data obtained from Canada’s National Climate Archive for the Pierre Elliott Trudeau Airport station, which is located approximately at the centre of the scene delimited around Montreal (Table 1). The effect of Hr (lower Tgas for higher Hr, and vice versa) can be seen in the NIR band because of the water vapour absorption lines between the 800–840 nm wavelengths [20]. Gaseous transmittances are close to one in the short wavelengths (blue band), and decline with increasing proximity to the MIR (mid-infrared) wavelengths, mainly owing to water vapour absorption.

5.1.2. AOD from Dark Targets Method and Diffusion Parameters

The AOD at 550 nm is estimated using the dark targets method for all of the images (Table 2). An average value was taken because the AOD spatial distribution did not show any significant differences. We can see that most of the images have a clear atmosphere. The atmospheric parameters also show the magnitude of the additive effect in the short wavelengths (blue band). However, owing to the new interpolation rule AOD(λ) = f (AOD550, λ) given by Equation (6) of Section 3.1.2, the diffusion parameters (tdir, tdir, tdiff, tdiff, ρatm and Salb) are no longer overestimated in the blue band compared to tests carried out with the original rule of the 6S code (data not shown).

5.2. Comparison between Low and High Values of DN and Ground Reflectance

The images listed in Table 1, which were all acquired in the summer, were used to demonstrate the homogenisation of dynamic ranges by the radiometric corrections. We examined the minima (fifth percentiles) and maxima (95th percentiles) of the DNs and ground reflectance values (ρground) calculated by REFLECT. The minima provide an indication of the correction of the additive atmospheric effects and sensor offsets, while the maxima show the correction of the multiplicative atmospheric effects and sensor gain. Table 3 shows that the differences between the DNs are significant for the different dates (atmospheric conditions) and different sensors (gains).
The DNs of the SPOT-1 HRV image of 1 August 1987 are particularly low, and those of the SPOT-5 HRG image of 29 July 2006 are particularly high. The reason for this is that SPOT has variable gain settings for each spectral band. The gain values in the green, red, and NIR bands for the HRV image of 1 August 1987 were 0.852, 0.794, and 0.888, respectively, versus 2.33, 2.8, and 1.31 in the same bands for the HRG image of 29 July 2006. The fifth and 95th percentiles of ρground are much less variable; differences are around 0.01 for the low reflectance values (visible bands), and around 0.02 for the high reflectance values (NIR and MIR).

5.3. Comparison between Values of DN and Ground Reflectance for Selected Materials

Pixels corresponding to clear water (lakes), asphalt (highways), and vegetation (deciduous trees) surfaces were selected (Table 4) on all of the previously used Landsat and SPOT images. The purpose of this test was to compare the ground reflectance values calculated by REFLECT for these materials from the various images and compare these reflectance values with those obtained from the USGS ASTER library. The latter values are calculated as indicated by Equation (16).
As we can see from Table 5, there is significant variability in the DNs by type of sensor.
The variability by date is relatively less significant for the TM and ETM+ sensors. The reason for this is that the images were acquired under the same conditions of: (i) sensors gain per spectral band (unique for TM and fixed to high in visible bands and low in infrared bands for ETM+); (ii) season and time (similar illumination and observation conditions); and (iii) clear skies. However, the levels of the DNs from the SPOT HR images are not only different from those of the TM and ETM+ sensors; they also differ among themselves. As mentioned above, this is due to the difference in gains for the two acquisitions. Hence, the variation in DNs is affected more by the sensor than by the atmospheric conditions corresponding to clear skies in our case. The ground reflectance values derived from the various images show significantly lower variability compared to the DNs (Table 6).
For each material and each spectral band, the ground reflectance values are very similar for all of the sensors and all of the acquisition dates. In most cases, the variation does not exceed ±0.01 in the reflectance units. In the case of vegetation (deciduous trees), the reflectance values derived from the images are lower than those given in the ASTER library, which was probably because of mixed pixels.
Note that for the different sensors, we assigned unique values for the reflectance values of the materials from the ASTER library (Table 6). The reason for this is that the differences between the integrated reflectance values by spectral band of the sensors used (TM, ETM+, HRV, and HRG) were very small (less than 0.01).

6. Reflectance Retrieval on Vegetation Targets

6.1. Comparison with ASD Measurements for Agricultural Targets

A series of reflectance measurements were made with the ASD FieldSpec HH Pro spectroradiometer in agricultural fields (corn at several growth stages) located in the Montérégie, Quebec, Canada, concurrently with the acquisition of the Landsat-7 ETM+ images (5 June 2000, 26 July 2001, 11 August 2001, and 14 August 2002) and WorldView-2 images (5 July 2011) (Figure 15). These measurements were compared with the ground reflectance values that were obtained with REFLECT.
- Landsat-7 ETM+ images
The ASD measurements corresponding to the Landsat-7 ETM+ images showed significant variability by date and measurement point (Figure 13), indicating a certain variability in vegetation density.
The comparison of the measured reflectance values and the reflectance values that are estimated using REFLECT (Figure 14) shows that the estimates are fairly dispersed relative to the measurements with mean absolute errors (|estimate-measurement|) of 0.012 in the blue band, 0.014 in the green band, 0.019 in the red band, and 0.046 in the NIR, which corresponds to a mean relative absolute error (|estimate-measurement|/measurement) of approximately 10%. These deviations seem to be acceptable, since the images were obtained under variable atmospheric conditions and have spatial resolutions (30 m) that are very different from those of the ASD measurements (1 circular metre). In the case of the image of 5 June 2000, a systematic overestimation is observed, which was probably due to the presence of a fine cloud layer that led to a high additive effect in all of the the spectral bands (non-selective diffusion). If the AOD is increased to correct the additive effect of the blue and green bands, the NIR reflectance is increased, since the transmittances (T↓ and T↑) are reduced without increasing the atmospheric reflectance (ρatm).
- WordView-2 image
The WorldView-2 has eight multispectral bands: coastal blue (400–450 nm), blue (450–510 nm), green (510–580 nm), yellow (585–625 nm), red (630–690 nm), red-edge (705–745 nm), NIR-1 (770–895 nm), and NIR-2 (860–1040 nm, including water vapour absorption lines between 900–990 nm). Differences between the estimated reflectance from the WorldView-2 image and ASD measurements are slightly less than those found for ETM+ images (mean relative absolute error of less that 10%) (Figure 15). The good correction of the coastal blue band indicated that the AOD spectral extrapolation formulae that was proposed in REFLECT gave correct AOD values for the lower wavelengths (close to 400 nm). The original power law used in 6S gave over estimation of AOD and thus overcorrection in this band (data not shown). The reflectances that were estimated for NIR-1 and NIR-2 (particularly affected by water vapour absorption) bands were very close, indicating a good correction of water vapour effect in the NIR domain.

6.2. Formosat-2 Images of a Vegetal Surface: Atmospheric and Adjacency Effects

This test shows the atmospheric corrections—including the adjacency effect—for nine Formosat-2 images that were provided in apparent reflectance values (16-bit format). The AOD was determined on dense vegetation dark targets only, since water was too affected by specular reflection (Figure 16a). In order to show the effect of atmospheric conditions as well as the adjacency effect on agricultural applications, we selected a vegetated surface that was located in the Boucherville area, Quebec, Canada, with very bright edges in the red band (Figure 16b).
For the nine Formosat-2 images used (acquired on 5 June, 19 June, 21 June, 22 June, 25 June, 27 June, 28 June, 2 July, and 3 July 2005), we calculated the apparent NDVIs and extracted the profiles by field width (Figure 17a). Figure 16b shows that the apparent NDVIs are very different, even for dates separated by only one or two days; this is due to the difference in atmospheric conditions. The adjacency effects can also be seen on the NDVI profiles: except for the mixed pixels (containing soil and vegetation), the pixels at the edge of the field have lower NDVIs than the pixels in the centre (Figure 17b).
Figure 18 shows the red and NIR apparent reflectance values of the study area, as well as the environment reflectance values that were calculated as explained previously.
Figure 19 illustrates the NDVIs calculated from ground reflectance values retrieved by REFLECT.
Correction of the atmospheric effects reduces the NDVI underestimation obtained with apparent reflectances [64]. In addition, with correction of the adjacency effect, the edge pixels have values closer to those of the centre. By excluding the ground pixels (numbered 1 to 4 and 14 to 18) and the mixed pixels (numbered 5, 6, and 13), the comparison between the apparent and corrected NDVIs of the field pixels (7 to 12) with the centre pixels (9 and 10) shows that the NDVIs that were calculated from the corrected images are less affected by the adjacency effect than the apparent NDVI. In addition, the adjacency effect is greater for pixels 11 and 12, which were located on the side of the field where the reflectance values in the red band from the neighbouring surface are higher (northeast direction, pixels 14 to 18). Correction of the adjacency effect for pixels 11 and 12 is clearly expressed by the reduction in the difference between their NDVIs and those of the pixels in the centre of the field.

7. Reflectance Retrieval for Inclined Surfaces

We conducted two validation tests for the topographic corrections that were used in REFLECT: the first for flat inclined planes (building roofs) in a high northern latitude, and the second for a high-elevation mountainous terrain covered by forests. In the first case, the choice of high latitudes results in a high solar zenith angle and very low atmospheric effects because of the clear sky. The choice of high elevation in the second case also minimises the incidence of aerosol particles, which are fairly rare at these elevations.

7.1. Flat Inclined Surfaces on an Ikonos Image

On an Ikonos image dating from 29 August 2002 covering an area of the Paulatuk region in the Northwest Territories, Canada, we identified large buildings with inclined symmetrical roofs. Since solar elevation was very low, the effect of the difference in irradiance on the two sides of the roof is quite significant (Figure 20).
The roof slopes can be approximately estimated (ignoring the difference in azimuth between the sensor and the surface) from the ratio of the apparent lengths of each of the half roofs by means of a trigonometric calculation (Figure 21).
The ratio between the half roof with a large apparent width and the half roof with a small apparent width can be expressed as follows: A B = L + d L d = k ; therefore: d = L k 1 k + 1 .
From the view incidence angle, we obtain: tan ( i v ) = d h ; therefore: h = d tan ( i v ) .
The slope of the roof can be expressed as follows: tan ( β ) = h L = d tan ( i v ) L = L ( k 1 ) / ( k + 1 ) tan ( i v ) L
Finally: β = arctan ( ( k 1 ) ( k + 1 ) tan ( i v ) )
Angle iv represents the view incidence angle, and is calculated by taking into account the curvature of the earth; the result is iv = 28.44°. For buildings 1, 2, and 3 (Figure 19), the k values are 8/6, 7/5, and 4/2.5 pixels, which gives β slopes of 15°, 18°, and 24°, respectively.
By considering these slopes with the buildings’ orientations, the illumination (solar zenith of 60.06° and solar azimuth of 182.75°) and viewing parameters (incidence angle of 28.44°), the REFLECT radiometric corrections that were used with and without accounting for the topographic effects gives the results provided in Table 7.
From this table, we can see that there are significant differences between the DNs of the light and dark sides of the three buildings. The same is true for the ground reflectance values without topographic corrections. We can also see the effectiveness of the topographic corrections in finding reflectance values that are very close for the two sides of the same building. Hence, the proposed topographic correction method works well for flat surfaces.

7.2. Forest Canopy in Inclined Terrain on SPOT-1 Images

For the validation (and adaptation) of the topographic corrections in forested environments, we used three SPOT-1 images (5 September 1988, 24 September 1989, and 3 September 1990) of the Tarn region, France (43°39′N; 02°42′E), with a highly varied relief (elevations between 300–1200 m). The result shows that the distribution of the DNs of the NIR band for low slopes (β < 8°) is not much affected by the orientation of the pixels α relative to the solar azimuth φs. On the other hand, for slopes of more than 18°, surfaces oriented towards the sun (|αφs| <20°) have DNs that are clearly higher than those of surfaces with low slopes. However, surfaces with slopes of more than 18° and oriented in the direction opposite the sun (|αφs| > 160°) have DNs that are barely lower than the surfaces with a low slope. Hence, in a forest environment, the magnitude of the slope effect varies depending on the orientation of the surfaces relative to the sun. This observation has also been reported by Ekstrand, S. [65].
For atmospheric corrections, AOD values were estimated using the dark targets method on the lakes in the scene. Values obtained were 0.05, 0.05, and 0.07 for the 1988, 1989, and 1990 images, respectively. First, we used REFLECT to calculate the ground reflectance values without considering the topographic effects. The overestimation of the reflectance values in inclined terrains oriented toward the sun is observed. However, surfaces oriented in the direction opposite to the sun are not significantly underestimated relative to the flat terrains. This could be explained by the geotropic nature of trees and the amount of shading within the canopy that varies depending on the slope exposition to the sun [66]. Hence, by applying the topographic correction of the model that we developed for flat inclined surfaces, a very high overcorrection is obtained for the surfaces that are oriented in the direction opposite the solar azimuth. The model that was developed for inclined planes is not applicable to a forest canopy. Therefore, one would expect that the equivalent slope of a forest canopy would be gentler than that of the terrain in which it is located. We therefore propose a slope correction that depends on the orientation of the pixel relative to the solar azimuth given by:
β = κ ( f s - α ) β  
The parameters β and β’ are actual and corrected slopes, respectively. The function κ ( φ s α ) , which is less than or equal to 1, is determined by analysing the data derived from the SPOT images.
In order to find the equivalent slope that adapts our topographic correction model to the forest canopy, we looked for the multiplicative factor κ (by varying it from 1 to 0), which minimises the RMSE between the ground reflectance of inclined terrains (three slope ranges: 6° < β < 12°, 12° < β < 18°, and β > 18°) and of flat terrains (β < 6°) for all of the forested area of the three SPOT-1 images used, and all of the forest stands combined. The obtained values of κ can be smoothed based on the difference between the solar azimuth φ s and the orientation of the surface α , which was noted as δ φ , by the function :
κ ( δ φ ) = A + B ( 1 exp ( 2 π | δ φ | ) ) C  
where A = 0.1, B = 0.9, and C is an exponent that defines the shape of the curve. All that remains to determine is whether these values are of the same order of magnitude under other conditions (other images, other forested environments, etc.).
The application of topographic correction with the equivalent slopes principle yields very good results for the NIR band (the most affected by topographic effects in our example) of the SPOT-1 images used; the mean error is 0.02 reflectance units. This error is very small compared to that obtained with the Minnaert model [67], which is frequently used for this type of application [68,69], and yields an error in the order of 0.1 reflectance units. For example, we can see on the NIR band of the image of 5 September 1988 (Figure 22), in the circled areas, that the ground reflectance values of the inclined surfaces (steep slope) oriented toward the sun (low |orientation-azimuth|) have been homogenised relative to the digital numbers.

8. Conclusions

The objective of this research was to develop a tool (REFLECT software package) for ground reflectance restitution, which considers all of the radiometric effects affecting multispectral images, including illumination and viewing conditions, sensor properties, atmospheric conditions, and topography. REFLECT uses the formulation and the routines of the 6S code, and introduces several new features such as: the incorporation of a more recent aerosol model (OPAC); the proposal of a law of extrapolation of the AOD at 550 nm to the other wavelengths that are better adapted to the studied area (Quebec, Canada); accounting for the adjacency effect by means of a model inspired by the Phong model; improvement of the dark targets method that is used for estimating the AOD with the possibility of subdividing the image and the semi-automatic choice of the target “darkness” criterion; topographic corrections for smooth surfaces by means of a model that separates the direct and diffuse components of solar radiation and for forest canopy by means of the “equivalent slopes” principle; and the reduction of processing time in the case of large images by generating internal look-up tables (LUTs) of atmospheric parameters.
Validation tests showed that the ground reflectance is retrieved with a good accuracy: the variability of this reflectance for the same surfaces was approximately ±0.01 among the atmospheric conditions and sensors for all of the spectral bands. The topographic correction that was proposed in REFLECT allowed reducing the error of ground reflectance retrieval for flat inclined surfaces (roofs) observed by an Ikonos image from 0.1 without topographic corrections to 0.01. This topographic correction model adapted to forest canopy reduced the error of NIR reflectance retrieval from three SPOT images from 0.08 (without topographic correction) to less than 0.02 (with the “equivalent slopes” principle). These performances are very satisfactory compared to ATCOR (used in PCI Geomatics) or FLAASH (used in ENVI), according to the comparative study of Fuyi, T. et al. [70].

Author Contributions

Conceptualization, Y.B. and F.C.; Methodology, Y.B.; Software, Y.B.; Validation, Y.B., F.C. and N.Y.; Formal Analysis, Y.B.; Investigation, Y.B.; Resources, N.T.; Writing-Original Draft Preparation, Y.B. and W.B.; Writing-Review & Editing, W.B.; Supervision, F.C.; Funding Acquisition, N.T.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

AOD550aerosols optical depth at 550 nm
E0(λ)exoatmospheric solar spectral irradiance
Edirdirect solar irradiance
E diff, skydiffuse sky irradiance
Hrrelative humidity
Hsatheight of the satellite
Lsatapparent luminance at the satellite level
Pspartial pressure of water vapour in saturated air
Salbspherical albedo of the atmosphere
Tgastotal gas transmittance (descending and ascending path)
Taambient temperature
isangle of incidence of the solar radiation on an inclined terrain
tdir, tdirdiffusion transmittances of the direct solar radiation in the ascending and descending paths
tdiff, tdiffdiffusion transmittances of the diffuse solar radiation in the ascending and descending paths
wtotal water vapour content in the atmospheric column
αorientation of inclined surface
βslope of inclined surface
βcorrected slope (for vegetation)
λwavelength
θssolar zenith angle
θvsensor viewing angle (relative to the nadir)
ρatmatmospheric reflectance
ρsatreflectance of the target at the satellite level
ρenvreflectance of the environment
ρBbidirectional components of the target reflectance ρtar
ρHDhemispherical-directional components of the target reflectance ρtar
τtotal optical thickness of the atmosphere
φssolar azimuth angle

References

  1. Stratoulias, D.; Tolpekin, V.; de By, R.A.; Zurita-Milla, R.; Vasilios Retsios, V.; Bijker, W.; Alfi Hasan, M.; Vermote, E.A. Workflow for Automated Satellite Image Processing: From Raw VHSR Data to Object-Based Spectral Information for Smallholder Agriculture. Remote Sens. 2017, 9, 1048. [Google Scholar] [CrossRef]
  2. Dodge, R.L.; Congalton, R.G. Meeting Environmental Challenges with Remote Sensing Imagery; American Geosciences Institute: Alexandria, VA, USA, 2013; pp. 6–77. ISBN 978-0-922152-94-0. [Google Scholar]
  3. Zhu, S.; Lei, B.; Wu, Y. Retrieval of Hyperspectral Surface Reflectance Based on Machine Learning. Remote Sens. 2018, 10, 323. [Google Scholar] [CrossRef]
  4. Richardson, A.J. Relating Landsat digital count values to ground reflectance for optically thin atmospheric conditions. Appl. Opt. 1982, 21, 1457–1464. [Google Scholar] [CrossRef] [PubMed]
  5. Putsay, M.A. Simple atmospheric method for the short-wave satellite images. Int. J. Remote Sens. 2007, 13, 1549–1558. [Google Scholar] [CrossRef]
  6. Moran, M.S.; Bryant, R.; Holifield, C.D.; McElroy, S. A refined empirical line approach for retrieving surface refletance from EO-1 ALI images. Remote Sens. Environ. 2003, 78, 71–82. [Google Scholar] [CrossRef]
  7. Lavoie, A.; Cavayas, F.J.M.; Dubois, J.M. Algorithme de simulation du signal des masses d’eau côtières au niveau des capteurs satellite à haute résolution spatiale fondé sur le code atmosphérique 6S. Int. J. Remote Sens. 2001, 22, 1683–1708. [Google Scholar] [CrossRef]
  8. Cavayas, F.; Bouroubi, M.Y.; Vigneault, P.; Tremblay, N. Algorithme de correction d’images ETM+ de Landsat-7 fondé sur le code atmosphérique 6S et la méthode des cibles obscures. In Proceedings of the 25e Symposium Canadien sur la Télédétection: «De L’image à L’information», Montréal, QC, Canada, 14–17 October 2003. [Google Scholar]
  9. Bouroubi, M.Y.; Vigneault, P.; Cavayas, F.; Tremblay, N. Le Progiciel «EFLECT»pour la correction atmosphérique d’images satellites: Validation sur la Montérégie, Québec. Télédétection 2006, 6, 1–8. [Google Scholar]
  10. Vermote, E.F.; Tanré, D.J.L.; Deuzé, J.L.; Herman, M.; Morcrette, J.J. Second Simulation of the Satellite Signal in the Solar Spectrum: 6S User Guide Version 3; University of Maryland: College Park, MD, USA; Laboratoire d’optique atmosphérique CNRS: Paris, France, 2006. [Google Scholar]
  11. Chandrasekhar, S. Radiative Transfer; Dover publication Inc.: New York, NY, USA, 1960; pp. 14–327. [Google Scholar]
  12. Teillet, P.M. Rayleigh Optical Depth Comparisons from Various Sources. Appl. Opt. 1990, 29, 1897–1900. [Google Scholar] [CrossRef] [PubMed]
  13. Petty, G.W. A First Course in Atmospheric Radiation, 2nd ed.; Sundog Publishing: Madison, WI, USA, 2006; p. 459. ISBN 09729033-1-3. [Google Scholar]
  14. Leckner, B. The spectral distribution of solar radiation at the earth’s surface—Elements of a model. Solar Energy 1987, 20, 143–150. [Google Scholar] [CrossRef]
  15. Iqbal, M. An Introduction to Solar Radiation; Academic Press Inc.: Vancouver, BC, Canada, 1983; p. 408. [Google Scholar]
  16. Van De Hulst, H.C. Light Scattering by Small Particles, 1st ed.; Dover Publications Inc.: Meniola, NY, USA, 1981; p. 470. [Google Scholar]
  17. Shettle, E.P.; Fenn, R.W. Models of the Aerosols of the Lower Atmosphere and the Effects of Humidity Variations on their Optical Properties; Air Force Geophysics Lab.: Saugus, MA, USA, 1979; pp. 12–87. [Google Scholar]
  18. Hess, M.; Koepke, P.; Schult, I. Optical Properties of Aerosols and Cloud: The Software Package OPAC. B Am. Meteorol. Soc. 1998, 79, 831–844. [Google Scholar] [CrossRef]
  19. Aubé, M. Modélisation de L’évolution Spatiale et Temporelle de L’épaisseur Optique des Aérosols à L’échelle Régionale. Ph.D. Thesis, Département de Géographie et Télédétection, Faculté des Lettres et Sciences Humaines, Université de Sherbrooke, Sherbrooke, QC, Canada, 2003. [Google Scholar]
  20. Bouroubi, M.Y. REFLECT: Logiciel de Restitution des Réflectances au sol pour L’amélioration de la Qualité de L’information Extraite des Images Satellitales à haute Résolution Spatiale. Ph.D. Thesis, Département de Géographie, Faculté des Arts et des Sciences, Université de Montréal, Montréal, QC, Canada, 2009. [Google Scholar]
  21. Yu, H.; Kaufman, Y.J.; Chin, M.; Feingold, G.; Remer, L.A.; Anderson, T.L.; Balkanski, Y.; Bellouin, N.; Boucher, O.; Christopher, S.; et al. A Review of Measurement-based Assessment of Aerosol Direct Radiative Effect and Forcing. Atmos. Chem. Phys. 2006, 6, 613–666. [Google Scholar] [CrossRef]
  22. Vermeulen, A.C.; Devaux, C.; Herman, M. Retrieval of the scattering and microphysical properties of aerosols from ground-based optical measurements including polarisation. I—Method. Appl. Opt. 2000, 39, 6207–6220. [Google Scholar] [CrossRef]
  23. Dubovik, O.; King, M.D. A flexible inversion algorithm for retrieval of aerosol optical properties from sun and sky radiance measurements. J. Geophys. Res. 2000, 105, 20673–20696. [Google Scholar] [CrossRef]
  24. Holben, B.N.; Eck, T.F.; Slutsker, L. AERONET—A federated instrument network and data archive for aerosol characterization. Remote Sens. Environ. 1998, 66, 1–16. [Google Scholar] [CrossRef]
  25. Holben, B.N.; Vermote, E.; Kaufman, Y.J.; Tanré, D.; Kalb, V. Aerosol Retrieval over Land from AVHRR Data-Application for Atmospheric Correction. IEEE Trans. Geosci. Remote Sens. 1992, 30, 212–222. [Google Scholar] [CrossRef]
  26. Soufflet, V.; Tanré, D.; Royer, A.; O’Neill, N.T. Remote sensing of aerosols over boreal forest and lake water from AVHRR data. Remote Sens. Environ. 1997, 60, 22–34. [Google Scholar] [CrossRef]
  27. Ignatov, A.; Stowe, L. Aerosol retrievals from individual AVHRR channels. Part I: Retrieval algorithm and transition from Dave to 6S Radiative Transfer Model; Part II: Quality control, probability distribution functions, information contents and consistency checks of retrievals. J. Atmos. Sci. 2002, 59, 313–334, 335–362. [Google Scholar] [CrossRef]
  28. Schmechting, C.; Carrere, V.; Dubuisson, P.; Roger, J.C.; Santer, R. Sensitivity analysis for the aerosol retrieval over land for MERIS. Int. J. Remote Sens. 2003, 24, 2921–2944. [Google Scholar] [CrossRef]
  29. Kaufman, Y.J.; Tanré, T.; Remer, L.A.; Vermote, E.F.; Chu, A.; Holben, B.N. Operational remote sensing of tropospheric aerosol over land from EOS moderate resolution imaging spectroradiometer. J. Geophys. Res. 1997, 102, 17051–17067. [Google Scholar] [CrossRef] [Green Version]
  30. Kaufman, Y.J.; Boucher, O.; Tanré, D.; Chin, M.; Remer, M.L.; Takemura, T. Aerosol anthropogenic component estimated from satellite data. Geophys. Res. Lett. 2005, 32, L17804. [Google Scholar] [CrossRef]
  31. Kaufman, Y.J.; Tanré, D. Algorithm for Remote Sensing of Tropospheric Aerosol from MODIS Products; Algorithm Theoretical Basis Document, ATBD-MOD-02; NASA Goddard Space Flight Center: Greenbelt, MD, USA, 1998; 85p. [Google Scholar]
  32. Vermote, E.F.; El Saleous, N.; Justice, C.O. Atmospheric correction of MODIS data in the visible to middle infrared: First results. Remote Sens. Environ. 2002, 83, 97–111. [Google Scholar] [CrossRef]
  33. Chu, D.A.; Kaufman, Y.J.; Ichoku, C.; Remer, L.A.; Tanré, D.; Holben, B. Validation of MODIS aerosol optical depth retrieval over land. Geophys. Res. Lett. 2002, 29, MOD2-1–MOD2-4. [Google Scholar] [CrossRef]
  34. Kaufman, Y.J.; Tanré, D.; Boucher, O. A satellite view of aerosols in the climate system. Nature 2002, 419, 215–223. [Google Scholar] [CrossRef] [PubMed]
  35. Hsu, N.C.; Tsay, S.C.; King, M.D.; Herman, R. Aerosol properties over bright reflecting source regions. IEEE Trans. Geosci. Remote Sens. 2004, 42, 557–569. [Google Scholar] [CrossRef]
  36. Costa, M.J.; Silva, A.M.; Levizzani, V. Aerosol characterization and direct radiative forcing assessment over the ocean. Part I: Methodology and sensitivity analysis. J. Appl. Meteorol. 2004, 43, 1799–1817. [Google Scholar] [CrossRef]
  37. Remer, L.A.; Kaufman, Y.J.; Tanré, D.; Mattoo, S.; Chu, D.A.; Martins, J.V.; Li, L.L.; Ichoku, C.; Levy, R.C.; Kleidman, R.G.; et al. The MODIS aerosol algorithm, products, and validation. J. Atmos. Sci. 2005, 62, 947–973. [Google Scholar] [CrossRef]
  38. Tripathi, S.N.; Dey, S.; Chandel, A.; Srivastva, S.; Singh, R.P.; Holben, B. Comparison of MODIS and AERONET derived aerosol optical depth over the Ganga basin, India. Ann. Geophys. 2005, 23, 1093–1101. [Google Scholar] [CrossRef]
  39. So, C.K.; Cheng, C.M.; Tsui, K.C. Weather and Environmental Monitoring Using MODIS AOD Data in Hong Kong, China. In Proceedings of the First International Symposium on Cloud-prone & Rainy Areas Remote Sensing, Hong Kong, China, 6–8 October 2005. [Google Scholar]
  40. Qiu, J. Broadband Extinction Method to Determine Aerosol Optical Depth from Accumulated Direct Solar Radiation. J. Appl. Meteorol. 2003, 42, 1611–1625. [Google Scholar] [CrossRef]
  41. Vermote, E.F.; Tanré, D.; Deuzé, J.L.; Herman, M.; Morcrette, J.J. Second Simulation of the Satellite Signal in the Solar Spectrum: 6S User Guide Version 2; University of Maryland: College Park, MD, USA; Laboratoire d’optique atmosphérique CNRS: Paris, France, 1997. [Google Scholar]
  42. Tanré, D.; Herman, M.; Deschamps, P.Y. Influence of the background contribution upon space measurements of ground reflectance. Appl. Opt. 1981, 20, 3676–3684. [Google Scholar] [CrossRef] [PubMed]
  43. Vermote, E.F.; Tanré, D.; Deuzé, J.L.; Herman, M.; Morcette, J.J. Second Simulation of the Satellite Signal in the Solar Spectrum, 6S: An Overview. IEEE Trans. Geosci. Remote Sens. 1997, 35, 675–686. [Google Scholar] [CrossRef]
  44. Liang, S. Quantitative Remote Sensing of Land Surfaces; Wiley Series in Remote Sensing; John Wiley Kr Sons, Inc.: Hoboken, NJ, USA, 2004; pp. 76–177. ISBN 9780471723721. [Google Scholar]
  45. Vermote, E.F.; Tanré, D.; Deuzé, J.L.; Herman, M.; Morcrette, J.J. Second Simulation of the Satellite Signal in the Solar Spectrum: User Manual; University of Maryland: College Park, MD, USA; Laboratoire D’optique Atmosphérique CNRS: Paris, France, 1994. [Google Scholar]
  46. Phong, B.T. Illumination for computer generated pictures. Commun. ACM 1975, 18, 311–317. [Google Scholar] [CrossRef] [Green Version]
  47. CUReT. Available online: http://www1.cs.columbia.edu/CAVE//exclude/curet/.index.html (accessed on 20 August 2018).
  48. Loutzenhiser, P.G.; Manz, H.; Felsmann, C.; Strachan, P.A.; Maxwell, G.M. An empirical validation of modeling solar gain through a glazing unit with external and internal shading screens. Appl. Thermal Eng. 2007, 27, 528–538. [Google Scholar] [CrossRef]
  49. Cavayas, F. Modelling and correction of topographic effect using multi-temporal satellite images. Can. J. Remote Sens. 1987, 13, 49–67. [Google Scholar] [CrossRef]
  50. Sandmeier, S.; Klaus, I. A Physically-Based Model to Correct Atmospheric and Illumination Effects in Optical Satellite Data of Rugged Terrain. IEEE Trans. Geosci. Remote Sens. 1997, 35, 708–717. [Google Scholar] [CrossRef]
  51. Shepherd, J.D.; Dymond, J.R. Correcting satellite imagery for the variance of reflectance and illumination with topography. Int. J. Remote Sens. 2003, 24, 3503–3514. [Google Scholar] [CrossRef]
  52. Richter, R. Correction of atmospheric and topographic effects for high spatial resolution satellite imagery. Int. J. Remote Sens. 1997, 18, 1099–1111. [Google Scholar] [CrossRef]
  53. Richter, R. ATCOR: Atmospheric and Topographic Correction; German Aerospace Center, Mars: Oberpfaffenhofen, Germany, 2004. [Google Scholar]
  54. Temps, R.C.; Coulson, K.L. Solar radiation incident upon slopes of different orientations. Sol. Energy 1977, 19, 331–333. [Google Scholar] [CrossRef]
  55. Mefti, A.; Bouroubi, M.Y.; Adane, A. Generation of hourly solar radiation for inclined surfaces using monthly mean sunshine duration in Algeria. Energy Convers. Manag. 2003, 44, 3125–3141. [Google Scholar] [CrossRef]
  56. Ouaidrari, H.; Vermote, E.F. Operational Atmospheric Correction of Landsat TM Data. Remote Sens. Environ. 1999, 70, 4–15. [Google Scholar] [CrossRef]
  57. Ahern, F.J.; Teillet, P.M.; Goodenough, D.G. Transformation of Atmospheric and Solar Illumination Conditions on the CCRS Image Analysis System. In Proceedings of the 5th Purdue Symposium on Machine Processings of Remotely Sensed Data, West Lafayette, IN, USA, 21–23 June 1979; pp. 34–52. [Google Scholar]
  58. Teillet, P.M.; O’Neill, N.T.; Kalinauskas, A.; Sturgeon, D.; Fedosejevs, G. A Dynamic Regression Algorithm for Incorporating Atmospheric Models into Image Correction Procedures. In Proceedings of the 1987 International Geoscience and Remote Sensing Symposium (IGARSS’87), Ann Arbor, MI, USA, 18–21 May 1987; p. 913918. [Google Scholar]
  59. Teillet, P.M. A status overview of earth observation calibration/validation for terrestrial applications. Canad. J. Remote Sens. 1997, 23, 291–298. [Google Scholar] [CrossRef]
  60. Chander, G.; Markham, B.L.; Helder, D.L. Summary of current radiometric calibration coefficients for Landsat MSS, TM, ETM+, and EO-1 ALI sensors. Remote Sens. Environ. 2009, 113, 893–903. [Google Scholar] [CrossRef] [Green Version]
  61. Bouroubi, Y.; Cavayas, F.; Tremblay, N. REFLECT: Software for Ground Reflectance Restitution to Enhance the Accuracy of the Information Extracted from Satellite Images. In Proceedings of the Conference of the Canadian Remote Sensing Society, the Prairie Summit, Regina, SK, Canada, 1–5 June 2010. [Google Scholar]
  62. Canadian Weather. Available online: http://www.weatheroffice.gc.ca/canada_e.html (accessed on 21 August 2018).
  63. Hagolle, O.; Dedieu, G.; Mougenot, B.; Debaecker, V.; Duchemin, B.; Meygret, A. Correction of aerosol effects on multi-temporal images acquired with constant viewing angles: Application to Formosat-2 images. Remote Sens. Environ. 2008, 112, 1689–1701. [Google Scholar] [CrossRef] [Green Version]
  64. Song, C.; Woodcock, E.C.; Seto, K.C.; Lenney, M.P.; Scott, A.M. Classification and Change Detection Using Landsat TM Data: When and How to Correct Atmospheric Effects? Remote Sens. Environ. 2001, 75, 230–244. [Google Scholar] [CrossRef]
  65. Ekstrand, S. Landsat TM-Based Forest Damage Assessment: Correction for Topographic Effects. Photogramm. Eng. Remote Sens. 1996, 62, 151–161. [Google Scholar]
  66. Cavayas, F.; Teillet, P.M. Geometric model simulations of conifer canopy reflectance. In Proceedings of the 3rd International Colloquium on Spectral Signatures of Objects in Remote Sensing, Les Arcs, France, 16–20 December 1985; pp. 183–189. [Google Scholar]
  67. Minnaert, N. The reciprocity principle in lunar photometry. Astrophys. J. 1941, 93, 403–410. [Google Scholar] [CrossRef]
  68. Murakami, T. Minnaert constant of several forest types from SPOT/HRV data. J. Jpn. Soc. Photogramm. Remote Sens. 2002, 41, 47–55. [Google Scholar] [CrossRef]
  69. Blesius, L.; Weirich, F. The use of the Minnaert correction for land-cover classification in mountainous terrain. Int. J. Remote Sens. 2005, 26, 3831–3851. [Google Scholar] [CrossRef]
  70. Fuyi, T.; Mohammed, S.K.; Abdullah, K.; Lim, H.S.; Ishola, K.S. A comparison of Atmospheric Correction Techniques for Environmental Applications. In Proceedings of the International Conference on Space Science and Communication (IconSpace 2013), Melaka, Malaysia, 1–3 July 2013; Available online: https://0-ieeexplore-ieee-org.brum.beds.ac.uk/document/6599471 (accessed on 16 September 2013).
Figure 1. Solar radiation components received by a remote sensing sensor according to the Second Simulation of the Satellite Signal in the Solar Spectrum (6S) code.
Figure 1. Solar radiation components received by a remote sensing sensor according to the Second Simulation of the Satellite Signal in the Solar Spectrum (6S) code.
Remotesensing 10 01638 g001
Figure 2. Relation between aerosol optical depth (AOD) at 550 nm measured by Microtops-II instrument in Montérégie, Québec, Canada, and relative humidity given by the closest meteorological station.
Figure 2. Relation between aerosol optical depth (AOD) at 550 nm measured by Microtops-II instrument in Montérégie, Québec, Canada, and relative humidity given by the closest meteorological station.
Remotesensing 10 01638 g002
Figure 3. Comparison between aerosol properties extinction coefficient (EX), scattering coefficient (SC), and asymmetry parameter (ASY) of the Optical Properties of Aerosols and Clouds (OPAC) and 6S models.
Figure 3. Comparison between aerosol properties extinction coefficient (EX), scattering coefficient (SC), and asymmetry parameter (ASY) of the Optical Properties of Aerosols and Clouds (OPAC) and 6S models.
Remotesensing 10 01638 g003
Figure 4. Comparison between phase function (PH) of water-soluble aerosol in 6S (WATE) and OPAC (WASO).
Figure 4. Comparison between phase function (PH) of water-soluble aerosol in 6S (WATE) and OPAC (WASO).
Remotesensing 10 01638 g004
Figure 5. Comparison of average relative AOD to its value at 550 nm with 6S and proposed models.
Figure 5. Comparison of average relative AOD to its value at 550 nm with 6S and proposed models.
Remotesensing 10 01638 g005
Figure 6. Adjacency effect with a specular component of the reflectance of the environment.
Figure 6. Adjacency effect with a specular component of the reflectance of the environment.
Remotesensing 10 01638 g006
Figure 7. Function g ( θ s , θ v , Δ φ ) with a Lambertian and a specular component of the reflectance of the environment for different values of θs, kd, ks, and n.
Figure 7. Function g ( θ s , θ v , Δ φ ) with a Lambertian and a specular component of the reflectance of the environment for different values of θs, kd, ks, and n.
Remotesensing 10 01638 g007
Figure 8. Specular effect of the reflectance of the environment in the NIR band of a Landsat-7 ETM+ image acquired on 8 June 2001 and covering Hertel lake in Mont St-Hilaire, Quebec, Canada.
Figure 8. Specular effect of the reflectance of the environment in the NIR band of a Landsat-7 ETM+ image acquired on 8 June 2001 and covering Hertel lake in Mont St-Hilaire, Quebec, Canada.
Remotesensing 10 01638 g008
Figure 9. Specular adjacency effect of a bright pixel on a dark pixel according to the position relative to sun azimuth.
Figure 9. Specular adjacency effect of a bright pixel on a dark pixel according to the position relative to sun azimuth.
Remotesensing 10 01638 g009
Figure 10. Histogram thresholding of blue, green, red, and near-infrared (NIR) bands (e.g., image ETM+ of 8 June 2001) for clear deep-water dark targets detection.
Figure 10. Histogram thresholding of blue, green, red, and near-infrared (NIR) bands (e.g., image ETM+ of 8 June 2001) for clear deep-water dark targets detection.
Remotesensing 10 01638 g010
Figure 11. Histogram thresholding of digital number (DN) differences of the NIR and red bands (e.g., image ETM+ of 8 June 2001) for the detection of dense forest dark targets.
Figure 11. Histogram thresholding of digital number (DN) differences of the NIR and red bands (e.g., image ETM+ of 8 June 2001) for the detection of dense forest dark targets.
Remotesensing 10 01638 g011
Figure 12. Scene used to compare DNs and ground-based reflectances for five Landsat-5 TM images, two Landsat-7 ETM + images, 1 HRV SPOT image and 1 SPOT 5 HRG images.
Figure 12. Scene used to compare DNs and ground-based reflectances for five Landsat-5 TM images, two Landsat-7 ETM + images, 1 HRV SPOT image and 1 SPOT 5 HRG images.
Remotesensing 10 01638 g012
Figure 13. ASD spectroradiometer measurements in agricultural fields located in Montérégie, Quebec, Canada. The blue, green, red, and NIR spectral bands of ETM+ are indicated by blue, green, red, and purple lines, respectively.
Figure 13. ASD spectroradiometer measurements in agricultural fields located in Montérégie, Quebec, Canada. The blue, green, red, and NIR spectral bands of ETM+ are indicated by blue, green, red, and purple lines, respectively.
Remotesensing 10 01638 g013
Figure 14. Comparison of reflectances measured (by ASD) and estimated (by REFLECT) for Landsat-7 ETM+ images.
Figure 14. Comparison of reflectances measured (by ASD) and estimated (by REFLECT) for Landsat-7 ETM+ images.
Remotesensing 10 01638 g014
Figure 15. Comparison of reflectances measured (ASD) and estimated (REFLECT) for the World-View-2 image of 5 July 2011.
Figure 15. Comparison of reflectances measured (ASD) and estimated (REFLECT) for the World-View-2 image of 5 July 2011.
Remotesensing 10 01638 g015
Figure 16. Example of an agricultural field bordered by glossy surfaces in the red band on the Formosat-2 image dated 5 June 2005. Location coordinates: 45°33′30″N; 73°25′30″W; (a) global location; (b) zoom on the field.
Figure 16. Example of an agricultural field bordered by glossy surfaces in the red band on the Formosat-2 image dated 5 June 2005. Location coordinates: 45°33′30″N; 73°25′30″W; (a) global location; (b) zoom on the field.
Remotesensing 10 01638 g016
Figure 17. Profile of apparent NDVI (normalized difference vegetation index) extracted over the width of the plant surface selected on the Formosat-2 images that were used.
Figure 17. Profile of apparent NDVI (normalized difference vegetation index) extracted over the width of the plant surface selected on the Formosat-2 images that were used.
Remotesensing 10 01638 g017
Figure 18. Example of apparent reflectance and environmental reflectance in the red and NIR bands for the selected farm field on the Formosat-2 images used.
Figure 18. Example of apparent reflectance and environmental reflectance in the red and NIR bands for the selected farm field on the Formosat-2 images used.
Remotesensing 10 01638 g018
Figure 19. Corrected NDVI profile extracted on the width of the selected field on the Formosat-2 images.
Figure 19. Corrected NDVI profile extracted on the width of the selected field on the Formosat-2 images.
Remotesensing 10 01638 g019
Figure 20. Ikonos image of Paulatuk region de, Nordwest Territories, Canada, showing the illumination difference of tilted rooftops. Scene center is 69°20′30”N, 124°04′00”O; sun zenith angle is 60.06°; sun azimuth is 182.75°; viewing angle is 25.49°; viewing azimuth is 287.78°; AOD = 0.05 (very clear sky).
Figure 20. Ikonos image of Paulatuk region de, Nordwest Territories, Canada, showing the illumination difference of tilted rooftops. Scene center is 69°20′30”N, 124°04′00”O; sun zenith angle is 60.06°; sun azimuth is 182.75°; viewing angle is 25.49°; viewing azimuth is 287.78°; AOD = 0.05 (very clear sky).
Remotesensing 10 01638 g020
Figure 21. Calculation of the slope of the rooftop parts from viewing geometry of the Ikonos image.
Figure 21. Calculation of the slope of the rooftop parts from viewing geometry of the Ikonos image.
Remotesensing 10 01638 g021
Figure 22. Example of topographic effects reduction by the equivalent slopes model applied to an NIR band of the SPOT 1 image acquired on 5 September 1988 in the Tarn region, France.
Figure 22. Example of topographic effects reduction by the equivalent slopes model applied to an NIR band of the SPOT 1 image acquired on 5 September 1988 in the Tarn region, France.
Remotesensing 10 01638 g022
Table 1. Gaseous transmittances calculated for Landsat and SPOT images used for validation. Temperature and relative humidity were used for atmospheric water vapour estimation. MIR: mid-infrared.
Table 1. Gaseous transmittances calculated for Landsat and SPOT images used for validation. Temperature and relative humidity were used for atmospheric water vapour estimation. MIR: mid-infrared.
ImagesTemperature (°C)Relative Humidity (%)Total Gaseous Transmittances (Tgas)
BlueGreenRedNIRMIR
TM17 June 198422600.99310.95080.94920.91590.8723
25 July 199223650.99320.94940.94730.90810.8647
18 June 199622460.99290.95220.95120.92680.8824
27 August 199824700.99290.94630.94390.89910.8560
27 June 200527560.99330.94950.94730.90610.8636
ETM+8 June 200120420.99460.95950.95650.94700.9462
11 August 200123430.99480.95880.95510.93910.9407
HRV1 August 19872142-0.97250.95880.9339-
HRG29 July 20062760-0.96750.94780.90690.9574
Table 2. AOD values calculated with REFLECT dark objects method and corresponding atmospheric parameters (tdir, tdiff, tdir, tdiff, ρatm, and Salb) for Landsat and SPOT images used for validation.
Table 2. AOD values calculated with REFLECT dark objects method and corresponding atmospheric parameters (tdir, tdiff, tdir, tdiff, ρatm, and Salb) for Landsat and SPOT images used for validation.
BandTMTMTMTMTMETM+ETM+HRVHRG
17 June 198425 July 199218 June 199627 August 199827 June 20058 June 200111 August 20011 August 198729 July 2006
AOD5500.120.110.060.230.110.060.070.100.22
tdir
Blue0.72360.71140.76690.59770.72920.76700.7011--
Green0.81240.80350.85360.69930.81640.85380.79620.79530.7079
Red0.86970.86330.90500.77510.87260.91030.86450.86260.7971
NIR0.93270.92930.95440.87450.93420.95610.92860.93020.8902
MIR0.99810.99800.99840.99720.99820.99830.9978-0.9971
tdiff
Blue0.10090.10370.08400.13180.09960.08170.1035--
Green0.08310.08580.06400.11990.08190.06400.08820.08850.1193
Red0.06270.06480.04560.09640.06170.04350.06450.06610.0919
NIR0.02790.02870.01970.04250.02750.01930.02930.03010.0405
MIR0.00000.00000.00000.00000.00000.00000.0000-0.0000
tdir
Blue0.75840.75840.80120.67950.75840.79180.7488--
Green0.83720.83720.87620.76450.83720.87010.83060.81440.7425
Red0.88750.88750.92000.82590.88750.92060.88820.87600.8225
NIR0.94220.94220.96180.90420.94220.96130.94140.93720.9046
MIR0.99840.99840.99860.99790.99840.99850.9982-0.9975
tdiff
Blue0.09240.09240.07450.11900.09240.07520.0926--
Green0.07510.07510.05600.10510.07510.05820.07740.08260.1114
Red0.05630.05630.03960.08340.05630.03940.05610.06130.0849
NIR0.02530.02530.01720.03800.02530.01750.02580.02810.0379
MIR0.00000.00000.00000.00000.00000.00000.0000-0.0000
ρatm
Blue0.07580.07270.07000.08590.07530.07280.0770--
Green0.04260.04290.03780.05080.04230.03980.04600.04570.0535
Red0.02580.02620.02190.03350.02570.02150.02600.02460.0306
NIR0.01250.01260.00990.01830.01240.00990.01290.01050.0161
MIR0.00170.00170.00110.00300.00170.00120.0019-0.0027
Salb
Blue0.08630.08630.08980.08090.08630.09270.0887--
Green0.06260.06260.06110.06420.06260.06370.06480.06910.0686
Red0.04590.04590.04160.05170.04590.04130.04570.04920.0524
NIR0.02770.02770.02220.03610.02770.02240.02800.02900.0360
MIR0.00720.00720.00420.01220.00720.00440.0075-0.0129
Table 3. The fifth and 95th percentiles of DNs and ground reflectances for Landsat and SPOT images.
Table 3. The fifth and 95th percentiles of DNs and ground reflectances for Landsat and SPOT images.
BandTMTMTMTMTMETM+ETM+HRVHRG
17 June 198425 July 199218 June 199627 August 199827 June 20058 June 200111 August 20011 August 198729 July 2006
Fifth percentile of DN
Blue70676561666661--
Green282626232648433080
Red202019182032301650
NIR121913121317171417
MIR7108771513-16
Fifth percentile of ρground
Blue0.0390.0410.0370.0310.0340.0290.031--
Green0.0470.0450.0460.0390.0410.0330.0360.0340.039
Red0.0290.0310.0300.0270.0310.0270.0260.0270.028
NIR0.0300.0450.0320.0270.0290.0340.0360.0410.038
MIR0.0080.0110.0100.0070.0080.0110.009-0.010
95th percentile of DN
Blue118981058510511797--
Green58455237511068357205
Red65486141601239344175
NIR132133134118132141122101171
MIR13710313288131161131-185
95th percentile de ρground
Blue0.1320.1150.1210.1120.1110.1200.010--
Green0.1460.1230.1410.1180.1330.1420.1190.1210.136
Red0.1510.1240.1490.1190.1410.1530.1350.1280.143
NIR0.4940.5140.4950.5120.4870.5110.4950.4840.501
MIR0.3080.2830.3050.2750.2990.3090.295-0.301
Table 4. Sites used for comparing ground reflectances of selected materials with their values given by the USGS ASTER library. All of these sites are around Quebec, Canada.
Table 4. Sites used for comparing ground reflectances of selected materials with their values given by the USGS ASTER library. All of these sites are around Quebec, Canada.
MaterialsSiteCoordinates
Water (Lakes)Hertel lake, Mont St-Hilaire45°32′40″N; 73°09′00″W
Seigneurie lake, Mont St-Bruno45°32′50″N; 73°19′35″W
L’Achigan lake, Laurentides45°56′30″N; 73°58′00″W
Connelly lake, Laurentides45°53′50″N; 73°58′00″W
Asphalt (Roads)Highways arround Montreal
Deciduous treesSte-Thérèse-de-Blainville45°43′15″N; 73°49′00″W
Verchères45°43′16″N; 73°18′36″W
Table 5. Averages of DNs of the selected materials on Landsat and SPOT images.
Table 5. Averages of DNs of the selected materials on Landsat and SPOT images.
BandTMTMTMTMTMETM+ETM+HRVHRG
17 June 198425 July 199218 June 199627 August 199827 June 20058 June 200111 August 20011 August 198729 July 2006
DN of water (lakes)
Blue66656459636260--
Green2223212021363640101
Red151816141624252263
NIR132015141513121619
MIR897671110-21
DN of asphalt (roads)
Blue101889078919489--
Green4640413441767550151
Red4640433542737837136
NIR696969607174645968
MIR88778169758694-114
DN of vegetation (deciduous trees)
Blue73706763676863--
Green322929242854473290
Red222221192135321758
NIR1161121228712312410088118
MIR82677457758271-109
Table 6. Averages of ground reflectances of the selected materials on Landsat and SPOT images.
Table 6. Averages of ground reflectances of the selected materials on Landsat and SPOT images.
BandASTER Lib.TMTMTMTMTMETM+ETM+HRVHRG
17 June 198425 July 199218 June 199627 August 199827 June 20058 June 200111 August 20011 August 198729 July 2006
ρground of water (lakes)
Blue0.020.030.030.030.020.020.010.02--
Green0.040.030.040.030.020.020.020.010.040.05
Red0.030.020.030.020.030.030.020.010.040.03
NIR0.000.010.020.020.000.010.020.010.030.02
MIR0.000.000.010.000.000.000.000.00-0.01
ρground of asphalt (roads)
Blue0.150.100.090.090.080.080.080.08--
Green0.190.120.110.110.100.100.090.100.100.12
Red0.220.110.100.100.090.090.090.100.100.11
NIR0.260.250.260.250.250.250.250.240.250.23
MIR0.370.200.190.190.180.170.180.20-0.19
ρground of vegetation (deciduous trees)
Blue0.060.050.050.040.030.040.030.02--
Green0.090.070.070.060.050.060.060.050.050.05
Red0.060.040.040.040.030.030.040.030.030.03
NIR0.550.420.430.440.400.440.450.400.390.40
MIR0.310.180.160.170.150.170.170.16-0.18
Table 7. DNs and ground reflectances without and with topographic corrections for bright and dark sides of the rooftops in Ikonos Paulatuk.
Table 7. DNs and ground reflectances without and with topographic corrections for bright and dark sides of the rooftops in Ikonos Paulatuk.
BandsRoof 1Roof 2Roof 3
Bright SideDark SideBright SideDark SideBright SideDark Side
Orientation (°)2254522545135315
DNsBlue352340263825
Green261731183317
Red161022102511
Reflectance without topographic correctionsBlue0.1660.0650.20870.09030.1920.082
Green0.1070.0510.13770.05760.1500.051
Red0.0970.0560.13770.05600.1580.063
Reflectance with topographic correctionsBlue0.1230.1090.1520.1410.1440.145
Green0.0820.0810.1020.0930.1090.095
Red0.0760.0810.1040.0970.1180.122

Share and Cite

MDPI and ACS Style

Bouroubi, Y.; Batita, W.; Cavayas, F.; Tremblay, N. Ground Reflectance Retrieval on Horizontal and Inclined Terrains Using the Software Package REFLECT. Remote Sens. 2018, 10, 1638. https://0-doi-org.brum.beds.ac.uk/10.3390/rs10101638

AMA Style

Bouroubi Y, Batita W, Cavayas F, Tremblay N. Ground Reflectance Retrieval on Horizontal and Inclined Terrains Using the Software Package REFLECT. Remote Sensing. 2018; 10(10):1638. https://0-doi-org.brum.beds.ac.uk/10.3390/rs10101638

Chicago/Turabian Style

Bouroubi, Yacine, Wided Batita, François Cavayas, and Nicolas Tremblay. 2018. "Ground Reflectance Retrieval on Horizontal and Inclined Terrains Using the Software Package REFLECT" Remote Sensing 10, no. 10: 1638. https://0-doi-org.brum.beds.ac.uk/10.3390/rs10101638

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop