Atmospheric Correction Guide / Conversion to Reflectance: COST, DOS, and TOA Models
For Landsat 7, 5, or 4 the best way to get surface reflectance data is to order it for free from USGS EarthExplorer; see the Update below. If you cannot order the surface reflectance data, you can follow the steps in the guide below. For Landsat 8, there are no surface reflectance products available for free download; however, the easier Landsat 8 DOS Method can be applied; learn how in the Landsat 8 Atmospheric Correction page.
Update: The USGS has recently made atmospherically corrected surface reflectance based on the LEDAPS method (pdf) from Masek et al. (2006) available for Landsat 4, 5, and 7 imagery through EarthExplorer (can be accessed in the Data Sources folder to the left; from the EarthExplorer homepage, open the Landsat CDR folder to download the surface reflectance data). The LEDAPS method is a more involved process and is the preferred surface reflectance data for Landsat 7, 5, and 4; the surface reflectance is not available for Landsat 8. This USGS page has details about the atmospheric correction method and different ways to acquire reflectance the data; (outside page; opens in new tab).
Important: the COST model has not been thoroughly assessed for accuracy for Landsat 8 which has refined bandwidths. Based on the following, DOS should be used for Landsat 8 reflectance:
Reflectance can be calculated by the models and steps included below; examples of calculated reflectances with histograms are shown in the Reflectance Examples page (can be accessed above) and can be viewed in GIS through the Landsat Crop Reflectance & Indices data download in the Landsat Custom Downloads folder. The cosine of the solar zenith correction (COST) model (Chavez, 1996; pdf) was tested and developed mainly under the parameters of 1) visible band reflectance being < 0.2, 2) near infrared (NIR) reflectance being > 0.2, and 3) the solar elevation being ≥ 45°. The COST model adds a cosine of the solar zenith factor in the denominator (TAUz) to account for atmospheric transmittance. If the above three parameters are not met, TAUz should be omitted (which is then the DOS model). The COST model has not yet been tested for accuracy for Landsat 8.
The cosine of the solar zenith angle correction (COST) model (Chavez, 1996) can be written as:
(Components of equation are described below; "CO" is from cosine; "S" is from solar; the "T" is from "TAUz".)
It is important to realize that when converting digital number (DN) to radiance (Lsatrad) that there are different gain states (high or low) for Landsat 7; the different states are a function of land surface and/or solar elevation.
The dark-object subtraction (DOS) model (Chavez, 1988) is the above COST model with TAUz omitted and is the correct equation for certain situations; for example, this can be used when solar elevation is < 45° and/or the surface reflectance is not within vegetation and soil ranges.
The top of atmosphere (TOA) model (also known as apparent reflectance or in-band planetary albedo) is the above COST model with the deduction for scatter and TAUz omitted and is the correct equation for certain situations; for example, this can be used when solar elevation is < 45° and/or the surface reflectance is not within vegetation and soil ranges and if the band is mid-infrared or larger because scatter is insignificant enough.
A precipitable water factor (Wu et al., 2005) described below, needs to be applied to band 4 (NIR) for Landsat 4, 5, and 7 (should not be applied to Landsat 8 due to refined bandwidth) to account for the absorption effect of total atmospheric precipitable water. It is important to note that the precipitable water factor was developed based on Quickbird imagery which is different from Landsat in many ways (including much finer pixel resolution) and that the factor has not been tested for accuracy with Landsat data.
For the above equation:
ρλ = Spectral reflectance of the surface.
Lsatrad = Radiance at sensor.
Lhaze1%rad = Haze radiance (path radiance or scatter) at sensor; haze is determined and is lessened by one percent (see the one percent equation below) so in theory the darkest object has one percent reflectance.
π (pi) = 3.14159265
d = Earth-sun distance in AUs calculated from the following equation from (ESA, 2007): d = (1-0.016729 cos[0.9856(DOY-4)]); where DOY = day of year. Values for d for all days of the year that are virtually the same as those derived by this equation can be found here (pdf). Landsat 8 .MTL files include this value. This value is squared in conversion to reflectance processes.
Eoλ (also referred to as ESUNλ) = Mean exoatmospheric spectral solar irradiance on a surface perpendicular to the sun’s ray (W x mˉ² x μmˉ¹). Eoλ (also know as ESUNλ) values for Landsat 4, 5, and 7 can be found here (pdf; opens in new tab). Landsat 8 ESUN values used by GIS Ag Maps are described here.
cosθs = Cosine of solar zenith angle ([solar zenith angle = 90 - solar elevation]; solar elevation for scene center is found in .MTL file; there may be online tool to acquire more local solar elevation value)
TAUv = Atmospheric transmittance along the path from the ground surface to the sensor (1.0 for Landsat).
TAUz = The transmittance factor determined by Chavez (1996) to account for atmospheric transmittance along the path from the sun to the ground surface; the TAUz value = cosθs (Chavez  also developed default values that can be used).
One-percent (1%) reflectance equation (ARSC, 2002). The one-percent equation is applied in order that the established darkest object has one-percent reflectance because of " the fact that very few targets on the Earth's surface are absolute black, so an assumed one-percent minimum reflectance is better than zero percent" (Chavez, 1996). For COST, the one-percent equation = 0.01 x ([Eoλ x cosθs²] / [d² x pi]); for DOS (TAUz is omitted), the one-percent equation = 0.01 x ([Eoλ x cosθs] / [d² x pi]). When radiance based on the one-percent equations is deducted from the established haze value (dark object subtraction value), the resulting reflectance of the established dark pixel will be one percent (.01 in reflectance units or 1% if converted to percent by multiplying by 100) instead of zero when the atmospheric is complete.
NIR total atmospheric precipitable water correction factor (Wu et al. 2005); not necessary for Landsat 8
It is important to note that the precipitable water factor was developed based on Quickbird imagery which is different from Landsat in many ways (including much finer pixel resolution) and that the factor has not been tested for accuracy with Landsat data.
The COST model was tested in a semi-arid to arid conditions (Chavez, 1996). It is important to realize that a meaningful amount of Landsat 5 and 7 band 4 (NIR) radiation can be absorbed by atmospheric precipitable water (mainly a function of relative humidity, temperature, and solar elevation) and that this affects the accuracy of retrieved reflectance (causes reflectance values to be too low). Wu et al. (2005) found that retrieved NIR reflectance was over 20 percent too low without the factor and "acceptable for agricultural applications" with it. (This should be less of a factor with Landsat 8 as the bandwidth for the NIR has been modified.) Wu et al. (2005) developed a factor based on Quickbird NIR data that helps correct for the absorption effect of atmospheric water on NIR radiation; the factor significantly improves the accuracy of retrieve NIR reflectance as is shown below. The graphic below from Wu et al. (2005) shows the improved accuracy of retrieved NIR reflectance when the precipitable water factor is applied (circles in graph).
For GIS Ag Maps atmospheric correction of Landsat 5 and 7, the correction factor replaces TAUz (cosine of solar zenith angle), lessening the denominator and increasing reflectance. The first step in determining the correction factor is to estimate atmospheric precipitable water. Wu et al. provide a chart that can be used to derive an estimation based on relative humidity and temperature as is shown below.
After the precipitable water amount is determined from the chart above it is applied along with the solar zenith to determine the factor from the contour map below. The value from the contour map is used instead of the TAUz (cosine of solar zenith value) in the denominator of the COST equation.
Haze (also referred to as path radiance) Selection Methods
Landsat 8 Dark Object Subtraction based on visible band dark object; Example 1
The following graphic shows an area that has low values in the red, green, and blue bands (center-right); the dots are centered on corresponding pixels, the colors represent the corresponding bands.
Zoomed in to dark object area you can see there are many low values and that they correspond to a shadow in an area of vegetation.
The low value pixel area can be viewed with higher resolution imagery below (the LandsatLook imagery above is based on Landsat 30 x 30 meter pixels). The low value pixels above correspond to the center area of the image below (extent of the images are similar but not exactly the same). It can be seen by viewing both images that the low values in the visible bands at the time of imaging correspond to green vegetation which is a high absorber of visible radiation (green vegetation reflects more green light than blue and red, but green reflectance is still low). The already low visible band reflectances are lowered by the shadows from clouds, and possibly from trees and/or topography. For these reasons, it makes sense that this is an area of minimal visible light reflectance. NIR reflectance may not be the lowest here because of the associated vegetation; NIR can transmit through shadows and reflect the surface.
From the cluster of points in the dark object area above, the lowest DN in each band should be selected as the haze value (unless there is an obvious erroneous value) and converted to radiance (lowest values do not need to be from the same pixel but in the same dark object area). There should be a power relationship between band haze (path) radiance (Chavez, 1988). Longer wavelength will have relatively less haze. Clearer atmospheres will have a greater path radiance difference between bands. Hazier atmospheres will have more path radiance for all bands and the haze amounts will be more similar between bands. The haze radiance below shows that a strong power relationship exists. Haze deduction does not have to be applied to bands longer that NIR.
Band Band Center Haze DN Haze Radiance
2 (blue) .480 8289 41.88
3 (green) .560 6993 23.23
4 (red) .655 6140 11.26
The band center and haze radiance values from above are plotted below with a power line.
If it is difficult to find pixels in dark object areas with your software, keep in mind that the values shown below (in the tables after the next paragraph) are very low. To find haze DNs, simply fit low values from the blue, green, and red bands to a power line with a high correlation, avoiding low values that have a large gap between them and the rest of the data (as the lowest values in the blue and green data below). This gap can commonly exist somewhere at the very low end of the histogram, it does not always exist between the lowest value and the rest of the data, but could be between a very low value and the rest of the data. Because values are low, make sure that the DNs you select for different band haze amounts correspond to radiances greater than zero; there may be some digital number near the low end of the histogram that correspond to negative radiance values. Negative radiance mainly applies to the NIR band however; is it unlikely for visible bands.
Landsat 8 low end of histogram data that show selected haze DNs; from left to right DNs are, band 2 (blue), 8289; band 3 (green), 6993; and band 4 (red), 6140. Landsat 8 data is based on a 12-bit dynamic range but are delivered as 16-bit images; previous Landsats were based on 8-bit for both. Landsat 8 DNs have a maximum value of 65,535 as opposed to previous Landsat which had a maximum value of 255. The number in the left column of each table is the ranking of the DN value (1 is the lowest DN in the entire scene); the number in the middle column is the DN value; and the number in the right column is the amount of the DN values there are in the entire scene. There are about 40 million pixels in a Landsat image.
If the histogram for band 5 (NIR) shows DN values that correspond to zero or negative radiance (there can be times where NIR or larger bands have some pixels that correspond to zero or negative radiance) such as it does in this example (5000 equals zero radiance), there is no need for haze deduction. Sometimes NIR should have a haze deduction and sometimes it should not.
Band 5 low end of histogram (5000 equals zero radiance)
As previously mentioned, bands with wavelength longer than NIR do not need haze deducted.
Landsat 8 Dark Object Subtraction based on visible band dark object; Example 2
The following describes the method used to establish haze (path) radiance; this example differs with Example 1 in that NIR has haze removed. As was the case in the Example 1, an area that has very dark visible band values is used for dark object subtraction (DOS) amounts. As with Example 1, the dark object area corresponded to vegetation under a cloud shadow; the LandsatLook imagery with the cloud shadow is not included in this example, but the corresponding high resolution imagery is included below. The graphic below shows the tree stand (in center; without the shadow that existed at time of imagery) where the low visible band values were located.
(Google Earth imagery)
Landsat 8 low end of histogram data that show selected haze DNs that are all from the tree stand are shown below. From left to right DNs are, band 2 (blue), 8075; band 3 (green), 6914; and band 4 (red), 6136. Landsat 8 data is based on a 12-bit dynamic range but are delivered as 16-bit images; previous Landsats were based on 8-bit for both. Avoid using a low value if there is a large gap somewhere between it and the remainder of the histogram. Landsat 8 DNs have a maximum value of 65,535 as opposed to previous Landsat which had a maximum value of 255. The number in the left column of each table is the ranking of the DN value (1 is the lowest DN in the entire scene); the number in the middle column is the DN value; and the number in the right column is the amount of the DN values there are in the entire scene. There are about 40 million pixels in a Landsat image.
Values for the band 5 (NIR) low end of the histogram are shown below. The NIR imagery for this scene should have a haze deduction (sometimes NIR should have haze deducted and sometimes it should not). Notice there is a large break between the lowest value and the second lowest (as there are breaks for visible bands above); the lowest value should not be used as the haze or dark object subtraction value. The DN value 5000 is zero reflectance (the maximum value is 65,535).
As previously mentioned, there should be a power relationship between haze values (Chavez, 1988). To establish an NIR haze radiance, a DN value was fit into a power line with the blue, green, and red values and was plotted against the band center wavelength for the band 5 value that derived the highest correlation. The value used is 5609. The difference in resulting reflectance units is minimal whether you select 5609 or a smaller value; reflectance will be proportionally more minimal if applying to vegetation which typically has NIR reflectance > .50. As previously explained on this page, haze should have one-percent deducted (with the 1% equation shown above); the lowest value in the histogram above that corresponds to radiance more than the radiance derived from the 1% equation is 5469 for both COST and DOS. If 5469 is used as the haze DN to calculate Lhaze1%rad, reflectance for COST and DOS would only increase by .0034 and .0031 units, respectively; this amounts to less than 1% of .50 reflectance units, which vegetation commonly exceeds. Establishing a haze value by using the power line method is less subjective and more systematic than by establishing a value based on looking at values in the low end of the histogram and focuses on keeping the power line relationship correlation high. Below is the plot of haze values for the visible bands and the NIR value that has been fit into the power line.
Dark Object Subtraction (DOS) based on histogram method
The Lhaze (path radiance or scatter) value can be based on the dark object subtraction (DOS) method where the assumption is made that within an image there are some pixels that should not have any radiance and the radiance received at the sensor is a function of atmospheric scatter (Chavez ). This assumption is then combined with "the fact that very few targets on the Earth's surface are absolute black, so an assumed one-percent minimum reflectance is better than zero percent" (Chavez, 1996). To apply the one-percent method, the scatter amount is selected then converted to radiance, then one-percent is deducted from the radiance amount. Chavez (1988) did not select the lowest value in the image but the value at the base of the low end of the histogram as is shown below (where the histogram abruptly increases, which is near the lowest value). Important: for Landsat 5 and 7 images downloaded from GloVis there will be erroneous low (nonzero) values along and near the dentil eastern and western edges of images; these should not be factored in the decision-making process for a haze value. These areas can be excluded different ways in GIS so proper haze value can be selected.
Landsat TM band 1; DN 40 (frequency of 18) is selected for the scatter value (Chavez, 1988).
Landsat TM band 2; DN 13 (frequency of 20) is selected for the scatter value (Chavez, 1988).
Landsat TM band 3; DN 12 (frequency of 69) is selected for the scatter value (Chavez, 1988).
Landsat TM band 4; DN 8 (frequency of 5) is selected for the scatter value (Chavez, 1988).
Scatter for bandwidths longer than NIR (band 4 for Landsats 5 and 7) are insignificant enough to ignore. Overall when there is more haze, the scatter amount becomes more proportionally equal between bands but the shorter wavelength band always have more scatter.
Other haze selection methods
In addition to the scatter amount being determined by viewing the histogram (Chavez, 1988), haze has also been determined by using a value corresponding to clear and deep water (Wu et al., 2005), or the lowest value in the histogram with at least 1,000 pixels (Teillet, 1995).
To calculate the correct relative scatter for Landsat 5 or 7 band 2 (green) and band 4 (NIR) based on band 3 (red) scatter, equations that follow from Hollinger (2009) can be applied; the equations were suggested by Chavez (1988) and the data to develop the equation are from Chavez (1988). Equation 1 calculates the relative scatter in band 4 based on band 3 and Equation 2 calculates the relative scatter in band 2 based on band 3.
Landsat atmospheric correction background and accuracy
Landsat imagery is free, has a relatively long historical archive, and a resolution (pixel size) of 30 meters (4.5 pixels per acre) for all bands except thermal. Currently, there are two operational Landsat satellite (Landsat 8 and Landsat 7 ETM+); Landsat 5 TM stopped providing imagery in the latter part of 2011. Landsat 8 was launched in February of 2013. There are details about Landsat in the Imaging / Landsat folder. Landsat pixels have a calibrated digital number (DN) value that, depending on the application, it may need to be converted to reflectance. Atmospheric correction derives surface reflectance by converting DNs to a ratio of the amount of radiance reflected from a surface to the total incoming amount so comparisons can be made over time. For example, it accounts for the fact that, with all else being equal, 1) DNs will be lower for the same surface when the solar elevation is lower, 2) DNs will be higher for the same surface when there are more particles in the atmosphere that scatter light back to the sensor before it reaches the ground, and 3) the slightly varying distance between the earth and the sun. Reflectance data for band 2 (green), band 3 (red), and NIR (band 4) are included in the Landsat dataset that can download and viewed in free Quantum GIS.
The COST model (Chavez, 1996) was designed and tested with Landsat 5 for soil and vegetation reflectance accuracy and should only be applied for these reflectance ranges for Landsat 5 and 7 (have virtually the same bandwidths) unless research shows otherwise. The range of values tested were generally < 0.2 for visible bands and > 0.2 for the NIR band (Moran et al., 1992). The COST model includes a factor to account for atmospheric transmittance along the path from the sun to the ground surface that the Chavez (1988) DOS model did not have; this improved accuracy for the data tested. The accuracy listed in Chavez (1996) for vegetation for bands 1, 2, 3, and 4 is 0.0107, 0.0043, 0.0127, and 0.0213, respectively.
As previously mentioned, the COST model includes a factor that accounts for atmospheric transmittance along the path from the sun to the ground surface (TAUz) which improved the accuracy compared to Chavez' (1988) DOS model (Chavez ). The following excerpt from Chavez (1996) describes the factor:
Atmospheric correction of vegetation should use the cosine of the solar zenith angle for TAUz if in the range of 0.70 (solar zenith of approximately 45°) to 0.91 (or higher), as is shown above excerpt from (Chavez, 1996). TAUz values that are less have not been tested and and should not applied.
ARSC. 2002. Arizona Remote Sensing Center: Landsat 5 atmospheric and radiometric correction. Information on website adapted from Skirvin, S (2000). Cited at: http://arsc.arid.arizona.edu/resources/image_processing/landsat/ls5-atmo.html. Last accessed: July, 2011.
Chavez, P.S., Jr. 1996. Image-based atmospheric corrections–revisited and improved. Photogrammetric Engineering and Remote Sensing 62(9): pp.1025-1036.
Chavez, P.S., Jr. 1988. An improved dark-object subtraction technique for atmospheric scattering correction of multispectral data. Remote Sensing of Environment 24: pp.459-479.
ESA. 2007. European Space Agency. Earth Observation Quality Control: Landsat frequently asked questions. Copyright 2000-2007. Cited at: http://earth.esa.int/pub/ESA_DOC/landsat_FAQ/. Last accessed: July, 2011.
Hollinger, D. 2009. A GIS-based method to predict county corn yield based on retrieved Landsat reflectance variability in western Ohio. Papers of the Applied Geography Conference 32: pp.281-290.
Masek, J. G., Vermote, E. F., Saleous, N. E., Wolfe, R., Hall, F. G., Huemmrich, K. F., Gao, F., Kutler, J., and T.K. Lim. 2006. A Landsat surface reflectance data set for North America, 1990-2000. IEEE Geosci. Rem. Sens. Lett., vol. 3, no.1, pp. 69-72, Jan. 2006.
Moran, M.S., R.D. Jackson, P.N. Slater, and P.M. Teillet. 1992. Evaluation of simplified procedures for retrieval of land surface reflectance factors from satellite sensor output. Remote Sensing of Environment 41: pp.169-184.
Teillet, P. M., and G. Fedosejevs. 1995. On the dark target approach to atmospheric correction of remotely sensed data. Canadian Journal of Remote Sensing 21: pp. 373–387.
Wu, J., Wang, D., and M.E. Bauer. 2005. Image-based atmospheric correction of Quickbird imagery of Minnesota cropland. Remote Sensing of Environment 99: pp.315-325.