DOWNLOAD FREE QGIS **HERE**; ACCESS QGIS USER MANUAL **HERE**.

**COURSE 2A - Image-Based Atmospheric Correction / Conversion of Landsat 8 Imagery to Surface Reflectance with QGIS (click here for Landsat 8 Conversion to Surface Reflectance with ArcGIS)**

*** IMAGERY FOR THIS TUTORIAL CAN BE DOWNLOADED HERE** *

*Landsat surface reflectance can be ordered for free from the USGS. This tutorial shows how to convert to SR independently* - please read:** ****Background information about USGS & GIS Ag Maps surface reflectance.**

* *

**Relevant Pages: ****Landsat Tiers & RMSE** (offsite page)** Landsat 8 / Sentinel-2 Rare Comparison & Data**

*This course shows how to convert Landsat 8 imagery to surface reflectance (radiance reflected from a surface / total radiation striking a surface) in free QGIS software. BEING ABLE TO CONVERT IMAGERY TO SURFACE REFLECTANCE WILL ENABLE YOU TO INDEPENDENTLY PROPERLY APPLY IMAGERY. Landsat 8 can be downloaded directly from the USGS GloVis Next website in digital number (DN) integer raster form. *

**IMPORTANT: ***This Landsat 8 DOS SR method has only been tested here for accuracy for solar elevations greater than 50⁰. We have found that visible bands, particularly the blue band, SR values become much less accurate at solar elevations less than 30⁰, but have not completed research for solar elevations from 30 to 50⁰. It is very likely that visible band retrieved SR, particularly the blue band, becomes increasingly less accurate (too high) as solar elevations decrease from about 45⁰ - NIR AND SWIR BANDS REMAIN ACCURATE AT VERY LOW SOLAR ELEVATIONS. 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 DOWNLOADS.*

**Landsat 8 Conversion to Surface Reflectance with QGIS **(this course uses QGIS Version 2.18.7, but you can follow along if using ArcGIS or another GIS software):

**1)** Add Landsat 8 band 4 (red band) image to QGIS by accessing the imagery download link near the top of the page. Scenes are about 185 x 180 kilometers miles (slightly more in the northerly direction; N↑), but scene sizes can vary somewhat, as the amount of nonzero pixels (black border area) can vary. The image below has 41,573,559 nonzero pixels. The image below is of the Southern California 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 base of 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.

**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. 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 currently be used.

**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 scatter value. 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 for a dataset like Landsat. For purposes here, the values ranging from about 6,190 to 6,220 can approximate the base of the histogram - they represent low values and are not in the statistical tail too much (though it is on the tail somewhat; this process can be somewhat subjective) - for purposes here, **DN 6,220 will be used as the base of the histogram**. 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 graphic below, 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; however, if you feel the band 4 histogram is too incoherent and are not comfortable using it, use another band for the starting scatter.

**RECOMMENDED SCATTER CALCULATION METHODS**

(Three Recommended Methods; all methods results in similar scatter amounts - be consistent with method of choice.

**IMPORTANT:** *We 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), but it is not a requirement (another visible band can be used). Consistency in the method used is important, however.*

**BASE OF THE HISTOGRAM REFLECTANCE - .01 (CHAVEZ, 1988 & 1996)**

Chavez (1988) established a scatter 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). **This Chavez scatter method (which is also known as 1% dark object reflectance) IS ONE RECOMMENDED METHOD if using the QGIS histogram**; however, 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 (scenes have similar amount of pixels, but Landsat 8 has thousands more values), which could contribute to the longer Landsat 8 tail. **Another recommended scatter method is described next.*

*BIN 5* SCATTER REFLECTANCE - .008

GIS AG MAPS has developed a systematic scatter method called the *BIN 5* METHOD, where a QGIS histogram value to base scatter on 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; easier to understand by viewing example below and **examples here**).

**In the 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. Advantages of establishing scatter with this method are that it is systematic and that the scatter reflectance derived will be similar to the ArcGIS

**(as is described in the ArcGIS tutorial).**

*Frequency 50*methodThe Bin 5 scatter DN is then converted to reflectance and .008 is deducted to calculate starting scatter. This is a slightly different deduction than the Chavez (1996) .01 amount. The .008 amount is based on research here with Landsat 8 data that shows the average and median difference between *Frequency 50 **reflectance (which is similar to the Bin 5 reflectance) and GIS Ag Maps*

*for 34 images throughout the year (imagery was for 10 of 12 months; included a significant amount of different surfaces, including agriculture, mountain, desert, and water) were both .008. The*

**Lowest Valid Value Attribute Table Method**reflectance*Bin 5*reflectance with a .008 deduction is sought to estimate the GIS Ag Maps Lowest Valid Value Attribute Table method reflectance

*because the method is*

**cited in**

*Remote Sensing of Environment*as producing quite accurate results for ND**VI calculation**.

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

**15)** The value of 6,220 is close to 6,191. You can see where the 6,220 and 6,191 QGIS histogram values from above 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 .01 amount that Chavez [1996] deducted from the base of the histogram and the .008 value deducted here; though there is variability in this amount 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 be computed with the calculator 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):**

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. The red band is recommended starting scatter band here; however, if the red band low values seem very incoherent and you feel another band has a more valid starting value, use that preferred band and the Lookup Table, as any band can be the starting band when using the Lookup Table. *(The calculator assures a proper perfect power relationship as shown in second graphic; Landsat 8 band 2 and 3 center wavelenghths have change very slightly on the current calculator based on recent information.)*

**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, deduct the scatter reflectance from Step 17 from the TOA raster for each band you are interested in using the Raster Calculator. *Bin 5* minus .008 scatter reflectance from Step 17 for band 2 (blue), 3 (green), 4 (red), and 5 (NIR), is .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.* There are tens of different equations/indices you may want to apply. *Contact us if you have questions about indices.*

**Please Read the Following if you are Applying Imagery to Crops:**

** ***There will be a course posted in the future that shows how to use imagery for crop condition and yield assessment purposes (will also include extracting imagery to the extent of a field). *Until then, if you are planning on calculating and vegetation spectral index, such as the Wide Dynamic Range Index (WDRI) (Gitelson, 2004), where α=.1 (the recommended index here), or others such as NDVI or MSAVI2, keep in mind it should only be used when red and NIR surface reflectance (SR) negatively correlate (areas of higher NIR values correspond to areas of lower red values) – the healthier the vegetation the higher the NIR SR and the lower the red (and all other visible band) SR. If there is a positive correlation between NIR and red, too much soil is visible, as there is a positive correlation between NIR and all visible band SR, including the red band, for soil. The point in the season when soil is covered enough, is when the correlation be NIR and visible bands flip from positive to negative. There have been indices and methods developed to help account for the soil background effect, but it is better to avoid, if possible, using imagery where too much soil is visible. Also, for certain crops, such as soybeans, there may be a point in the season where red SR saturates low (meaning that is becomes insensitive to crop differences throughout a field, as all values have virtually the same low reflectance). If this occurs, the NIR band, solely, is likely to continue to show crop condition patterns. Also, keep in mind that there are vegetation indices, that include other visible band, mainly green, and NIR SR, but the same principles described here about the red band apply to any other visible band. *Contact GIS Ag Maps if you have questions about imagery and crop or yield patterns. *

* *

**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.

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.