An inter-sensor calibration and atmospheric correction system for long-term time series of AVHRR imagery for coastal waters

Abstract

The Advanced Very High Resolution Radiometer (AVHRR) data series has more than 30 years of unique and valuable Earth observation imagery available. To use the long-term remote-sensor data, this article presents a system that produces consistent inter-sensor calibration and atmospheric correction for coastal waters. The system can process all five High Resolution Picture Transmission (HRPT) file formats from the twelve (12) AVHRR sensors that operated from the 1980s to the present. The system has been used to process AVHRR data of three Texas estuaries from 1985 to 2010 to document changes in suspended sediment patterns.

Keywords:

• AVHRR; suspended sediments; coastal water; atmospheric correction

1. Introduction

The amount of suspended sediment plays an important role in determining the distribution of habitats in estuaries (Montagna, Kalke, and Ritter 2002; White 2002). Observations from the Advanced Very High Resolution Radiometer (AVHRR) can document the spatial distribution of suspended sediment concentration in coastal waters over large geographic areas (Curran and Novo 1988; Stumpf 1988a, 1988b; Woodruff et al. 1999; Ruhl et al. 2001). In addition, the AVHRR data series has more than 30 years of unique and valuable Earth observation imagery acquired since 1979. Currently, other sensors, such as the Moderate Resolution Imaging Spectroradiometer (MODIS) and Sea-viewing Wide Field-of-view Sensor (SeaWiFS), are acquiring daily Earth imagery and also have been used to investigate suspended sediment concentration (Myint and Walker 2002; Miller and McKee 2004; Siegel et al. 2004; Fettweis 2007). However, these sensor sytems have relatively short-service histories compared to the AVHRR. The AVHRR data are valuable, and there is no alternative data source, in particular for years before 2000. The historical AVHRR data archive is used extensively by remote-sensing scientists (Price et al. 2004; Mennis 2006; Reed 2006; Liu et al. 2011). To take advantage of this valuable data set, a system of inter-sensor calibration and atmospheric correction was developed to process long-term AVHRR imagery time series data.

Radiometric calibration and atmospheric correction are required to obtain surface reflectance data from digital numbers (DN) in AVHRR imagery. Pertinent algorithms and documentation have appeared in various papers and software documents over an extended period due to the long history of AVHRR. This situation, however, does not facilitate comprehension of the algorithms and development of associated software. The objective of this article was to present a comprehensive procedure to perform an inter-sensor calibration and atmospheric correction procedure for NOAA AVHRR imagery of coastal waters. It includes not only a complete mathematical derivation from DN to water-leaving reflectance, but also a workflow, framework, and implementation of the computational procedure for five file formats of AVHRR data. The algorithms adopted by this article are well respected (Stumpf and Pennock 1989; Heidinger et al. 2010). Initially they appeared either in journal articles or in legacy software codes. The new system performs radiometric calibration and inter-sensor calibration at the same time following Heidinger's equations and parameter values for the inter-sensor calibration (Heidinger et al. 2010). The parameter values are the latest for AVHRR inter-sensor calibration, which were obtained by comparing archived AVHRR data with MODIS data, though Heidinger's approach is technically similar to earlier methods by Rao and Chen (1995, 1996). Although MODTRAN (Berk et al. 2006) may yield a more accurate atmospheric correction, it is complicated and patented. The algorithm used for atmospheric correction in this article covers the main components using a relatively simple approach. It is easy to comprehend, to implement, and to use on various computational platforms. To build an efficient computational system, implementation of the computational procedure focuses on streamlining computational tasks and integrating components that were developed with different languages and software tools. The comprehensive procedures provided in this article should help research scientists comprehend the remote-sensing data processing of AVHRR data for coastal waters (Miller 1993; Miller and McKee 2004).

2. Algorithms

This section includes a mathematical derivation of the four algorithms used in the system. The four algorithms are (1) inter-sensor consistent calibration including radiometric calibration, (2) atmospheric correction, (3) calculation of solar direction angles, and (4) satellite direction angles.

2.1. Radiometric calibration and inter-sensor consistent calibration

AVHRR sensors convert measured radiance to voltage and store it as a 10-bit digital count, which varies from 0 to 1023. Conversely, radiometric calibration converts the counts to radiance (Figure 1). In addition, an inter-sensor consistent calibration is required to compensate for degradation of AVHRR sensors.

Figure 1. Radiometric measurement and calibration.

An approach that performs radiometric calibration and inter-sensor calibration at the same time uses the following equations (Heidinger et al. 2010):

(1)

where R _Cal is a calibrated planet reflectance; C _dark is the dark count, which is what the instrument would measure under dark conditions; S is the calibration slope, which is also known as the inverse-gain; and C ₁₀ is the AVHRR signal in 10-bit counts. A second-order polynomial is used to calibrate the sensor degradation:

(2)

where S ₀ is the calibration slope at time t = 0 and t is the time after launch expressed in years. The values of S ₀, a, and b are provided by Heidinger et al. (2010) for AVHRR sensors on the TIROS-N to Metop-A satellites. For example, NOAA-12 has a launch time of 14 May 1991; S ₀, a, and b are 0.123, 2.624, and –0.116 and 0.147, 1.191, and –0.041 for channels 1 and 2, respectively.

R _cal must be converted to top-of-atmosphere (TOA) radiance, . For a diffuse surface, . F ₀ is the solar irradiance at the top of the atmosphere. The intensity of solar radiation received outside the Earth's atmosphere varies as the square of the Earth–Sun distance. is the correction for the Earth–Sun distance. This variation is approximated by

(3)

where E ₀ is the in-band solar irradiance constant, which depends on its equivalent band width and spectral response function. The solar irradiance is integrated over the spectral response of AVHRR channels. For NOAA-12, the integrated solar spectral irradiance of channels 1 and 2 is 200.1 and 229.9 W/m², respectively. Further, the equivalent width of the spectral response functions of channels 1 and 2 is 0.124 and 0.219 μm (Table 3.3.2-2 on website http://www.ncdc.noaa.gov/oa/pod-guide/ncdc/docs/podug/html/c3/sec3-3.htm, accessed 4 April 2013). So for NOAA-12, the solar constant is 1614 and 1050 W/(m² μm) for channels 1 and 2, respectively. Thus, we have

(4)

The is the in-band upwelling radiance in W/(m² sr μm) received by AVHRR sensors. It can also be written as L _sensor. Finally, Equation (5)(5) is used to convert a count to an in-band upwelling radiance at the top of the atmosphere for the 12 AVHRR sensors from 1985 to 2010.

(5)

2.2. Atmospheric correction

Based on analysis of legacy FORTRAN codes developed by Stumpf (personal communication) and a review of related literature (Gordon et al. 1983, 1988; Gordon and Castano 1987; Stumpf 1988a, 1988b; Stumpf and Pennock 1989), to our knowledge, this is the first complete presentation of a mathematical derivation for atmospheric correction.

Assuming that the sea surface is flat, the water reflectance (R) is determined as

(6)

where L _water is water-leaving radiance in W/(m^2. μm^.sr) from the water column, λ is the wavelength, is calculated by Equation (3)(3), and θ₀ is the solar zenith angle. The exponential term is transmission losses of the incident solar irradiance wherein is the optical thickness for Rayleigh scattering and is the ozone optical thickness. The Rayleigh attenuation is reduced by a factor of 2 because the received radiance depends not only on radiance emitted within the sensor field of view, but also on radiance scattered into the sensor solid angle from the surrounding area.

L _water can be found from the radiance received at the sensor (L _sensor) using Equation (7)(7):

(7)

where L _sensor is obtained as described in Section 2.1, and L _glint is the radiance reflected from the water surface (i.e. sun glint). Outside of regions containing sun glint, L _glint can be ignored. L _path is the atmospheric path radiance, and T ₁, as calculated in Equation (8)(8), is the atmospheric diffuse transmission coefficient from the Earth to the satellite.

(8)

where τ_aerosol is the aerosol optical thickness. L _path can be found using Equation (9)(9):

(9)

where is the Rayleigh radiance, whereas L _aerosol is the aerosol (haze) radiance. Now let L _aerosol be zero and consider only here. L _aerosol can be estimated using radiance over clear water nearby.

Thus, the water reflectance is written by

(10)

(11)

(12)

Here, is the diffuse transmittance of the atmosphere on two trips. One trip is between the sea surface and the sensor, and the other is between the sea surface and the Sun.

Under the single scattering approximation, the Rayleigh radiance is given by

(13)

The single scattering approximation divides the atmosphere into two layers: the stratosphere, where only ozone attenuation occurs with no scattering, and the troposphere, where Rayleigh and aerosol scattering occur.

is the instantaneous solar irradiance, which is equal to reduced by two trips through the ozone layer,

(14)

and are the zenith and azimuth angles, respectively, of a vector from the point on the sea surface under a pixel to the Sun.

and are the zenith and azimuth angles, respectively, of a vector from the pixel to the sensor.

(15)

is the scattering phase function of Rayleigh at wavelength . The term involving is the contribution due to photons which are backscattered from the atmosphere without interacting with the sea surface. In the term involving , term accounts for sky radiance that is specularly reflected from the sea surface into the field of view of the sensor; term accounts for photons that are first specularly reflected from the sea surface into the atmosphere and then scattered by the atmosphere into the field of view of the sensor. is the Fresnel reflectance of the interface between the atmosphere and the water body for an incident angle . It is not relative to wavelength.

and can be calculated by Equations (16) and (17), respectively:

(16)

(17)

In a summary, the water reflectance is calculated by

(18)

where and are assigned empirical constants. For example, for NOAA-12 channels 1 and 2, is set to 0.051 and 0.022, is set to 0.035 and 0.090, respectively.

2.3. Solar angles calculation

When the imaging date and time are known, the zenith and azimuth angles of the Sun can be obtained as follows:

(19)

(20)

where is the solar zenith angle, is the solar azimuth angle (calculated clockwise from the North), is the latitude of the pixel, is the declination of the Sun (related to the imaging date), and is the hour angle of the Sun (related to the imaging time). , calculated by Equation (21)(21), changes from day to day and changes with time through the day.

(21)

 = 15(12 – the current hour of the day)

The hour angle, , is the angular distance that the Earth has rotated in a day. It is equal to 15 degrees multiplied by the number of hours from local solar noon. It is positive during the morning, reduces to zero at solar noon and becomes increasingly negative as the afternoon progresses.

2.4. Sensor viewing angles calculation

The sensor viewing zenith angle, , can be calculated by the sensor viewing geometry (Niu et al. 2001):

(22)

where β is the scan angle of a pixel, a is the Earth's radius, and h is the satellite altitude. The scan angle, β, is given by

(23)

where c is the pixel's column, and V is a half (one side) of the field-of-view, which is about 55°.

Sensor viewing azimuth angle, , is calculated using Mercator projection, which is a cylindrical map projection that is conformal. Equations (24) and (25) determine the x and y coordinates of a point on a Mercator map from its latitude and longitude (with being the longitude in the center of map):

(24)

(25)

Roughly, a scan line can be assumed to have one azimuth angle. The latitudes and longitudes of the 10th and 30th Earth location points of each scan line are used to calculate points (x ₁₀, y ₁₀) and (x ₃₀, y ₃₀) on a Mercator projection plane. The viewing azimuth is first given by Equation (26)(26):

(26)

must be adjusted according to the scan direction so that it represents an angle measured clockwise from the North.

3. System implementation

The AVHRR sensor began acquiring data in 1978 with the data stored as AVHRR High Resolution Picture Transmission (HRPT) files. There are a total of five different HRPT formats of 12 AVHRR sensors from the 1980s to present (Table 1 ) (NOAA 2012a, 2012b). The new system first decodes five different HRPT formats and then transforms counts to water-leaving reflectance using the algorithms described in the previous section. Finally, the generated reflectance images are transformed onto a projected map. Figures 2 and 3 show the workflow and architecture of the system, respectively.

Figure 2. Workflow of AVHRR data processing.

Figure 3. Implementation of the AVHRR data processing system.

Table 1. Five HRPT formats of 12 AVHRR sensors from 1980s.

Sensor	HRPT format
NA = NOAA-6	* * * POD Pre-8 September 1992 * * * * *
NC = NOAA-7
NE = NOAA-8
NF = NOAA-9		* * POD between 8 September 1992 and 15 November 1994 * *	* * POD Post-15 November 1994 * *
NG = NOAA-10
NH = NOAA-11
ND = NOAA-12
NJ = NOAA-14
NK = NOAA-15				* KLM Pre-28 April 2005 *	* * Post-14 November 2006 * *
NL = NOAA-16
NM = NOAA-17
NN = NOAA-18
Note: *Signifies the sensor used the HRPT format.

The system was developed using ESRI ArcGIS using C#.NET, Python GDAL (GDAL 2011), and ArcObjects C# Add-ins. Portions of the legacy FORTRAN codes were re-written for the up-to-date computing environments and integrated with new codes to form an efficient system. C#.NET was used to implement decoding AVHRR HRPT files into ArcGIS ASCII Grid files, inter-sensors consistent calibration, calculation and interpolation of solar and satellite angles, and atmospheric correction. The inter-sensors' consistent calibration and atmospheric correction were based on the ArcGIS ASCII Grid files. C#.NET also was used to generate ASCII 3D files from the ArcGIS ASCII Grid files by combining geographic coordinates and reflectance by pixels in the format of (latitude, longitude, and reflectance). ArcGIS Python was used to import 3D points to a multipoint feature class and then interpolate those points to a raster. Python GDAL was used to combine the reflectance of channels 1 and 2 and the reflectance difference between channels 1 and 2 into a GeoTIFF file and to conduct statistical calculation on these layers. Finally, ArcObjects C# Add-in was used to populate a geodatabase of long-term water-leave reflectance for Texas bays and estuaries. This system has been used to process 25 years of AVHRR imagery from 1985 to 2010 to derive a geospatial time series of the concentration of total suspended matter of Texas bays and estuaries.

An AVHRR scene of NSS.HRPT.NM.D03364.S1717.E1727.B0788080.WI was processed by this data processing system and displayed in Figure 4. The final image is displayed using a Mercator map projection. Figure 5 shows a subset of the AVHRR scene (108 × 118 pixels) around the City of Corpus Christi, Texas. The three histograms show the change of values from the digital numbers (unitless), to the calibrated planet reflectance (unitless, in percentage), to the surface reflectance (unitless, in percentage) although the sample images of 108 × 118 pixels did not show large visual differences. The water color is different in the bay, in the gulf near barrier islands, and in the gulf far from barrier islands. The water color was used to estimate concentration of total suspended matter.

Figure 4. Raw digital number AVHRR image (top) and final projected reflectance image (bottom) of channel 1 of image NSS.HRPT.NM.D03364.S1717.E1727.B0788080.WI.

Figure 5. The Corpus Christi, Texas, area of image NSS.HRPT.NM.D03364.S1717.E1727.B0788080.WI. The area of 108 × 118 pixels are from channel 1 of this AVHRR scene. The top is the raw digital number (DN) image and its histogram. The middle is the calibrated planet reflectance image and its histogram. The bottom is the surface reflectance image and its histogram.

4. Discussion and conclusion

This article presents the algorithms of an AVHRR inter-sensors calibration and atmospheric correction data processing system for coastal waters. In general, various inter-sensor consistent calibrations for AVHRR sensors follow the same principle and similar approaches as earlier works. For example, Rao and Chen (1995, 1996) investigated the inter-sensor consistent calibration, and their approach was adopted by NOAA (NOAA 2012a, 2012b). For the atmospheric correction, Rayleigh scattering can be eliminated. Aerosol can be reduced partly via R _Diff, which is the water-leaving reflectance difference between channels 1 and 2 (Stumpf 1988a, 1988b), though the aerosol component is not considered for single bands. Since the algorithm of atmospheric correction had already produced documented results (Gordon and Castano 1987; Stumpf 1988a, 1988b; Stumpf and Pennock 1989), its applications and analysis are not included here.

The AVHRR file decoder, inter-sensors consistent calibration, and atmospheric correction represent a complete AVHRR data processing methodology. They are suitable for remote-sensing scientists to comprehend both physical theory and mathematical procedures, and how to apply the required calculations to convert radiances received by the instruments to surface reflectance. In addition, this methodology could be converted for use with Cloud Computing platforms. In a summary, this article provided a complete mathematical derivation from digital number values to water-leaving reflectance, as well as a processing system design and implementation. The comprehensive methodology should help scientists take advantage of long-term AVHHR data, in particular, the data before 2000.

Acknowledgement

This work was supported by NOAA grant “CAMEO: Building the Foundation – An Integrative Approach to Managing the Dewatering of Estuaries” (Grant Number: NA09NMF4720179).