GIS Ag Maps.org

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.

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:

_____________________________________

GIS Ag Maps DOS Surface Reflectance (SR) Tutorial

First Complete Steps 1 & 2, Converting from Digital Number (DN) to Top of Atmosphere (TOA) Reflectance:

1) Apply the following United States Geological Survey (USGS) equation to the entire raster for each band: ([DN x .00002] - 0.1). 

2) Also, per the USGS, divide the values from Step 1 by the cosine of the solar zenith angle (which is shown below) 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.

TOA is not surface reflectance but should be used as SR 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/

Then Complete Steps 3, 4, 5 & 6, Establishing and Deducting Atmospheric Scatter Reflectance from TOA Reflectance to Derive Surface Reflectance:

METHODS TO ESTABLISH ATMOSPHERIC SCATTER REFLECTANCE

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 METHODS 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). 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.

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

(ORIGINAL METHOD; ONE OF TWO RECOMMENDED LANDSAT 8 METHODS)

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 (OTHER 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 is 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 a recommended method 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).

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.

Jensen, J. R. (2007). Remote sensing of the environment: An earth resource perspective (2nd ed.). Upper Saddle River, NJ: Pearson Education, Inc.