Landsat 8 Imagery & Reflectance Tutorial  Free
Includes DOS, COST, and TOA reflectance
(Lowest Valid Value, Preferred Method [continuous relative scatter], and Dark Object haze removal; .MTL file text)
Data download includes blue, green, red, nearinfrared (NIR), shortwave infrared (SWIR), and LandsatLook imagery; guide and values to calculate reflectance are included below. See the information below about the easier Landsat 8 DOS surface reflectance method (the correct imagebased atmospheric correction method for Landsat 8).
File with RP includes recently reprocessed imagery (the USGS started reprocessing all imagery February 3, 2014, to improve accuracy; this had a minor impact for bands included in tutorial data as can be seen in the before and after values listed). Steps to calculate reprocessed imagery reflectance are included in tutorial.
Data upon opening in free Quantum GIS; LandsatLook imagery is visible while band layers are turned off.
The Landsat 8 imagery is for 7/11/13, path 22 / row 33 (located in Southern half of Illinois in Corn Belt, USA). Imagery includes LandsatLook natural color (shown above), and DNs for band 2 (blue), band 3 (green), band 4 (red), band 5 (near infrared), and band 6 (shortwave infrared). The images can be seen in Landsat 8 Color and Band Graphics page and is the same used for the example in the Landsat 8 Atmospheric Correction page. Data is saved as a free Quantum GIS map project. If opening in Quantum GIS the map will appear as above, with a LandsatLook natural color image visible. Band layers (turned off) are all grayscale with lighter shades representing higher values scaled ± 2 standard deviations from the mean. If you are using other GIS software, such as ArcGIS, the imagery can be opened as individual layers (.tif files). Surfaces represent a wide ranges of reflectance in each band and include crop and noncrop vegetation, water, clouds, shadows, and urban areas.
Landsat 8 DOS and COST Methods
(easier conversion to surface reflectance)
Whether using DOS or COST for Landsat 8, steps necessary in the atmospheric correction process can be reduced compared to previous Landsats because terms have been embedded in Landsat 8 DN values. Simply apply the steps in the following paragraph to convert to Landsat 8 DOS or COST Method surface reflectance (Landsat 8 DOS Method and Landsat 8 COST Method were published by GIS Ag Maps at www.gisagmaps.com on 9/5/2013 and are based on the DOS method from Chavez [1988] and COST method from Chavez [1996]).
To calculate Landsat 8 DOS or COST Method surface reflectance:
1) Apply the following equation to the entire raster for each band: ([DN x .00002]  0.1).
2) Divide the values from Step 1 by the cosine of the solar zenith if calculating DOS reflectance or the (cosine of the solar zenith)² for COST reflectance. The solar zenith = (90  solar elevation); the solar elevation given in the .MTL file can be used or a more local value can be used if available.
3) Establish the scatter DN (there are different methods for this as is explained in the tutorial and Landsat 8 Atmospheric Correction page; the Lowest Valid Value from band 4 [red] in conjunction with the Continuous Relative Scatter Lookup Table is the recommended method here and is described in the Preferred Method example next). If not applying relative scatter, the Lowest Valid Value for individual bands can be used.
4) Apply the following equation to the scatter DN from Step 4: ([DN x .00002]  0.1).
5) Divide the scatter value from Step 4 by the cosine of the solar zenith if calculating DOS reflectance or the (cosine of the solar zenith)² for COST reflectance (as explained in Step2).
6) Subtract .01 from scatter reflectance from all bands that have scatter more than .01 so the dark object (pixel value) has one percent reflectance (which is .01 in reflectance units) as "one percent minimum reflectance is better than zero percent" (Chavez, 1996).
7) Subtract the onepercent scatter reflectance values from corresponding band reflectances derived from Step 2 to derive DOS or COST surface reflectance. This is much simpler than the standard conversion to reflectance method (results may vary at the 100,000 place based on ESUN values listed in website).
Using dark objects (explained below) to base scatter on will produce similar values than if basing scatter on the Lowest Valid Value Method (explained below) for each band individually. However, the Preferred Method is to use the Lowest Vaild Value from band 4 (red) in conjunction with the Continuous Relative Scatter Table; an example of the Preferred Method follows (DOS is preferred overall, though both DOS and COST are calculated here):
* Preferred Method Steps  Landsat 8 DOS Method with Continuous Relative Scatter Lookup Table *
(Landsat 8 Preferred Method reflectance can be acquired as File 4 from the LEDAPS Reflectance & Index Downloads page [includes the values and steps used for reflectance calculation].) Click here to view a reflectance comparison between Landsat 8 Preferred Method and Landsat 7 LEDAPS.
The preferred method to convert Landsat 8 to surface reflectance is to base reflectance on the Landsat 8 Method steps above and establish scatter based on the Landsat 8 Continuous Relative Scatter Lookup Table; the table has been calibrated with Landsat 8 Lowest Valid Values (described later in tutorial). As far as whether DOS or COST should be used for Landsat 8, DOS is preferred here for Landsat 8 (COST is recommended over DOS for Landsat 7 and 5) mainly due to more accurate NIR values, though a thorough study is needed to confirm DOS is better for Landsat 8, such as the research that showed COST is better than DOS for Landsat 5 (Chavez, 1996). NIR values supporting DOS as the better method for Landsat 8 can be accessed through the following links: Landsat 8 DOS vs COST NIR Reflectance Landsat 8 DOS & COST vs. LEDAPS Reflectance
An importance difference when applying this preferred method is that though .01 is deducted (as is described in Step 2 above), it is deducted from values acquired from the lookup table (as is explained more below). Using the relative scatter lookup table assures that scatter values have a high power line correlation to Landsat 8 band center wavelength, which they should (see the Haze Removal and Lookup Table page for details of power correlations) and removes some uncertainty in establishing haze from different bands.
Steps to calculate reflectance for Landsat 8 DOS or COST Method with the Continuous Relative Scatter Lookup Table are as follows (DOS is the Preferred Method but steps for COST are also shown):
1) Apply the following equation to the entire raster for each band: ([DN x .00002]  0.1).
2) Divide the values from Step 1 by the cosine of the solar zenith if calculating DOS reflectance or the (cosine of the solar zenith)² for COST reflectance. The solar zenith = (90  solar elevation); the solar elevation given in the .MTL file can be used or a more local value can be used if available. The cosine of the solar zenith angle for the data in the tutorial is .90908487.
3) Establish the starting scatter value using the Lowest Valid Value (described in the previous link, below, and in the Haze Removal & Lookup Table page from band 4 (red). The lowest valid band 4 value for the tutorial data is 6022. (The Lowest Valid Value is a from an entire scene; the tutorial data is part of a scene. The Lowest Valid Value of 6022 is from an area of the scene outside the extent of the tutorial data.)
4) Apply the following equation to the scatter DN from Step 3: ([DN x .00002]  0.1).
5) Divide the scatter value from Step 4 by the cosine of the solar zenith (.90908487) to convert to reflectance; the reflectance value of the band 4 starting scatter for the tutorial data is 0.02248. Do not divide this scatter amount by the cosine of the solar zenith for COST at this point, you will divide the values from the lookup table in this next step.
6) Acquire relative scatter from the Landsat 8 Continuous Relative Scatter Lookup Table; the scatter values relative to the starting band 4 of .02248 are: 0.06975 (blue); 0.03971 (green), and 0.00766 (NIR); for COST scatter, divide the DOS values in the table by the cosine of the solar zenith. It is not necessary to deduct NIR scatter if applying onepercent dark object reflectance (which is recommended) unless the scatter is greater than .01, so in this case there will be no scatter deducted from NIR.
7) Subtract .01 from scatter reflectance for the blue, green, and red bands shown in Step 3 to calculate onepercent dark object DOS or COST surface reflectance.
8) Subtract the onepercent scatter reflectance values from the reflectance values from Step 2 to calculate surface reflectance; bands 5 or 6 have nothing subtracted.
(A Landsat 8 NIR scene commonly has many DN values that correspond to negative reflectance. Of the twentythree images assessed here, four have had any value at all above zero reflectance; none of those four have been greater than .01 reflectance [which is a realistic minimum reflectance].)
Following the steps and based on the values given, the onepercent dark object scatter for DOS for bands 2, 3, and 4 is 0.05975, 0.02971, and 0.01248, respectively. For COST, scatter for bands 2, 3, and 4 is 0.06673, 0.03368, and 0.01473, respectively. Applying the scatter values will result in the following average DOS and COST reflectance for the tutorial data:
Landsat 8 DOS Method Mean Reflectance for Tutorial Data with Lookup Table
Band 2 (blue)  0.0504  
Band 3 (green)  0.0703  
Band 4 (red)  0.0609  
Band 5 (NIR)  0.4136  
Band 6 (SWIR)  0.2055 
Landsat 8 Cost Method Mean Reflectance for Tutorial Data with Lookup Table
Band 2 (blue)  0.0544  
Band 3 (green)  0.0763  
Band 4 (red)  0.0660  
Band 5 (NIR)  0.4550  
Band 6 (SWIR)  0.2260 
* Reprocessed Landsat imagery *
Steps to calculate reflectance for Landsat 8 DOS or COST Method with the Continuous Relative Scatter Lookup Table are as follows (DOS is the Preferred Method but steps for COST are also shown):
1) Apply the following equation to the entire raster for each band: ([DN x .00002]  0.1).
2) Divide the values from Step 1 by the cosine of the solar zenith if calculating DOS reflectance or the (cosine of the solar zenith)² for COST reflectance. The solar zenith = (90  solar elevation); the solar elevation given in the .MTL file can be used or a more local value can be used if available. The cosine of the solar zenith angle for the data in the tutorial is .90908487.
3) Establish the starting scatter value using the Lowest Valid Value Method (described below and in the Haze Removal & Lookup Table page) from band 4 (red). The lowest valid band 4 value for the tutorial data is 6029. (The Lowest Valid Value is a from an entire scene; the tutorial data is part of a scene. The Lowest Valid Value of 6029 is from an area of the scene outside the extent of the tutorial data.)
4) Apply the following equation to the scatter DN from Step 3: ([DN x .00002]  0.1).
5) Divide the scatter value from Step 4 by the cosine of the solar zenith (.90908487) to convert to reflectance; the reflectance value of the band 4 starting scatter for the tutorial data is 0.02264. Do not divide this scatter amount by the cosine of the solar zenith for COST at this point, you will divide the values from the lookup table in this next step.
6) Acquire relative scatter from the Landsat 8 Continuous Relative Scatter Lookup Table; the scatter values relative to the starting band 4 of .02264 are: 0.06986 (blue); 0.03980 (green), and 0.00774 (NIR); for COST scatter, divide the DOS values in the table by the cosine of the solar zenith. It is not necessary to deduct NIR scatter if applying onepercent dark object reflectance (which is recommended) unless the scatter is greater than .01, so in this case there will be no scatter deducted from NIR.
7) Subtract .01 from scatter reflectance for the blue, green, and red bands shown in Step 3 to calculate onepercent dark object DOS or COST surface reflectance.
8) Subtract the onepercent scatter reflectance values from the reflectance values from Step 2 to calculate surface reflectance; bands 5 or 6 have nothing subtracted.
(A Landsat 8 NIR scene commonly has many DN values that correspond to negative reflectance. Of the twentythree images assessed here, four have had any value at all above zero reflectance; none of those four have been greater than .01 reflectance [which is a realistic minimum reflectance].)
Following the steps and based on the values given, the onepercent dark object scatter for DOS for bands 2, 3, and 4 is 0.05986, 0.02980, and 0.01264, respectively. For COST, scatter for bands 2, 3, and 4 is 0.06685, 0.03378, and 0.01490, respectively. Applying the scatter values will result in the following average DOS and COST reflectance for the tutorial data:
Landsat 8 DOS Method Mean Reflectance for Reprocessed Landsat Tutorial Data with Lookup Table
Band 2 (blue)  0.0508  
Band 3 (green)  0.0701  
Band 4 (red)  0.0610  
Band 5 (NIR)  0.4137  
Band 6 (SWIR)  0.2061 
Landsat 8 Cost Method Mean Reflectance for Reprocessed Landsat Tutorial Data with Lookup Table
Band 2 (blue)  0.0549  
Band 3 (green)  0.0762  
Band 4 (red)  0.0661  
Band 5 (NIR)  0.4551  
Band 6 (SWIR)  0.2267 
* For entire images from three scenes compared before and after reprocessing, TOA reflectance (not DOS or COST surface reflectance) average percent changes were the same as the USGS estimated changes prior to reprocessing of +.4 (band 2), 0 (band 3), +.3 (band 4), 0 (band 5), and +.4 (band 6) (for example, a .4 percent increase is [original reflectance x .004] added to the original reflectance).
Landsat 8 DOS and COST Traditional Method Surface Reflectance
ESUN values needed to be established for Landsat 8 in order to calculate surface reflectance with traditional methods. Correct ESUN values were calculated by GIS Ag Maps and are listed and described in the Landsat 8 ESUN Values section of the Landsat 8 ESUN, Radiance, and TOA Reflectance page. The COST model has been shown more accurate than the DOS or TOA (not surface reflectance) models with previous Landsats due to the inclusion of an additional cosine factor in the denominator (Chavez, 1996). However, Landsat has many differences compared to previous Landsats including a much different (and improved) NIR bandwidth. A comparison was made here between Landsat 8 DOS and COST with the objective of providing information to help determine the better atmospheric correction model for Landsat 8. To GIS Ag Maps knowledge, this is the first Landsat 8 reflectance comparison made between the different reflectance models. It is important to understand that more research needs to be done to make a general conclusion about whether DOS or COST DOS should be used for Landsat 8. The research here corresponds to a particular surface and condition; a conclusion as to whether DOS or COST is better for Landsat 8 should be based on research that encompasses more surfaces and conditions such as the Chavez (1996) research did that showed COST was a superior method. The research and information from GIS Ag Maps supports DOS as the better overall method for Landsat 8 surface reflectance (mainly due to more accurate NIR reflectance values). Landsat 8 COST versus DOS reflectance results follow: Landsat 8 COST vs. DOS NIR Reflectance Landsat 8 vs. Landsat 7 LEDAPS Field Reflectance Landsat 8 vs. Landsat 7 LEDAPS Reflectance.
Landsat 8 ESUN values
If calculating surface reflectance by traditional methods (though no longer necessary with Landsat 8), it is important to understand how the ESUN values included in the tutorial below were derived (the values are not necessary for Landsat 8 COST or DOS Methods); ESUN values are described in detail in the Landsat 8 ESUN, Radiance, and TOA Reflectance page.
_________________________________________________________
Tutorial for Converting Landat DNs to Reflectance
Steps to convert digital numbers (DNs) to COST, DOS, and TOA reflectance for data download
You can learn how to convert Landsat 8 digital numbers to imagebased atmospherically corrected surface reflectance by following the steps in the tutorial and converting the Landsat 8 digital numbers (DNs) in the download. Landsat 7, 5, and 4 surface reflectance data can now be ordered for free (based on the LEDAPS method; pdf opens in new tab), so only Landsat 8 needs to be converted to surface reflectance (Landsat LEDAPS surface reflectance data can be accessed through the EarthExplorer link in the Data Sources folder to the left and is the recommended reflectance data for Landsat). The tutorial involves real atmospheric correction with real Landsat 8 imagery. Landsat DN data in the download can be converted to reflectance based on different models as explained in the tutorial; the necessary values are included. You can use the raster calculator in free Quantum GIS or in other GIS software. It is a good idea to verify your calculations one way or another during the process. The steps to derive reflectance are based on the example in the Landsat 8 Atmospheric Correction page for the COST, DOS, and TOA models. Reflectance has been calculated by GIS Ag Maps and listed at the end of the tutorial; your values should be the same or very similar.
* Landsat 8 COST and DOS Methods (described above) will derive the same values (with many fewer steps) as the COST and DOS methods in the tutorial below.
The COST model (Chavez, 1996) for atmospheric correction can be written as:
(Eoλ is also referred to as ESUN; TAUv = 1.0 for Landsat; TAUz = cosθs for COST method)
The cosine of the solar zenith angle correction (COST) model (Chavez, 1996) is described in detail in the Atmospheric Correction Guide page. The "CO" is from cosine; the "S" is from solar; and the "T" is from "TAUz".
The darkobject subtraction (DOS) model (Chavez, 1988) is the above COST model with TAUz omitted.
The top of atmosphere (TOA) model (also known as apparent reflectance or inband planetary albedo) is the above COST model with the deduction for scatter and TAUz omitted. It is important to realize that TOA reflectance does not represent surface reflectance. However, for Landsat 8 the USGS (2013) has developed a particular equation based on data in the .MTL file (the text for the .MTL for the data in this tutorial is included near bottom of page) that can also be applied to derive TOA reflectance; the equation is included at the bottom of the page.
Double asterisks ( * * ) below represent points in the process where a calculation needs to be made
Step 1
Convert DNs to radiance (Lsatrad) for each band. Data for equations to convert to radiance are in .MTL file that is included with Landsat download from USGS. Equations used for the data download are shown at the end of Step 1.
Conversion to TOA Radiance per USGS (2013) is as follows:
OLI and TIRS band data can be converted to TOA spectral radiance using the radiance rescaling factors provided in the metadata file:
L_{λ} = M_{L}Q_{cal} + A_{L}
where:
L_{λ} = TOA spectral radiance (Watts/( m2 * srad * μm))
M_{L} = Bandspecific multiplicative rescaling factor from the metadata (RADIANCE_MULT_BAND_x, where x is the band number)
A_{L} = Bandspecific additive rescaling factor from the metadata (RADIANCE_ADD_BAND_x, where x is the band number)
Q_{cal} = Quantized and calibrated standard product pixel values (DN)
For the data download, DNs are converted to radiance based on the following from the .MTL file (need to check individual file for each scene as values may vary slightly):
M_{L}
RADIANCE_MULT_BAND_2 = 1.2732E02 which equals .012732
RADIANCE_MULT_BAND_3 = 1.1658E02 which equals .011658
RADIANCE_MULT_BAND_4 = 9.8736E03 which equals .0098736
RADIANCE_MULT_BAND_5 = 5.9914E03 which equals .0059914
RADIANCE_MULT_BAND_6 = 1.5095E03 which equals .0015095
A_{L}
RADIANCE_ADD_BAND_2 =  63.65864
RADIANCE_ADD_BAND_3 =  58.28984
RADIANCE_ADD_BAND_4 =  49.36793
RADIANCE_ADD_BAND_5 =  29.95701
RADIANCE_ADD_BAND_6 =  7.54767
* * Equations to covert to radiance (Lsatrad) based on the above follow (it is necessary to check the individual file for each scene as values may slightly vary). Calculation can be made with a raster calculator in GIS. Data can also be converted to a vector format in GIS and processed efficiently in a vector or spreadsheet environment if the file size is small enough. The data here have over a half million pixels. Free Quantum GIS and other GIS software, such as ArcGIS, have raster calculators.
Landsat 8 band 2 (blue): radiance = ( .012732 * DN )  63.65864
Landsat 8 band 3 (green): radiance = ( .011658 * DN )  58.28984
Landsat 8 band 4 (red): radiance = ( .0098736 * DN )  49.36793
Landsat 8 band 5 (NIR): radiance = ( .0059914 * DN )  29.95701
Landsat 8 band 6 (SWIR): radiance = ( .0015095 * DN )  7.54767
Step 2 (three parts)
Calculate Lhaze1%rad for bands for COST and DOS models. Part 3 shows amounts used for data download
Calculate DOS and COST haze radiance (also known as path radiance or scatter) and deduct the radiance that equals onepercent reflectance (based on specific equations listed below) for each band. This is explained in detail for this specific data below and in the Landat 8 Atmospheric Correction page.
Lhaze1%rad represents the amount radiance is incorrectly increased at the sensor due to the effect of atmospheric haze; this amount is commonly referred to as path radiance or scatter and needs to be deducted from total radiance of applicable bands. As solar radiation travels through space and enters Earth's atmosphere it can strike particles and reflect back to the satellite sensor, erroneously increasing radiance values. The theory behind imagebased correcting for the effect of haze is that within a Landsat scene (which contains tens of millions of pixels) there should be a surface that reflects no radiation; the amount this surface is above zero reflectance is to due to haze. However, because of " the fact that very few targets on the Earth's surface are absolute black...an assumed onepercent minimum reflectance is better than zero percent" (Chavez, 1996). Subtracting .01 reflectance units from the haze reflectance will result in the established haze (dark object) as having onepercent reflectance. Establishing the haze amount will be detailed later in this section. There are different methods to establish a haze amount; it is not a rule to use the very lowest pixel value.
Part 1
Calculate Haze Radiance (Lhaze) (also known as Path Radiance or Scatter)
There are different ways to calculate the effect of haze. Methods that applied to Landsat 4, 5, and 7 do not necessary apply to Landsat 8. Methods in this tutorial include Continuous Relative Scatter (values are listed in the Landsat 8 Haze Removal & Lookup Table page), Lowest Valid Value (as called the DN 100 Method as is described below), and Dark Object. Reflectance values based on Lowest Valid Value and Dark Object haze removal are included after last step of tutorial.
Continuous Relative Scatter vs. Lowest Valid Value
Continuous Relative scatter can be applied based on a starting onepercent reflectance value (using the Landsat 8 Cost or DOS Method to calculate surface reflectance); the corresponding bands can then have relative scatter applied (based on scientific principles and past research). Relative scatter values are listed and the method is described in the Landsat 8 Haze Removal & Lookup Table page. In regards to low values for dark object subtraction purposes, it is not a rule that the lowest pixel value in the histogram should be selected as the haze amount. A way to establish a DN scatter amount here, called the Lowest Valid Value Method, uses the lowest DN that is more than any break of values ≥ 100 in the low end of the histogram. Based on a typical solar zenith cosine value, a break of 100 corresponds to approximately .0025 reflectance units (a quarter percent); it is more likely that a DN value associated with a break such as this is affected by something nonphysical and is erroneous than a DN not associated with such a separation. A break of over 900 is shown in a histogram example in the Haze Removal and Lookup Table page. If using this method for the tables below, values would be 8189, 6942, and 6022 for the blue, green, and red bands, respectively.
Landsat 8 histogram values for, from left to right: band 2 (blue); band 3 (green); band 4 (red)
Dark Object Selection Method
Haze can be based on a true visible band dark object area or pixel. See more dark object examples here. Dark objectbased haze will result in similar values as the Lowest Valid Value method described above. For this tutorial, a dark object corresponds to an area of vegetation in cloud cover where there are very relatively low values in the red, green, and blue bands (centerright); the dots are centered on pixels where the colors represent the corresponding red, green and blue bands. The dark object method typically produces values that are very similar to the lowest valid value method described above.
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). 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 reflectance amounts 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.
(Bing imagery)
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.
The following are the haze radiance values before the onepercent reduction as detailed in the Landsat 8 Atmospheric Correction page:
Band Band Center Haze DN Haze Radiance
2 (blue) .480 8289 41.876908
3 (green) .560 6993 23.234554
4 (red) .655 6140 11.255974
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. If you prefer, to establish haze amounts simply fit low radiance values from the blue, green, and red bands (as well as the NIR band, if applicable, as is described in the main Atmospheric Correction Guide page) 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 12bit dynamic range but are delivered as 16bit images; previous Landsats were based on 8bit 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.
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. There is an example where NIR haze should be deducted in the Atmospheric Correction Guide.
Band 5 low end of histogram (5000 equals zero radiance)
As previously mentioned, bandwidths longer than NIR (such SWIR) do not need haze deducted.
As described in the Landsat 8 Atmospheric Correction page, and as is the case here, NIR may not need haze removed. (For bands longer than NIR, such as band 6 (SWIR) in this case, there is no haze deduction necessary.) For this scene, there are zero and negative radiance values for NIR in this case (there can be times where NIR or larger bands have some pixels that correspond to zero or negative radiance). If there are zero or negative radiance values, there is no haze deduction necessary. As previously mentioned, sometimes NIR will have radiance amounts over zero where a haze deduction should be made and sometimes it will not; an example of NIR haze deduction is shown in the Atmospheric Correction Guide page.
Part 2
Calculate Lhaze1%rad; this is haze radiance minus the radiance that derives onepercent reflectance
The following equations should be used to calculate the radiance that derives onepercent reflectance when input into the COST or DOS model. The amount should be deducted from the haze radiance of each band, that way the pixel value that has been established as the darkest pixel (dark object subtraction pixel) will have onepercent reflectance instead of zero when the atmospheric correction process is complete.
Onepercent reflectance equations:
COST model onepercent reflectance equation (ARSC, 2002) = 0.01 x ([Eoλ x cosθs²] / [d² x pi]);
DOS model onepercent reflectance equation (a cosθs value is omitted) = 0.01 x ([Eoλ x cosθs] / [d² x pi]).
The values to complete the onepercent reflectance equations follow:
Landsat 8 Eoλ, also known as ESUNλ, have not been made publicly available; however, they have been derived by GIS Ag Maps. The Landsat 8 ESUN, Radiance, & TOA Reflectance page describes in detail how the derived ESUN values were calculated; it is important to view the page to understand the basis of the values. The derived ESUN values for bands 2  6 (the data in the download) are as follows:
Landsat 8  
Band  ESUN  
2  2067  
3  1893  
4  1603  
5  972.6  
6  245.0  
The cosθs value is the cosine of the solar zenith; the solar zenith = 90⁰  solar elevation. The solar elevation of the scene center is included in the .MTL file (there may be an online tool to find a more local value). For the data here, cosθs = .90908487. This value needs to be squared for the COST equation but not for the DOS equation.
The d value represents the Earthsun distance and is included in the .MTL file; it should be squared in all cases for atmospheric correction. For the data here, d² = 1.0334727
pi = 3.14159265
For the COST and DOS models, the onepercent deduct amounts for bands 2, 3, and 4 are as follows:
COST  DOS  
Band  onepercent  onepercent  
2  5.261390  5.787567  
3  4.818486  5.300370  
4  4.080313  4.488374 
Part 3
For the COST and DOS models, the haze radiance from Part 1 minus the onepercent deduction amounts for bands 2 (blue), 3 (green), and 4 (NIR) from above need to be calculated and are as follows:
COST  DOS  
Band  Lhaze1%rad  Lhaze1%rad  
2  36.615518  36.089341  
3  18.416068  17.934184  
4  7.175661  6.767600 
Step 3
Calculate (Lsatrad  Lhaze1%rad)
* * To derive (Lsatrad  Lhaze1%rad), the onepercent haze amounts above (Step 2: Part 3) need to be subtracted from the radiance calculated in Step 1. Each raster of band radiance can have the Lhaze1%rad values subtracted for a COST or DOS model with a raster calculator in GIS.
Step 4
Calculate the numerator of COST or DOS model
* * Multiply the rasters from Step 3 by pi ( 3.14159265 ) and d² ( 1.0334727 [for this particular imagery]) to derived a raster that represents the numerator of the atmospheric correction equations (raster x pi x d²). For bands 5 and 6 there was no haze deduction so multiply the radiance rasters from Step 1.
For the TOA numerator ignore the haze deduction, and simply multiply radiance from Step 1 by pi and d².
Step 5
* * Calculate surface reflectance (ρλ) by dividing the numerator in Step 4 by the denominator of COST, DOS, or TOA model shown below (all values have been included below to make calculations)
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".)
The darkobject subtraction (DOS) model (Chavez, 1988) is the above COST model with TAUz omittedand 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 inband 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 midinfrared or larger because scatter is insignificant enough.
Divide the numerator raster from Step 4 by the denominator to calculate atmospherically corrected surface reflectance. This can be done with raster calculator. The difference between the COST and DOS denominators is that the COST model has (cosine of the solar zenith angle)².
For the COST model, the denominator = [ESUN value for each band shown in Step 2: Part 2] x cosθs² x d²
For the DOS and TOA models, the denominator = [ESUN value for each band shown in Step 2: Part 2] x cosθs x d²
ESUN values for bands 2  6 (shown in Step 2: Part 2) are as follows:
Landsat 8  
Band  ESUN  
2  2067  
3  1893  
4  1603  
5  972.6  
6  245.0  

The cosθs value is described and shown in Step 2: Part 2 and = .90908487 ; cosθs² [for COST] = .82643530
The d² value is described and shown in Step 2: Part 2 and = 1.0334727
Means have been calculated for the tutorial data and are as follows:
Values are in reflectance [ratio] units (multiply x 100 for percent). Visible band reflectance is low because much of the area is vegetation which absorbs over 90 percent of visible light, while NIR reflectance is much higher because of vegetation.
The values below show reflectance values that result from entering precisely the same values from the tutorial (same amount of values after decimal point). If your results are different when calculating with raster calculator, take the mean of the DN raster (in Quantum GIS, you can rightclick on raster layer, then click Properties then Metadata to find mean values) and input it into an Excel spreadsheet (or other spreadsheet) to check raster calculator results.
Mean (based on Lowest Valid Value path radiance)
COST  DOS  TOA  
Band 2  0.054  0.050  0.110  
Band 3  0.073  0.067  0.100  
Band 4  0.066  0.061  0.073  
Band 5  0.455  0.414 *  0.414 *  
Band 6  0.226  0.205 *  0.205 * 
Mean (based on Dark Object path radiance)
COST  DOS  TOA  
Band 2  0.052  0.048  0.110  
Band 3  0.072  0.066  0.100  
Band 4  0.063  0.058  0.073  
Band 5  0.455  0.414 *  0.414 *  
Band 6  0.226  0.205 *  0.205 * 
* Bands 5 (NIR) and 6 (SWIR) did not have a haze deduction as explained in the tutorial (NIR sometimes does; haze deduction is not necessary for SWIR), so the DOS and standard TOA reflectance equations are the same. TOA reflectance was calculated with the standard TOA equation instead of the USGS Landsat 8 equation (shown below). The Landsat 8 TOA equation may derive very slightly different reflectance than the standard TOA equation (values at the 10,000th place may be slightly different).
TOA Reflectance Equation for Landsat 8 (USGS, 2013)
The USGS has developed an equation (shown below) to derive TOA reflectance for Landsat 8 that can be used instead of the standard TOA equation described at top of this section.
Conversion to TOA Reflectance
OLI band data can also be converted to TOA planetary reflectance using reflectance rescaling coefficients provided in the product metadata file (MTL file). The following equation is used to convert DN values to TOA reflectance for OLI data as follows:
ρλ^{'} = M_{ρ}Q_{cal} + A_{ρ}
where:
ρλ^{'} = TOA planetary reflectance, without correction for solar angle. Note that ρλ' does not contain a correction for the sun angle.
M_{ρ} = Bandspecific multiplicative rescaling factor from the metadata (REFLECTANCE_MULT_BAND_x, where x is the band number)
A_{ρ} = Bandspecific additive rescaling factor from the metadata (REFLECTANCE_ADD_BAND_x, where x is the band number)
Q_{cal} = Quantized and calibrated standard product pixel values (DN)
TOA reflectance with a correction for the sun angle is then:
where:
ρλ = TOA planetary reflectance
θ_{SE} = Local sun elevation angle. The scene center sun elevation angle in degrees is provided in the metadata (SUN_ELEVATION).
θ_{SZ} = Local solar zenith angle; θ_{SZ} = 90°  θ_{SE}
For more accurate reflectance calculations, per pixel solar angles could be used instead of the scene center solar angle, but per pixel solar zenith angles are not currently provided with the Landsat 8 products.
Landsat 8 .MTL file text for imagery in tutorial
GROUP = L1_METADATA_FILE
GROUP = METADATA_FILE_INFO
ORIGIN = "Image courtesy of the U.S. Geological Survey"
REQUEST_ID = "0501307110212_00017"
LANDSAT_SCENE_ID = "LC80220332013192LGN00"
FILE_DATE = 20130711T19:40:19Z
STATION_ID = "LGN"
PROCESSING_SOFTWARE_VERSION = "LPGS_2.2.2"
END_GROUP = METADATA_FILE_INFO
GROUP = PRODUCT_METADATA
DATA_TYPE = "L1T"
ELEVATION_SOURCE = "GLS2000"
OUTPUT_FORMAT = "GEOTIFF"
SPACECRAFT_ID = "LANDSAT_8"
SENSOR_ID = "OLI_TIRS"
WRS_PATH = 22
WRS_ROW = 33
NADIR_OFFNADIR = "NADIR"
TARGET_WRS_PATH = 22
TARGET_WRS_ROW = 33
DATE_ACQUIRED = 20130711
SCENE_CENTER_TIME = 16:31:42.2302176Z
CORNER_UL_LAT_PRODUCT = 39.95144
CORNER_UL_LON_PRODUCT = 89.06033
CORNER_UR_LAT_PRODUCT = 39.96784
CORNER_UR_LON_PRODUCT = 86.33490
CORNER_LL_LAT_PRODUCT = 37.81705
CORNER_LL_LON_PRODUCT = 88.99956
CORNER_LR_LAT_PRODUCT = 37.83226
CORNER_LR_LON_PRODUCT = 86.35453
CORNER_UL_PROJECTION_X_PRODUCT = 324000.000
CORNER_UL_PROJECTION_Y_PRODUCT = 4424400.000
CORNER_UR_PROJECTION_X_PRODUCT = 556800.000
CORNER_UR_PROJECTION_Y_PRODUCT = 4424400.000
CORNER_LL_PROJECTION_X_PRODUCT = 324000.000
CORNER_LL_PROJECTION_Y_PRODUCT = 4187400.000
CORNER_LR_PROJECTION_X_PRODUCT = 556800.000
CORNER_LR_PROJECTION_Y_PRODUCT = 4187400.000
PANCHROMATIC_LINES = 15801
PANCHROMATIC_SAMPLES = 15521
REFLECTIVE_LINES = 7901
REFLECTIVE_SAMPLES = 7761
THERMAL_LINES = 7901
THERMAL_SAMPLES = 7761
FILE_NAME_BAND_1 = "LC80220332013192LGN00_B1.TIF"
FILE_NAME_BAND_2 = "LC80220332013192LGN00_B2.TIF"
FILE_NAME_BAND_3 = "LC80220332013192LGN00_B3.TIF"
FILE_NAME_BAND_4 = "LC80220332013192LGN00_B4.TIF"
FILE_NAME_BAND_5 = "LC80220332013192LGN00_B5.TIF"
FILE_NAME_BAND_6 = "LC80220332013192LGN00_B6.TIF"
FILE_NAME_BAND_7 = "LC80220332013192LGN00_B7.TIF"
FILE_NAME_BAND_8 = "LC80220332013192LGN00_B8.TIF"
FILE_NAME_BAND_9 = "LC80220332013192LGN00_B9.TIF"
FILE_NAME_BAND_10 = "LC80220332013192LGN00_B10.TIF"
FILE_NAME_BAND_11 = "LC80220332013192LGN00_B11.TIF"
FILE_NAME_BAND_QUALITY = "LC80220332013192LGN00_BQA.TIF"
METADATA_FILE_NAME = "LC80220332013192LGN00_MTL.txt"
BPF_NAME_OLI = "LO8BPF20130711161958_20130711163900.01"
BPF_NAME_TIRS = "LT8BPF20130711161604_20130711163953.01"
CPF_NAME = "L8CPF20130701_20130930.01"
RLUT_FILE_NAME = "L8RLUT20130211_20431231v06.h5"
END_GROUP = PRODUCT_METADATA
GROUP = IMAGE_ATTRIBUTES
CLOUD_COVER = 3.56
IMAGE_QUALITY_OLI = 9
IMAGE_QUALITY_TIRS = 9
ROLL_ANGLE = 0.001
SUN_AZIMUTH = 126.58978669
SUN_ELEVATION = 65.37919226
EARTH_SUN_DISTANCE = 1.0165986
GROUND_CONTROL_POINTS_MODEL = 340
GEOMETRIC_RMSE_MODEL = 7.228
GEOMETRIC_RMSE_MODEL_Y = 4.876
GEOMETRIC_RMSE_MODEL_X = 5.336
GROUND_CONTROL_POINTS_VERIFY = 101
GEOMETRIC_RMSE_VERIFY = 3.739
END_GROUP = IMAGE_ATTRIBUTES
GROUP = MIN_MAX_RADIANCE
RADIANCE_MAXIMUM_BAND_1 = 755.79810
RADIANCE_MINIMUM_BAND_1 = 62.41405
RADIANCE_MAXIMUM_BAND_2 = 770.71515
RADIANCE_MINIMUM_BAND_2 = 63.64591
RADIANCE_MAXIMUM_BAND_3 = 705.71509
RADIANCE_MINIMUM_BAND_3 = 58.27818
RADIANCE_MAXIMUM_BAND_4 = 597.69751
RADIANCE_MINIMUM_BAND_4 = 49.35805
RADIANCE_MAXIMUM_BAND_5 = 362.68948
RADIANCE_MINIMUM_BAND_5 = 29.95102
RADIANCE_MAXIMUM_BAND_6 = 91.37965
RADIANCE_MINIMUM_BAND_6 = 7.54616
RADIANCE_MAXIMUM_BAND_7 = 29.72563
RADIANCE_MINIMUM_BAND_7 = 2.45475
RADIANCE_MAXIMUM_BAND_8 = 673.26843
RADIANCE_MINIMUM_BAND_8 = 55.59872
RADIANCE_MAXIMUM_BAND_9 = 149.04416
RADIANCE_MINIMUM_BAND_9 = 12.30812
RADIANCE_MAXIMUM_BAND_10 = 22.00180
RADIANCE_MINIMUM_BAND_10 = 0.10033
RADIANCE_MAXIMUM_BAND_11 = 22.00180
RADIANCE_MINIMUM_BAND_11 = 0.10033
END_GROUP = MIN_MAX_RADIANCE
GROUP = MIN_MAX_REFLECTANCE
REFLECTANCE_MAXIMUM_BAND_1 = 1.210700
REFLECTANCE_MINIMUM_BAND_1 = 0.099980
REFLECTANCE_MAXIMUM_BAND_2 = 1.210700
REFLECTANCE_MINIMUM_BAND_2 = 0.099980
REFLECTANCE_MAXIMUM_BAND_3 = 1.210700
REFLECTANCE_MINIMUM_BAND_3 = 0.099980
REFLECTANCE_MAXIMUM_BAND_4 = 1.210700
REFLECTANCE_MINIMUM_BAND_4 = 0.099980
REFLECTANCE_MAXIMUM_BAND_5 = 1.210700
REFLECTANCE_MINIMUM_BAND_5 = 0.099980
REFLECTANCE_MAXIMUM_BAND_6 = 1.210700
REFLECTANCE_MINIMUM_BAND_6 = 0.099980
REFLECTANCE_MAXIMUM_BAND_7 = 1.210700
REFLECTANCE_MINIMUM_BAND_7 = 0.099980
REFLECTANCE_MAXIMUM_BAND_8 = 1.210700
REFLECTANCE_MINIMUM_BAND_8 = 0.099980
REFLECTANCE_MAXIMUM_BAND_9 = 1.210700
REFLECTANCE_MINIMUM_BAND_9 = 0.099980
END_GROUP = MIN_MAX_REFLECTANCE
GROUP = MIN_MAX_PIXEL_VALUE
QUANTIZE_CAL_MAX_BAND_1 = 65535
QUANTIZE_CAL_MIN_BAND_1 = 1
QUANTIZE_CAL_MAX_BAND_2 = 65535
QUANTIZE_CAL_MIN_BAND_2 = 1
QUANTIZE_CAL_MAX_BAND_3 = 65535
QUANTIZE_CAL_MIN_BAND_3 = 1
QUANTIZE_CAL_MAX_BAND_4 = 65535
QUANTIZE_CAL_MIN_BAND_4 = 1
QUANTIZE_CAL_MAX_BAND_5 = 65535
QUANTIZE_CAL_MIN_BAND_5 = 1
QUANTIZE_CAL_MAX_BAND_6 = 65535
QUANTIZE_CAL_MIN_BAND_6 = 1
QUANTIZE_CAL_MAX_BAND_7 = 65535
QUANTIZE_CAL_MIN_BAND_7 = 1
QUANTIZE_CAL_MAX_BAND_8 = 65535
QUANTIZE_CAL_MIN_BAND_8 = 1
QUANTIZE_CAL_MAX_BAND_9 = 65535
QUANTIZE_CAL_MIN_BAND_9 = 1
QUANTIZE_CAL_MAX_BAND_10 = 65535
QUANTIZE_CAL_MIN_BAND_10 = 1
QUANTIZE_CAL_MAX_BAND_11 = 65535
QUANTIZE_CAL_MIN_BAND_11 = 1
END_GROUP = MIN_MAX_PIXEL_VALUE
GROUP = RADIOMETRIC_RESCALING
RADIANCE_MULT_BAND_1 = 1.2485E02
RADIANCE_MULT_BAND_2 = 1.2732E02
RADIANCE_MULT_BAND_3 = 1.1658E02
RADIANCE_MULT_BAND_4 = 9.8736E03
RADIANCE_MULT_BAND_5 = 5.9914E03
RADIANCE_MULT_BAND_6 = 1.5095E03
RADIANCE_MULT_BAND_7 = 4.9105E04
RADIANCE_MULT_BAND_8 = 1.1122E02
RADIANCE_MULT_BAND_9 = 2.4621E03
RADIANCE_MULT_BAND_10 = 3.3420E04
RADIANCE_MULT_BAND_11 = 3.3420E04
RADIANCE_ADD_BAND_1 = 62.42654
RADIANCE_ADD_BAND_2 = 63.65864
RADIANCE_ADD_BAND_3 = 58.28984
RADIANCE_ADD_BAND_4 = 49.36793
RADIANCE_ADD_BAND_5 = 29.95701
RADIANCE_ADD_BAND_6 = 7.54767
RADIANCE_ADD_BAND_7 = 2.45524
RADIANCE_ADD_BAND_8 = 55.60985
RADIANCE_ADD_BAND_9 = 12.31058
RADIANCE_ADD_BAND_10 = 0.10000
RADIANCE_ADD_BAND_11 = 0.10000
REFLECTANCE_MULT_BAND_1 = 2.0000E05
REFLECTANCE_MULT_BAND_2 = 2.0000E05
REFLECTANCE_MULT_BAND_3 = 2.0000E05
REFLECTANCE_MULT_BAND_4 = 2.0000E05
REFLECTANCE_MULT_BAND_5 = 2.0000E05
REFLECTANCE_MULT_BAND_6 = 2.0000E05
REFLECTANCE_MULT_BAND_7 = 2.0000E05
REFLECTANCE_MULT_BAND_8 = 2.0000E05
REFLECTANCE_MULT_BAND_9 = 2.0000E05
REFLECTANCE_ADD_BAND_1 = 0.100000
REFLECTANCE_ADD_BAND_2 = 0.100000
REFLECTANCE_ADD_BAND_3 = 0.100000
REFLECTANCE_ADD_BAND_4 = 0.100000
REFLECTANCE_ADD_BAND_5 = 0.100000
REFLECTANCE_ADD_BAND_6 = 0.100000
REFLECTANCE_ADD_BAND_7 = 0.100000
REFLECTANCE_ADD_BAND_8 = 0.100000
REFLECTANCE_ADD_BAND_9 = 0.100000
END_GROUP = RADIOMETRIC_RESCALING
GROUP = TIRS_THERMAL_CONSTANTS
K1_CONSTANT_BAND_10 = 774.89
K1_CONSTANT_BAND_11 = 480.89
K2_CONSTANT_BAND_10 = 1321.08
K2_CONSTANT_BAND_11 = 1201.14
END_GROUP = TIRS_THERMAL_CONSTANTS
GROUP = PROJECTION_PARAMETERS
MAP_PROJECTION = "UTM"
DATUM = "WGS84"
ELLIPSOID = "WGS84"
UTM_ZONE = 16
GRID_CELL_SIZE_PANCHROMATIC = 15.00
GRID_CELL_SIZE_REFLECTIVE = 30.00
GRID_CELL_SIZE_THERMAL = 30.00
ORIENTATION = "NORTH_UP"
RESAMPLING_OPTION = "CUBIC_CONVOLUTION"
END_GROUP = PROJECTION_PARAMETERS
END_GROUP = L1_METADATA_FILE
END
References
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/ls5atmo.html. Last accessed: July, 2011.
Chavez, P.S., Jr. 1996. Imagebased atmospheric corrections–revisited and improved. Photogrammetric Engineering and Remote Sensing 62(9): pp.10251036.
Chavez, P.S., Jr. 1988. An improved darkobject subtraction technique for atmospheric scattering correction of multispectral data. Remote Sensing of Environment 24: pp.459479.
USGS. 2013. Landsat Missions: Frequently Asked Questions About the Landsat Missions. USGS. Last modified: 5/30/123. Cited at: http://landsat.usgs.gov/band_designations_landsat_satellites.php
USGS. 2013b. Using the USGS Landsat 8 Product. Cited at: http://landsat7.usgs.gov/Landsat8_Using_Product.php