Landsat 8 Surface Reflectance Tutorial for Free QGIS Software
(For Imagery from 2021 and Earlier; Click Here for 2022 Landsat 8 & 9 Imagery Surface Reflectance Steps.)
IMPORTANT: We now recommend using QGIS Version 3 (can be downloaded from this website here) because it can produce a raster attribute table for atmospheric correction, then applying the Frequency 50 - .008 Method for starting scatter in Step 17 below. To produce a raster attribute table in QGIS Version 3, open the Processing Toolbox, then open Raster Analysis, then open and run Raster Layer Unique Value Report (the Unique Value Report is and .html page, while the Unique Value Table is a .shp file that will open as a table in the QGIS Layers panel).
Other Surface Reflectance Guides: Landsat 8 w/ArcGIS Sentinel-2 w/ArcGIS Sentinel-2 w/ Free QGIS
FOR A LESS DETAILED TUTORIAL use: Simplified Landsat 8 Conversion to Surface Reflectance Steps
DOWNLOAD FREE QGIS HERE; ACCESS QGIS USER MANUAL HERE.
Use FREE QGIS Courses on this website to get started downloading and using QGIS software.
This easy tutorial (which is also Course 2A) 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 Free QGIS software (can be downloaded from this website by accessing link below). If you have ArcGIS software, you can use the tutorial for ArcGIS software (applies to any GIS software that can produce a raster attribute table). Landsat surface reflectance can be ordered for free from the USGS - this tutorial shows how to convert to SR independently in order to properly apply imagery. Please read: About Landsat & Sentinel-2 Surface Reflectance.
Related Page: Landsat 8 / Sentinel-2 Rare Imagery Comparison Download
* IMAGERY FOR THIS TUTORIAL CAN BE DOWNLOADED HERE *
QGIS Landsat 8 SR Tutorial Imagery Extent (LandsatLook download shown in free QGIS)
COURSE 2A - Converting FREE Landsat 8 Imagery to Surface Reflectance (SR) with FREE QGIS
* 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.
IT IS IMPORTANT TO NOTE THAT PAT CHAVEZ, JR. (RETIRED FROM THE USGS), DEVELOPED IMAGE-BASED ATMOSPHERIC CORRECTION OF SATELLITE IMAGERY AND THAT HE SHOULD BE GIVEN CREDIT AND RECOGNIZED FOR PUTTING FORTH A RELATIVELY SIMPLE METHOD (COMPARED TO DIFFERENT ALGORITHMS) OF CONVERTING SATELLITE IMAGERY TO SURFACE REFLECTANCE.
ALSO IMPORTANT TO NOTE is that while the original DOS Method (Chavez, 1988) includes the cosine of the solar zenith in the denominator when converting to SR, the COST Method (Chavez, 1996; link opens to tutorial on this website) should be applied to Landsat 5 & 7 and includes the square of the cosine of the solar zenith in the denominator as a way to increase retrieved SR (Chavez applied COST to Landsat 5, but it should also be used for Landsat 7 based on research here). While the COST Method is essential for Landsat 5 & 7, we have found that just the cosine of the solar zenith angle (not the square of the cosine) should be used for Landsat 8; reasons for this could have to due the fact that Landsat 8 has a much different sensor and much greater radiometric resolution than Landsat 5 (Chavez did not research Landsat 8 imagery). In addition to the previous link, you can download actual imagery comparing USGS Landsat 8 Algorithm SR, DOS SR, and TOA reflectance. Also, TOA is not SR but can be used as SR for bands 6 & 7 (SWIR) because there is virtually no scatter. This is explained more in the course/tutorial below:
_________________________________________________________
Landsat 8 Conversion to SR with QGIS Steps (this course/tutorial uses QGIS software; a QGIS version that can be used for this tutorial can be downloaded through the previous link):
1) If necessary, download QGIS from this website HERE or through the previous link (if you need help downloading and installing QGIS, access Free Course 1A on the Free Courses page). Add Landsat 8 band 4 (red band) image from the Tutorial Imagery Download Page. Landsat 8 scenes are about 185 x 180 kilometers (slightly more in the northerly direction; N↑), but scene sizes can vary somewhat. The image below has 41,573,559 nonzero pixels and is of the Southern California, USA, area. The San Andreas Fault is easily noticeable, running northwest to southeast in this area, essentially dividing the San Gabriel Mountains and Mojave Desert (and separates the North American and Pacific Plates); the Garlock Fault extends to the northeast from the San Andreas Fault and can also be easily seen. The Pacific Ocean is to the south (some of which is cloud-covered).
2) Click Processing (in top menu), then Toolbox to add the Processing Toolbox (if necessary). Then open GDAL/ODR > [GDAL] Conversion > Translate (convert format). THIS STEP IS VITAL.
3) The Translate (convert format) window will open. Enter 0 in the NoData window (this will eliminate the zero pixels) and select Int16 in Advanced parameters, then click Run. THIS STEP IS VITAL.
4) The new Converted layer will appear in Layers Panel. Check the box next to it and uncheck the original layer. The image will appear without the surrounding zero pixel values.
STEPS 5 THROUGH 7 ARE OPTIONAL
5) You need to determine the scatter amount in the process of converting to surface reflectance (will be explained more later). Scatter is based on the histogram. If can be helpful (but not necessary) to know the lowest value in the scene (though it is highly unlikely to be the base of the histogram). To do this open QGIS algorithms in the Toolbox, then open Raster Tools to open the Raster layer statistics tool.
6) Open the Raster layer statistics tool and click run.
7) When the tool has run, the minimum value in the raster will appear along with other values; only concern yourself with the minimum value. We have compared the minimum value generated with this tool for many scenes to ArcGIS raster attribute table low values, and they have been accurate. The minimum value here is 5,828. You will use this value later.
STEPS 8 THROUGH 19 (SHOWN AFTER THE ATMOSPHERIC SCATTER BACKGROUND) SHOW HOW TO CALCULATE ATMOSPHERIC SCATTER REFLECTANCE, WHICH WILL BE DEDUCTED FROM TOA REFLECTANCE TO DERIVE SURFACE REFLECTANCE IN LATER STEPS
ATMOSPHERIC SCATTER REFLECTANCE BACKGROUND
The path of solar radiation to Earth then to a satellite sensor involves many interactions. The following graphic (Jenson, 2007) shows various paths of incoming solar radiation (see descriptions below graphic).
Above: Ls is the total radiance at the sensor, which includes erroneous additive path radiance, Lp, of Path 2 and Path 4. Path 2 is diffuse irradiance that can include Rayleigh scattering and Mie scattering. Path 4 is reflectance from neighboring areas which is adjacency radiation. Lt is the total radiance from the actual surface/target of interest. Path 1 is irradiance that was reduce very little before striking the surface. Path 3 is irradiance that has been scattered downward onto the study area. Path 5 is radiation that has been reflected from a neighboring area, then scattered down onto the study area.
Overall, the result of the incoming solar radiation interactions is to erroneously add to Top of Atmosphere (TOA) reflectance, predominantly due to Path 2 radiance (as shown in the simplified diagram below; TOA reflectance is the ratio of total radiance at the sensor [Ls above] to total incoming irradiance at the top of the atmosphere where a satellite sensor is located). Path 2 radiance is the predominant path radiance, and it erroneously increases visible and near infrared (NIR) reflectance (though it increases NIR much less than visible reflectance). IMAGE-BASED ATMOSPHERIC CORRECTION DEDUCTS PATH 2 RADIANCE (WHICH IS THE PREDOMINANT PATH RADIANCE AND IS TYPE OF ATMOSPHERIC SCATTER) TO DERIVE SURFACE REFLECTANCE.
Path 2 can include Rayleigh scatter which is the scattering of light by molecules (much smaller than wavelengths of light) and results in smaller wavelengths scattering much more than larger wavelengths, which is why the sky is blue on a clear day (blue light has smaller wavelengths than green or red). Path 2 can also include Mie scatter (which can includes haze) and is the scattering of light by particles with a similar size to the light wavelengths; Mie results in more larger wavelength scatter than Rayleigh. Clearer days have less overall scatter but a higher ratio of blue light to red light scatter because clearer days have more Rayleigh scatter; while, conversely, hazier days have more overall scatter but a smaller ratio of blue light to red light scatter because hazier days have more Mie scatter. Near infrared radiation (not shown below) scatters a very small amount, though is should still be deducted (it is the longest wavelengths that should be deducted; see Relative Scatter Calculator page for specific values).
Fortunately, even with all the atmospheric radiation paths and interactions, a relatively simple image-based method to calculate surface reflectance quite accurately was developed, the Dark Object Subtraction (DOS) COST method, (Chavez, 1996; opens to tutorial on this website). The DOS COST method was developed with Landsat 5, and can also be applied to Landsat 7. GIS AG MAPS DEVELOPED CUSTOM DOS METHOD FOR LANDSAT 8 AND SENTINEL-2 (Chavez did not research Landsat 8 or Sentinel-2 imagery) that were shown to be accurate when applied to NDVI and compared well to USGS Landsat 8 Surface Reflectance and to each other.
The theory behind DOS atmospheric correction is that in a satellite image scene with tens of millions of pixel (which satellite imagery commonly has; Landsat has about 40 million), there should be some pixels that have zero reflectance (in visible and NIR bands); the reason there is not, is due to atmospheric scattering erroneously increasing reflectance values (mainly caused Path 2 radiance above). (DOS does not consider secondary scattering into shadowed areas [Chavez, 1996]). Chavez (1988) calculated the scatter amount by establishing a base of the Landsat TM histogram (which is considered the Dark Object and is not necessarily the lowest value); the reflectance corresponding to the base was established as the scatter reflectance amount that should be subtracted (Dark Object Subtraction [DOS]). Chavez (1996) then deducted .01 from the established scatter reflectance value (so that the established dark object had a surface reflectance of .01) 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).
Chavez (1988) suggested that continual relative scatter could be applied. GIS Ag Maps followed that suggestion and developed custom continual relative scatter calculators and lookup tables for Landsat 8 and Sentinel-2 based on principles in Chavez (1988). Per Chavez (1988), scatter values should have a power relationship (shown below) between band center wavelength and amount of scatter (as the calculators and lookup tables on this website do). The point of applying relative scatter is to establish scatter for one band (we recommend the red band), then using 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).
8) Right-click the Converted layer, then click Properties.
9) For now, click Histogram on the left to view the image histogram. You are accessing the histogram to establish the scatter amount. The histogram may appear in line or bar format.
10) Enter the minimum value from Step 7 in the Min Box then click the Max box below to move the minimum value dotted line on the histogram.
11) You can use the line format to find the base of the histogram, but for this example the bar format will be used and is recommended. If the histogram is line format, click Prefs/Action and uncheck Draw as lines.
12) A histogram in bar format will appear. The QGIS histogram does not represent/show every value as the ArcGIS histogram does, but it does displays the shape well that is may be suitable to establish a scatter value that represents the base of the histogram (where there is a relatively abrupt increase in frequencies). The steps that follow show how the QGIS histogram can be used to establish scatter.
13) Hover the cursor over the histogram and it will turn to a magnifying glass. Drag the magnifying glass over the low end of the histogram to start zoom to the base of the histogram. If you zoom too much, or want to start over for any reason, click Prefs/Action and Recompute Histogram. The histogram is zoomed in below. You will see that the lowest bin is still quite higher than the low image value that you previously entered in the Min box (5,028); the QGIS histogram does not show many values in a statistical tail (Landsat 8 commonly has long statistical tails). The base of the histogram is the area used to establish scatter (which is described in detail below). For purposes here, the values ranging from about 6,190 to 6,220 can approximate the base of the histogram - they are very low values and are not on the statistical tail too much (this process can be somewhat subjective). For purposes here, DN 6,220 will be used as the base of the histogram (described more later). The histogram will be zoomed in once more for this example.
14) You can see values at the low end of the histogram more clearly at the scale in the next graphic, including the 6,220 value from above (the upper end of the bin is more precisely located at 6,219, but 6,220 can just as well be assumed to be the base; a more systematic method is described next). After completing the VITAL STEPS 2 AND 3, we did not view a histogram that had a pattern of low values that were too fragmented to establish a scatter value.
15) Steps 15 involves establishing the starting scatter amount (for relative scatter) so it can be deducted from TOA reflectance (calculated later; in Step 23) 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. 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).
METHODS TO ESTABLISH SCATTER REFLECTANCE FOR RELATIVE SCATTER
IMPORTANT BACKGROUND INFORMATION: Chavez (1988) developed image-based atmospheric correction and 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). The process of removing scatter based on the histogram is founded on the premise that in a satellite image scene containing millions of pixels there should be pixel reflectance values that are virtually zero; the reason there is not for certain bands, is due to atmospheric scatter, which erroneously increases reflectance to higher values (shifting the entire histogram to the right). 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 scatter band to base relative scatter on (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 THREE METHODS DESCRIBED BELOW is applied (Bin 5 - .008 method is recommended and is described in detail below). Differences between the scatter methods below typically result in minor to modest amounts.
BIN 5 SCATTER REFLECTANCE - .008 (RECOMMENDED METHOD)
GIS Ag Maps has developed a systematic and repeatable scatter method based on the QGIS histogram called the Bin 5 method, where a QGIS histogram scatter value is established as the lowest value (left side) of the lowest bin with a frequency of at least 5 that does not have a bin to the right of it (higher value bin) with a frequency less than 5 (helps keeps a value on statistical tail being selected), with a .008 deduction (Bin 5 reflectance - .008). The Bin 5 value is easier to understand by viewing the example below and the examples here. At the time of this tutorial (and as recently as July of 2018), as far as we know QGIS cannot produce a raster attribute table, as ArcGIS can, so a histogram-based method has been developed for QGIS users; if QGIS has gained that capability, you can use the Frequency 50 - .008 scatter (attribute table-based) method described in the ArcGIS Tutorial, or continue to use the Bin 5 method described here.
Bin 5 - .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 Bin 5 - .008 method which we found necessary for Sentinel-2 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 Bin 5 method is systematic and produces a value very similar to the Frequency 50 method for ArcGIS (see examples in previous link). 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 (which is very similar to the Bin 5 method) 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. Deriving surface reflectance values similar to GIS Ag Maps Lowest Valid Value attribute table method values is significant because the Lowest Valid Value method is cited in Remote Sensing of Environment as producing quite accurate results for NDVI calculation
Overall, Bin 5 - .008 is recommended for QGIS because it 1) is systematic and repeatable, 2) produces a value very similar to the Frequency 50 - .008 method for ArcGIS (see previous link), 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 (particularly vegetation). 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 (which produces very similar results to the Bin 5 - .008 method), have similar values to each other and USGS Landsat 8 Algorithm Surface Reflectance.
For the tutorial histogram below, the Bin 5 value is 6191 - the lower value bins with frequencies greater than 5 are ignored because there is a bin with a frequency less than 5 greater than them (the bin with frequency 4); this bin selection method helps keeps values off a statistical tail.
BIN5 SCATTER REFLECTANCE - .01
This is a hybrid Bin 5 - .008 and Base of the Histogram - .01 method. It offers the systematic Bin 5 method to establish the scatter value from the histogram and the .01 deduction from Chavez (1996).
BASE OF THE HISTOGRAM REFLECTANCE - .01 (CHAVEZ, 1988 & 1996)
Chavez (1988) established a scatter reflectance value at the base of Landsat TM histograms (which not necessary the lowest value; Landsat TM has minimal histogram tails, unlike Landsat 8 histograms) and then deducted .01 from that value 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" (Chavez, 1996). It is very noteworthy that Landsat 8 has much longer statistical tails than Landsat TM and it is vital to select a scatter value that is not on a tail. The QGIS histogram algorithm does not included all values but is suitable for this purposes, and can help make a histogram more coherent by eliminating values on a tail. The Landsat TM images used by Chavez (1988) had a maximum value of 255, while Landsat 8 has a maximum value of 65,535 (Landsat TM & 8 scenes have similar amount of pixels, but Landsat 8 has thousands more values), which could contribute to the longer Landsat 8 tail.
15) The value of 6,220 (which was previously established as the base of the histogram) is close to 6,191. You can see where the 6,220 and 6,191 QGIS histogram values are located in the ArcGIS attribute table to the right and below (which shows the 5,828 low value), as well as differences that exist between the QGIS histogram and actual values from the ArcGIS attribute table (and the Landsat 8 statistical tail). The ArcGIS histogram (below) shows that the 6,220 value represents the base of the histogram well (where frequencies start to abruptly increase), as would the 6,191 Bin 5 value.
(Continues below)
Calculating scatter reflectance for DN integers will be explained in the next step. For now, know that the value in reflectance units for the QGIS Bin 5 DN of 6,191 is .02922, while the value for the Lowest Valid Value based on the ArcGIS attribute table DN of 5828 is .02032. In this case, the reflectance difference between the 6220 and 5828 DNs is .0089 (quite similar to the .008 value deducted here and the .01 amount that Chavez [1996] deducted from the base of the histogram, though there is variability in this amount between different images, of course).
(Continues below)
16) After establishing the Landsat 8 scatter DN, it can be converted to scatter reflectance with the following equation: ([DN x .00002] - .1) / cosine of the solar zenith; where you: 1) multiply the DN value by 0.00002, then 2) subtract 0.1 from the product, then 3) divide the result by the cosine of the solar zenith (90 - solar elevation [in degrees]). The cosine of the solar zenith can easily be computed with the calculator on the top menu or below (or from a website listed after the calculator). An advantage of using the calculator is that you just need to enter the solar (sun) elevation listed in the MTL file (included in the compressed downloaded file); if using one of the websites listed, you need to convert the solar elevation to the solar zenith. (The cosine of the solar zenith is the same as the sine of the solar elevation; click this USGS link for their explanation of Landsat 8 TOA reflectance. The Mρ value mentioned on the USGS website is .00002; the Aρ value is -.1.)
To find the solar (sun) elevation, open the uncompressed downloaded folder, then open the text .mtl file, locate the Sun Elevation, then copy it.
______________________________________________
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 can also be calculated at the following website (solar zenith = [90 - solar elevation]): http://www.rapidtables.com/calc/math/Cos_Calculator.htm; https://web2.0calc.com/
If using the Bin 5 Scatter Method where .008 is deducted from the Bin 5 reflectance value, the Bin 5 6,191 DN value is calculated to a starting red band scatter reflectance of .02122, as follows: 6,191 x .00002 -.1 = .02382; then after inputting the sun elevation of 54.60235787 into the Cosine of the Solar Zenith Calculator above, the value .81515163 is computed and .02382 is divided by .81515163 to calculate a reflectance of .02922; (3) .008 is subtracting from .02922 to calculate the starting scatter from band 4 (red band) of .02122.
If using the Chavez (1988 and 1996) Base of the Histogram minus .01 method, the 6,220 DN base of the histogram value is calculated to a starting red band scatter reflectance of .01993, as follows: 6,220 x .00002 -.1 = .02440; then after inputting the sun elevation of 54.60235787 into the Cosine of the Solar Zenith Calculator above, the value .81515163 is computed and .02440 is divided by .81515163 to calculate a reflectance of .02993; (3) .01 is subtracting from .02993 to calculate the starting scatter from band 4 (red band) of .01993.
If using the Bin 5 Scatter Method where .01 is deducted from the Bin 5 reflectance value, the Bin 5 6,191 DN value is calculated to a starting red band scatter reflectance of .01922, as follows: 6,191 x .00002 -.1 = .02382; then after inputting the sun elevation of 54.60235787 into the Cosine of the Solar Zenith Calculator above, the value .81515163 is computed and .02382 is divided by .81515163 to calculate a reflectance of .02922; (3) .01 is subtracting from .02922 to calculate the starting scatter from band 4 (red band) of .01922.
After calculating scatter reflectance, copy the value (or write in down); it will be entered in the Relative Scatter Calculator in the next step to determine relative scatter for other bands.
17) If using the Bin 5 value from Step 16, access the Landsat 8 Relative Scatter Calculator and enter the band 4 (red band) scatter amount .02122 in the calculator. After entering the starting band 4 scatter value (may need to press enter), the relative values will appear. (The calculator has values for Landsat 8 and Sentinel-2.) The Calculator is cloud-based, if it is not loading or working, use the Relative Scatter Lookup Table below it.
Using Relative Scatter Calculator Provides a Proper Perfect Power Correlation
You can copy all band center wavelengths and scatter values, paste them in a spreadsheet like Excel, and plot a power line and see that there is a perfect correlation, as well as the power exponent. The power exponent represent how clear the atmosphere is: Very Clear is -4.0; Clear is -2.0.
20) After you have determined scatter amounts, the imagery can be converted to surface reflectance quite easily. Open the Raster Calculator in QGIS by clicking Raster then Raster Calculator.
21) The Raster Calculator will appear. Enter the equation shown on the Raster Calculator image below (same equation as in Step 18; DNs x .00002 - .1). Something noteworthy we have found when working with the Raster Calculator, is that we have not been able to find a function or way to have calculated layers automatically save to the last saved location - it always saves to the same default location. If a lot of data is being saved to the default location, be sure to find that location on your hard drive and delete data you do not need.
22) The calculated layer will be added to Layers Panel and window.
23) Open the Raster Calculator again and divide the calculated raster (Course_2A_Eq_Part1 above) by the cosines of the solar zenith (.81515163) from Step 18. This will calculate Top of Atmosphere (TOA) reflectance. (The cosine of the solar zenith is the same as the sine of the solar elevation; click this USGS link for their explanation of Landsat 8 TOA reflectance.)
24) The TOA layer will appear in the window. Once you deduct scatter from this raster, you will have calculated surface reflectance for the band 4. Repeat the same steps to calculate surface reflectance for a different band (using the corresponding band scatter from Step 17.
25) To calculate surface reflectance, use the Raster Calculator to deduct the relative scatter reflectance from Step 17 from the TOA raster for each band you are interested in. Bin 5 - .008 Method scatter reflectance from Step 17 for band 2 (blue), 3 (green), 4 (red), and 5 (NIR), is as follows: .06653, .03770, .02122, and .00762, respectively.
If you are using the Chavez (1996 and 1988) based of the histogram with a .01 deduction method, surface reflectance will be slightly higher because the scatter deduction is slightly less. Scatter reflectance for this method for band 2 (blue), 3 (green), 4 (red), and 5 (NIR), is .06483, .03607, .01993, and .00692, respectively.
If you are using the Bin 5 with a .01 deduction method, surface reflectance will be slightly higher because the scatter deduction is slightly less. Scatter reflectance for this method for band 2 (blue), 3 (green), 4 (red), and 5 (NIR), is .06387, .03516, .01992, and .00655, respectively.
__________________________________________________________________________
After calculating surface reflectance it is very simple to start applying the imagery with the Raster Calculator by calculating the proper index or analyzing surface reflectance from a single band. There are tens of different equations/indices you may want to apply. Contact us if you have questions about indices.
* An important aspect to be aware of, in regards to Landsat 8 DOS SR versus Landsat 8 USGS SR, is that the USGS Landsat 8 SR algorithm is associated with apparent scatter variability (as opposed to DOS), whereby, (USGS SR - TOA) will produce varying values. Relationships in the Relative Scatter Table are based on numerous images representing mainly vegetated surfaces in agricultural and mountainous areas in the United States, but also soil surfaces, in addition to published information in Chavez (1988). *
CALCULATING SPECTRAL INDICES (Equations)
Surface reflectance is necessary to calculate indices. After imagery has been converted to surfaces reflectance, calculating indices is a simple task using the Raster Calculator.
VEGETATION
(It is important to complete or review Course 3A to understand vegetation spectral indices.)
Plant Biomass and Vigor Indices
We recommend Wide Dynamic Randge Index (WDRI) (Gitelson, 2004) for crops and other vegetation. An advantage of WDRI is that it tends to more equally weight red and NIR surface reflectance. For healthy green vegetation, NIR surface reflectance is roughly 10 times greater than red surface reflectance (multiplying NIR by a 0.1 factor helps produce more equal NIR and red values).
WDRI can be written as follows (all bands are in surface reflectance): ([NIR * 0.1] - Red) / ([NIR * 0.1] + Red)
For Sentinel-2, use Band 8a NIR.
_________________________________
Normalized Difference Vegetation Index (NDVI) (Rouse et al., 1973) is the same as WDRI, except the 0.1 value is not applied. An advantage of NDVI is that it is the most researched and well-documented.
NDVI is as follows (all bands are in surface reflectance): (NIR - Red) / (NIR + Red)
For Sentinel-2, use Band 8a NIR.
_________________________________
Plant Water Content Index
Normalized Difference Water Index (NDWI) (Gao, 1996) represents plant water content and is written as follows (all bands are in surface reflectance):
(SWIR - NIR) / (SWIR + NIR)
For Landsat 8, use Band 6 SWIR; for Landsat 7, 5, and 4, use Band 5 SWIR; for Sentinel-2, use Band 8a NIR and Band 11 SWIR.
_________________________________
Red Edge Vegetation Indices (Applies to Sentinel-2 Imagery)
The red edge is the spectral region where vegetation reflectance abruptly increases from red to NIR. Course 3A graphical shows that Sentinel-2 Band 5 (lower region of red edge) negatively correlates to Band 6 and Band 7 red edge. Clevers and Gitelson (2013) found there would be a high correlation between total crop and grass chlorophyll and nitrogen content based on the Sentinel-2 red edge Band 6:Band 5 ratio (negatively correlated; research was completed prior to Sentinel-2 acquiring imagery by deriving bands similar to Sentinel-2 based on another satellite platform).
Use an equation from the publication accessed through the previous link or get started by simply dividing (in surface reflectance) Sentinel-2 Band 6 by Band 5 with the Raster Calculator (higher values correlate to higher nitrogen and chlorophyll content); also try Band 7 divided by Band 5 (Course 3A shows there is also a high negative correlation between Band 7 and Band 5). The value of red edge applied to vegetation is quite well-documented.
WILDFIRE
The extent of a wildfire is mapped with the Differenced Normalized Burn Ratio (dNBR) which is as follows (in surface reflectance):
dNBR = NBRprefire - NBRpostfire
where, NBR = (NIR - SWIR) / (NIR + SWIR)
(For Landsat 5, 7, and 8 use Band 7 SWIR; for Sentinel-2, use Band 8a NIR and Band 12 SWIR)
SNOW
The extent of snow can be mapped with the Normalized Difference Snow Index (NDSI) (Dozier, 1989). NDSI is as follows (in surface reflectance):
(Green - SWIR) / (Green + SWIR)
(For Landsat 4, 5, and 7; and 5 use Band 5 SWIR. For Landsat 8, use Band 6 SWIR; For Sentinel-2, use Band 11 SWIR)
Though NDSI has been used to map snow, for small scale-areas we recommend using Sentinel-2 green band solely (though the blue and red band also work well) because of the fine 10-meter resolution. See the snow mapping page on this website for an example.
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.
Clevers, J.G.P.W. and A.A. Gitelson. 2013. Remote estimation of crop and grass chlorophyll and nitrogen content using red-edge bands on Sentinel-2 and -3. International Journal of Applied Earth Observation and Geoinformation 23: pp. 344–351.
Dozier, J. 1989. Spectral signature of alpine snow cover from the Landsat Thematic Mapper. Remote Sensing of Environment 28: pp. 9-22.
Gao, B.C. 1996. NDWI - A normalized difference water index for remote sensing of vegetation liquid water from space. Remote Sensing of Environment 58: pp. 257-266.
Gitelson, A.A. 2004. Wide dynamic range vegetation index for remote quantification of biophysical characteristics of vegetation. Journal of Plant Physiology 161; pp. 165–173.
Jensen, J. R. (2007). Remote sensing of the environment: An earth resource perspective (2nd ed.). Upper Saddle River, NJ: Pearson Education, Inc.
Rouse, J. W., R. H. Haas, J. A. Schell, and D. W. Deering (1973). Monitoring vegetation systems in the Great Plains with ERTS, Third ERTS Symposium, NASA SP-351 I, 309-317.