Landsat 8 Surface Reflectance Guide/Tutorial
Other Surface Reflectance Guides: Landsat 8 w/Free QGIS Sentinel-2 w/ArcGIS Sentinel-2 w/ Free QGIS
For a less detailed tutorial use Simplified Landsat 8 & Sentinel-2 Conversion to Surface Reflectance Steps
This easy tutorial applies Image-Based Atmospheric Correction to convert Landsat 8 imagery to surface reflectance (SR), which is essentially a process of establishing and deducting atmospheric scatter reflectance from Top of Atmosphere reflectance. This particular guide is designed to be used with ArcGIS, but can be applied to any GIS software that can produce a raster attribute table. If you do not have the appropriate software, you can use the tutorial for Free QGIS Software. Landsat surface reflectance can be ordered for free from the USGS - this tutorial shows how to convert to SR independently. Please read: About Landsat & Sentinel-2 Surface Reflectance.
Related Pages: Landsat Tiers & RMSE (offsite page) Landsat 8 / Sentinel-2 Rare Comparison & Data
The Tutorial Data Follows:
The Above Data Download for Tutorial Includes: 1) Landsat 8 digital number (DN) imagery for bands 2, 3, 4, 5, and 6 (blue, green, red, NIR, and SWIR, respectively), 2) corresponding USGS algorithm surface reflectance imagery for comparison (in native integer raster format; divide by 10,000 for surface reflectance units), 3) composite bands image (click on the composite bands image with an Identify Tool and visible band tutorial surface reflectance can be seen), and 4) LandsatLook image, as well as thermal infrared and quality image, are included with Revised Tutorial data (keep in mind that the LandsatLook image can have different pixel extent than individual band imagery [if you zoom in, you will find that pixels boundaries do not precisely line up]).
Data Extent: July 30th, 2014; East Central Illinois, USA (orange extent; 23.7 km x 15.0 km)
GIS Ag Maps offers a simple and fast dark-object subtraction (DOS; Chavez [1988]) method customized for Landsat 8. Results published in Remote Sensing of Environment based on using GIS Ag Maps DOS SR for NDVI can be viewed here.
* IMPORTANT : You can use imagery as it is upon download (no processing necessary) to show relatively higher and lower areas of reflectance all year long. However, if converting to SR for visible bands, use imagery with a solar elevation > 45⁰ (for mid-latitudes in the northern hemisphere, this is from about the middle of March to the middle of September for Landsat and Sentinel-2 [opposite in southern hemisphere). THIS DOES NOT APPLY TO NIR AND SWIR BANDS, AS THEY REMAIN ACCURATE AT VERY LOW SOLAR ELEVATIONS. TOA reflectance for visible bands starts becoming too high as the solar elevation becomes lower than about 45⁰ (the shorter the wavelength, the more inaccurately high; so blue is the least accurate, followed by green then red). Chavez (1996) stated further research for the COST model is needed for solar elevations less than 35⁰ (solar zenith angles greater than 55⁰). CLICK HERE FOR LOW SOLAR ELEVATION IMAGERY SR COMPARISONS & DOWNLOADS.
_____________________________________
GIS Ag Maps DOS Surface Reflectance Tutorial
GIS Ag Maps DOS Surface Reflectance Method is outlined in the steps below and uses relative scatter based on band 4 (red) starting scatter:
STEPS TO CALCULATE LANDSAT 8 TO SURFACE REFLECTANCE FOLLOW.
First, Convert Digital Numbers to Top of Atmosphere (TOA) Reflectance in Steps 1 & 2:
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 to calculate Top of the Atmosphere (TOA) reflectance. An easy way to calculate this is to use the Cosine of Solar Zenith Calculator below or on the top menu (solar zenith = 90 - solar elevation; the scene center solar elevation is given in the .MTL file can be used or a more local value can be used if available). (The cosine of the solar zenith is the same as the sine of the solar elevation). Click this USGS link for the USGS explanation of Landsat 8 TOA reflectance (M_{ρ }value mentioned on the USGS website is .00002; the A_{ρ} value is -.1).
The Landsat scene center solar elevation for the tutorial imagery is 61.44223464. The cosine of the solar zenith angle for the data in the tutorial is 0.87833560. (In the COST Method for Landsat 5 & 7 [Chavez, 1996], the square of the cosine of the solar zenith is applied in the denominator, while this is essential for Landsat 5 & 7, we have found that the square of the cosine should not be applied for Landsat 8). TOA is not surface reflectance but should be used for bands 6 & 7 (SWIR) because there is virtually no scatter. The following calculator can be used to calculate the cosine of the solar zenith from the solar (sun) elevation (simply enter a solar elevation; the calculator will compute the zenith and cosine):
Cosine of Solar Zenith (Embedded Excel) Calculator (Cloud-Based; MAY TAKE A FEW MOMENTS TO LOAD); You Can Also Use the Lookup Table on this Website (opens in new tab).
PRESS ENTER AFTER INPUTTING SUN ELEVATION
The cosine of the solar zenith (90 - solar elevation) can also be calculated at the following website: http://www.rapidtables.com/calc/math/Cos_Calculator.htm; https://web2.0calc.com/
METHODS TO ESTABLISH SCATTER REFLECTANCE
IMPORTANT BACKGROUND INFORMATION: Chavez (1988) developed image-based atmospheric correction (Dark Object Subtraction [DOS]) and established scatter reflectance as the base of the Landsat TM histogram (which is not necessarily the lowest value) minus .01 because of the "fact that very few targets on Earth's surface are absolute black, so an assumed one-percent minimum reflectance is better than zero percent" (Chavez, 1996; the Chavez method is described more below). Chavez (1988) suggested that continual scatter could be applied. GIS Ag Maps developed Landsat 8 and Sentinel-2 scatter methods, as well as continual relative scatter calculators and tables based on principles in Chavez (1988) that are customized for Landsat 8 and Sentinel-2 (which have much greater radiometric resolution than Landsat TM). The point is to establish scatter for one band, then use the relative scatter for the other bands (explained more below).
We highly recommend the red band as the starting band for relative scatter (especially if you are using the red band in an index with NIR). Consistency in the method used is important, regardless of which of the FOUR METHODS DESCRIBED BELOW is applied (Frequency 50 - .008 method is recommended and is described in detail below). Differences between the scatter methods below typically result in minor to modest amounts.
Steps 3, 4, 5 & 6 Describe Methods to Establish and Deduct Scatter to Derive Surface Reflectance:
3) Step 3 involves establishing the starting scatter amount (for relative scatter) so it can be deducted from TOA reflectance (which was calculated in Step 2) in order to calculate surface reflectance. The red band will be used as the starting value; relative scatter for other bands should then be computed from the Landsat 8 Relative Scatter Calculator.
LOWEST VALID VALUE
The Lowest Valid Value from the attribute table establishes the lowest value that is not too disconnected from the histogram (explained in detail through previous link) - this method has retrieved good published results in Remote Sensing of Environment (though NDVI was based on a scatter table with slightly different NIR values than the current table). The Lowest Valid Value for the tutorial is 6,447 (this happens to be the lowest band 4 value in the entire scene; commonly the Lowest Valid Value is not the lowest value in the entire scene as is shown through the Lowest Valid Value link at the beginning of this paragraph).
FREQUENCY 50 - .008 (Recommended Method)
The Frequency 50 method with a .008 deduction is a simple attribute table-based method of establishing scatter (explained in detail through previous link; includes many examples). We recommend selecting a red band value from the attribute table as the starting value for relative scatter. The Frequency 50 method essentially establishes a base of the histogram (which is helpful for Landsat 8, in particular, due to long statistical tails) and also has the advantage of resulting in similar values as the QGIS Bin 5 Histogram Method (as shown in previous link and in QGIS Landsat 8 SR Tutorial). The Frequency 50 scatter value is visually interpolated or estimated from the attribute table to be the digital number where the frequency is 50 or the the closest value less than 50. There may be a small (but negligible) amount of ambiguity when selecting the value, but this will only result in an insignificant reflectance difference; a clearly suitable value will be obvious in nearly all cases.
Frequency 50 - .008 scatter is similar to the Chavez (1986 and 1996) Landsat TM histogram method where .01 is deducted from the base of the TM histogram reflectance to establish the low value as having a surface reflectance of .01, except slightly more scatter is deducted with the Frequency 50 - .008 method which we found necessary for Landsat 8 in order to lower visible band vegetation surface reflectance enough (more scatter corresponds to lower reflectance). (Chavez has established the accuracy of the COST method [1996] for Landsat TM.) The .008 deduction value has been validated by GIS Ag Maps, in part, based on research here with Landsat 8 data that shows the average and median difference between the Frequency 50 and Lowest Valid Value reflectance for 34 images throughout the year (10 of 12 months included; includes agriculture, mountain, desert, water, and other surfaces) were both .008. Establishing a base of the histogram by viewing the attribute table is challenging for Landsat 8 because there are commonly long gradual statistical tails - establishing a low value with the Frequency 50 method is more systematic and less challenging. The Frequency 50 value can be interpreted to be a base of the histogram or to be very close (see below). As shown below, the Frequency 50 value for the tutorial is 6,701.
Overall, the Frequency 50 - .008 method is recommended because it 1) is quite systematic, 2) produces a value very similar to the Bin 5 method with free QGIS, 3) derives SR values that are similar to USGS Landsat 8 Surface Reflectance Algorithm values, and 4) derives SR values that seem accurate for well-documented surfaces. A comparison here shows that Landsat 8 and Sentinel-2 imagery on the same date only 15 minutes apart both converted to surface reflectance with the Frequency 50 - .008 method, have similar values to each other and USGS Landsat 8 Algorithm Surface Reflectance.
FREQUENCY 50 SCATTER REFLECTANCE - .01
This method is a hybrid of the previous method and next method described (Base of the Histogram Reflectance - .01). It offers the systematic Frequency 50 rule to establish the scatter value from the attribute table and the .01 deduction from Chavez (1996).
BASE OF THE HISTOGRAM REFLECTANCE - .01 (CHAVEZ, 1988 & 1996)
This method is from Chavez (1988) where the starting scatter is established as the base of the histogram reflectance (Chavez, 1988; used for Landsat TM) minus .01; this established .01 as the scatter (dark object) reflectance because of the "fact that very few targets on Earth's surface are absolute black, so an assumed one-percent minimum reflectance is better than zero percent" (Chavez, 1996). Only use ArcGIS for this method if you feel confident that you can establish a base of the histogram by either viewing the attribute table or histogram. This can be challenging with Landsat 8 because there are typically much longer statistical tails than Landsat TM, which were used by Chavez (1988). We recommend using the QGIS histogram if you would like to apply this histogram method to establish a scatter value because QGIS generalizes values in a usefull way for this purpose, if processed correctly (which is explained in Course 2A; QGIS does not include every value as ArcGIS does). The QGIS histogram for the tutorial band 4 image is shown next; in graphics that follow, it is compared to the ArcGIS histogram and attribute table. For the QGIS histogram below, the base could reasonably be interpreted as anywhere from the beginning to the end of the lowest bin (6,670 to 6,696).
If viewing the histogram in ArcGIS (below), the base of the histogram can change depending on scale when viewing it. For the tutorial data, the base of the histogram value can be interpreted as 6,745 at the scale shown below (it can be a different value if zoomed in more, as is shown in the second histogram below). It can be challenging to establish Landsat 8 base of the histogram due to long statistical tails.
When the Landsat 8 histogram is zoomed into more in ArcGIS, the lowest visible bin changes, which could affect the interpretation of the base of the histogram as is shown in the next histogram (keep in mind, there are many values between the base of the histogram shown here and the lowest value in Landsat 8 histogram for the image represented by the histogram below, which is 6,447). Using the histogram to establish the Landsat 8 base of the histogram is ultimately somewhat subjective (much more than with Landsat TM in our experience). The Frequency 50 rule (previously described, and shown next) is much more systematic, which is beneficial.
FREQUENCY 50 VALUE SELECTION FOR TUTORIAL BAND 4 IMAGERY
Below is the ArcGIS attribute table and QGIS histogram for the band 4 tutorial data. Based on the attribute table, the Frequency 50 value is 6701 (though any value from 6701 to 6704 is understandable and fine; the reflectance difference is negligible and insignificant). If using the QGIS histogram (below right), the Bin 5 value is the lowest value (left side) of the bin with a frequency of at least 5 (access previous link for full explanation), which is a similar value to the ArcGIS Frequency 50 value (the Frequency 50 and QGIS Bin 5 scatter values are consistently similar). Keep in mind, when viewing the portion of the attribute table below, that the lowest value in the entire attribute table (entire scene) is 6,447 (this happens to be the Lowest Valid Value as previously mentioned - however, the Lowest Valid Value is commonly not the lowest value in the entire scene).
Applying relative scatter from the Relative Scatter Calculator assures a perfect correlation between band center wavelength and scatter (R² = 1.00) as there should be in theory (Chavez, 1988). For power correlations, power exponents range from about -4.0 to -2.0 (represents very clear and clear atmospheres, respectively [Chavez, 1988]). THE RELATIVE SCATTER CALCULATOR IS RECOMMENDED.
(While applying the Relative Scatter Calculator & Lookup Table ensures a perfect correlation between band center wavelength and scatter, the Revised Custom Landsat 8 Scatter Lookup Table ensures a high, but not perfect, power line scatter correlation between visible band center wavelength and scatter reflectance - the correlation is not perfect because it has been customized to help match USGS SR, which is not a pure DOS method.
Landsat 8 band center wavelengths for blue through NIR bands are [in micrometers]: blue=0.482; green=0.561; red=0.655; NIR=0.865). For red scatter of .01 and .05, power relationships between visible band center wavelength and scatter in the Revised Custom Landsat 8 Scatter Lookup Table vary between R² = 0.99 to R² = 0.92, respectively [correlations become lower when including NIR because the table was, in part, customized by decreasing NIR scatter to more closely match USGS scatter]. For red scatter of .01 and .05, power line exponents become larger as scatter increases [as they should], changing from -5.48 [exceptionally clear atmosphere] to -2.2 [a clear atmosphere is -2.0 (Chavez, 1988)], respectively. Important: USE RED BAND as the starting Lowest Valid Value for the Custom Scatter Table UNLESS red values are excessively disconnected [for example, if there are numerous breaks from 50 to 99]; in that case, use the blue band Lowest Valid Value as the starting value. DO NOT USE GREEN or NIR as the starting value. It is required to use the red band for the Calculator.)
The Lowest Valid Value for band 4 for the tutorial data is 6,447 (previously described). The Lowest Valid Value is a from an entire scene, while the tutorial data is just part of that scene - the Lowest Valid Value of 6,447 happens to be from an area of the scene outside the extent of the tutorial data. The Frequency 50 band 4 value is 6,701 (also previously described; it can be interpreted to be slightly different which will result in a negligible and insignificant reflectance difference).
4) For the tutorial, we will use the Lowest Valid Value and Frequency 50 - .008 band 4 values to base scatter on. (Apply the same steps that follow if you would like to calculate scatter based on one of the other methds.) Apply the same equation from Step 3 to the scatter DNs of 6,447 (Lowest Valid Value) and 6,701 (Frequency 50), which is ([DN x .00002] - 0.1). Deduct .008 from the Frequency 50 histogram (6,701 value) scatter to acquire .008 reflectance (so the value of the established scatter has .008 reflectance instead of zero reflectance [as previously explained]). (As previously mentioned, the deduction amount of .008 is based on research here; this amount is very similar to the .01 deduction value from Chavez [1996].)
5) Divide the scatter value from Step 4 by the cosine of the solar zenith (.878335597) to convert to reflectance; the reflectance value of the band 4 Lowest Valid Value Attribute Table starting scatter for the tutorial data is .032949. The band 4 Frequency 50 starting scatter is .030732 (.038732-.008). Both represent hazier than normal conditions.
6) Acquire Lowest Valid Value and Frequency 50 - .008 relative scatter values for the visible and NIR bands from the Revised Custom Landsat 8 Relative Scatter Lookup Table or the RECOMMENDED Landsat 8 Relative Scatter Calculator and deduct them from the TOA reflectance from Step 2 to calculate surface reflectance. (The composite included in the download still is based on reflectance values from the original Tutorial which is based on the Original Lookup Table.) The Calculator is cloud-based, if it is not loading or working, use the Relative Scatter Lookup Table below it. The red band is recommended as starting scatter. The Frequency 50 - .008 scatter values with Landsat 8 Relative Scatter Calculator values are shown in the next graphics. Others are listed in the Surface Reflectance Results Table near the end.
Relative Scatter Calculator Values for Starting Red Band Frequency 50 Value of .03073 (assures a proper perfect power relationship as shown in second graphic).
Power Line Plot for Relative Scatter Calculator Values
Scatter Values to Deduct from TOA Reflectance to Derived Surface Reflectance in Table Below
USGS SR is USGS Landsat 8 Algorithm Surface Reflectance. USGS SR values are in downloaded in a four-integer format (divide by 10,000 for reflectance units); for comparison purposes on the table below, Tutorial values have same amount of significant digits.
LVV SR1: Lowest Valid Value with the Landsat 8 Revised Custom Landast 8 Relative Scatter Table Deduction. Based on the Lowest Valid Value Attribute Table Method, the scatter values in the Revised Custom Lookup Table relative to the starting band 4 value of 0.032949 (.03295 falls directly between .03290 and .03300 in Lookup Table, in this case you need a 6th significant digit to determine scatter; this is not an issue with Relative Scatter Calculator) are: 0.07958 (Blue), 0.04539 (Green), and 0.00761 (NIR).
LVV SR2: Lowest Valid Value with the Relative Scatter Calculator Deduction - RECOMMENDED METHOD. Based on the Lowest Valid Value Attribute Table Method, the scatter values in the Relative Scatter Calculator relative to the starting band 4 value of 0.03295 are: 0.07999 (Blue), 0.05153 (Green), and 0.01490 (NIR).
F50 .8%: Frequency 50 Value Minus .008 with the Relative Scatter Calculator Deduction - RECOMMENDED METHOD. Based on the Frequency 50 Attribute Table Value - .008 Method, the scatter values in the Relative Scatter Calculator relative to the starting band 4 value of 0.03073: 0.07773 (Blue), 0.04905 (Green), and 0.01339 (NIR).
TOA is Top of Atmosphere reflectance (not surface reflectance).
Landsat 8 | USGS SR | LVV SR1 | LVV SR2 | F50 .8% | TOA | |
B2 (Blue) | 0.0195 | 0.0237 | 0.0233 | 0.0256 |
0.1033 | |
B3 (Green) | 0.0433 | 0.0398 | 0.0337 | 0.0362 |
0.0852 | |
B4 (Red) | 0.0258 | 0.0233 | 0.0233 | 0.0255 |
0.0562 | |
B5 (NIR) | 0.4981 | 0.4883 | 0.4810 | 0.4825 |
0.4959 | |
B6 (SWIR) | 0.1838 | 0.1794 | 0.1794 | 0.1794 | 0.1794 |
References
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.