**Landsat 8 Surface Reflectance Tutorial for Free QGIS Software (Course 2A)**

**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 & Sentinel-2 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 Pages: Landsat Tiers & RMSE (offsite page) Landsat 8 / Sentinel-2 Rare Comparison & Data*

*** I M A G E R Y F O R T H I S T U T O R I A L C A N B E D O W N L O A D E D H E R E** *

**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: ***If converting to SR for visible bands, try to 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). NIR and SWIR can be converted to SR at any solar elevation.* Retrieved Landsat 8 SR with the relative scatter

*been tested here for accuracy for solar elevations greater than 50⁰ and*

*with the methods described on this website have only**e have found that visible band, particularly the blue band, SR values become too high at solar elevations less than 30⁰ (because TOA becomes disproportionatly higher for shorter wavelengths at low solar elevations). We 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⁰*

*less than 30⁰; w**havez (1996) stated further research for the COST model is needed for solar elevations less than 35⁰ (solar zenith angles greater than 55⁰).*

**-**C*NIR AND SWIR BANDS REMAIN ACCURATE AT VERY LOW SOLAR ELEVATIONS.***CLICK HERE FOR LOW SOLAR ELEVATION IMAGERY SR COMPARISONS & DOWNLOADS**.

**Landsat 8 Conversion to SR with QGIS Steps **(this course uses QGIS Version 2.18.7 which can be downloaded from this website through the above link):

**1)** Download QGIS from this website **HERE** or through the link near the top of this page (if you need help downloading and installing QGIS, access **Course 1A here** or through the Course Listings link in the Free Courses drop-down menu). Add Landsat 8 band 4 (red band) image to QGIS by accessing the imagery download link near the top of the page and downloading the band 4 image. 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.

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

**method is systematic and produces a value very similar to the**

*Bin 5***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**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

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

*Frequency 50 - .008***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** = NBR

*prefire*- NBR

*postfire*

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.

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.