Contents

Search (over 200 pages above & below)

Crop Agriculture

Sentinel-2 Surface Reflectance Tutorial w/ Free QGIS Software (Course 2B)

Other Surface Reflectance Guides: Sentinel-2 w/ArcGIS  Landsat 8 w/ArcGIS  Landsat 8 w/FreeQGIS

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 2B) applies Image-Based Atmospheric Correction to convert Sentinel-2 imagery to surface reflectance (SR), which is essentially a process of establishing and deducting atmospheric scatter reflectance from Top of Atmosphere (TOA) reflectance. This particular guide is designed to be used with Free QGIS software (can be downloaded from page on 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). Sentinel-2 can only be downloaded in TOA integer format, but a free toolkit is available online that can be used to convert Sentinel-2 to SR - 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 Comparison & Data

* USE THE IMAGERY FOR ARCGIS TUTORIAL FOR THIS TUTORIAL; CLICK HERE TO ACCESS ARCGIS TUTORIAL AND SCROLL DOWN TO DOWNLOAD IMAGERY *

Sentinel-2 Download Extent - Northwest Ohio (Eastern Corn Belt, USA; 7/19/16 (6.2 x 6.2 miles)


 

COURSE 2B - Converting FREE Sentinel-2 Imagery to Surface Reflectance (SR) with FREE QGIS  

This course shows how to convert Sentinel-2 imagery to SR (radiance reflected from a surface / total radiance striking a surface) in free QGIS software. BEING ABLE TO CONVERT IMAGERY TO SR WILL ENABLE YOU TO INDEPENDENTLY PROPERLY ASSESS CROPS, OTHER VEGETATION, AND MUCH MORE. Sentinel-2 can be downloaded directly in digital number (DN) integer raster form that represents Top of Atmosphere reflectance.

* 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 Sentinel-2 SR with the methods described on this website have only been tested here for accuracy for solar elevations greater than 50⁰ and less than 30⁰; we 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⁰ - 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 SR COMPARISONS & DOWNLOADS.

 

For this course, Sentinel-2 imagery downloaded from this website (in Step 1 or above) will be used; however, if you prefer, you can download imagery from the Copernicus website from the Recommended Free Imagery Sources page.

 

1) Add Sentinel-2 band 4 (red band) image to QGIS - to do this, access the Sentinel-2 w/ ArcGIS SR Tutorial in the Landsat & Sentinel-2 SR Tutorials drop-down menu on the top menu (the same imagery is used for this QGIS tutorial as is used for the ArcGIS tutorial; scroll down when on page to view and access data downloads). An important aspect to realize about Sentinel-2 imagery is that they are downloaded in 110 x 110 kilometer tiles, that may or may not have NoData (zero) values (tile may be complete with imagery or have a sliver of imagery). A complete tile of 10-meter resolution imagery has 121,000,000 pixels, which is about 3 times more pixels than a Landsat scene (which has a much larger extent, but 30-meter pixel resolution). The tile downloaded from this website is full imagery (110 x 110 km; N↑). In general, try to only use tiles that are at least 1/3 complete. The imagery below is of northwest Ohio, USA, which is in the eastern Cornbelt.

2) When on the page, scroll down to the imagery downloads, and select the B4 image.

Image-Based Atmospheric Correction / Conversion of Sentinel-2 Imagery to Surface Reflectance with QGIS

 

3) The download window will appear. As stated in previous courses, make sure your internet browser setting asks you where to save downloads.

Image-Based Atmospheric Correction / Conversion of Sentinel-2 Imagery to Surface Reflectance with QGIS

 

4) Open the Sentinel-2 band 4 image in QGIS (as explained in Course 1B).

Image-Based Atmospheric Correction / Conversion of Sentinel-2 Imagery to Surface Reflectance with QGIS

 

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

Image-Based Atmospheric Correction / Conversion of Sentinel-2 Imagery to Surface Reflectance with QGIS

 

6) The Translate (convert format) window will open. Enter 0 in the NoData window (this will eliminate zero pixels if there are any; do this even if you do not think there are any NoData values) and select Int16 in Advanced parameters, then click Run. THIS STEP IS VITAL.

Image-Based Atmospheric Correction / Conversion of Sentinel-2 Imagery to Surface Reflectance with QGIS

 

7) 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/NoData pixel values (if there are any); and will be prepared to view the histogram.

Image-Based Atmospheric Correction / Conversion of Sentinel-2 Imagery to Surface Reflectance with QGIS

 

STEPS 8 THROUGH 10 ARE OPTIONAL

8) 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 low end (or base) of the histogram. If can be helpful (but not necessary) to know the lowest value in the scene (though it is 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.

Image-Based Atmospheric Correction / Conversion of Sentinel-2 Imagery to Surface Reflectance with QGIS

 

9) To do this open QGIS algorithms in the Toolbox, then open Raster Tools to open the Raster layer statistics tool.Open the Raster layer statistics tool and click run.

Image-Based Atmospheric Correction / Conversion of Sentinel-2 Imagery to Surface Reflectance with QGIS

 

10) 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 169. You will use this value later.

Image-Based Atmospheric Correction / Conversion of Sentinel-2 Imagery to Surface Reflectance with QGIS

 

11) Right-click the Converted layer, then click Properties.

Image-Based Atmospheric Correction / Conversion of Sentinel-2 Imagery to Surface Reflectance with QGIS

 

12) 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. It is recommended to use the bar histogram. If it is in line format, change it to a bar histogram as shown in the next step.

QGIS histogram for Sentinel-2 scatter for atmospheric correction and surface reflectance

 

13) Enter the minimum value from Step 10 in the Min box (in this case the value is 169), and click the cursor on the Max box to set the value. Then click Prefs/Actions and uncheck Draw as lines.

QGIS histogram for Sentinel-2 scatter for atmospheric correction and surface reflectance

 

14) A bar histogram will appear. A bar histogram would not have been possible without completing Steps 5 & 6. The QGIS histogram does not represent/show every value as the ArcGIS histogram does, but it does displays values well enough to establish a scatter value as described in the steps that follow.

QGIS histogram for Sentinel-2 scatter for atmospheric correction and surface reflectance

 

15) 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 at the low end 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, but not far enough. The QGIS histogram does not show many values in a statistical tail; Sentinel-2 is processed so there is usually not a significant tail (which is different than Landsat 8 which usually has very long statistical tails). Keep zooming to find the base of the histogram.

QGIS histogram for Sentinel-2 scatter for atmospheric correction and surface reflectance

 

16) The goal is to locate the base of the histogram. A more zoomed in view is shown in next step.

QGIS histogram low values for Sentinel-2 scatter for atmospheric correction and surface reflectance

 

 

17) Steps 17 involves establishing the starting scatter amount (for relative scatter) so it can be deducted from TOA reflectance (calculated later) 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

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 (preferably the red 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 TWO METHODS DESCRIBED BELOW is applied (Base of the Histogram - .008 scatter is recommended and is described in detail below). Differences between the scatter methods below typically result in minor to modest amounts, but you should be consistent with the method used.  

BASE OF QGIS HISTOGRAM - .008 SCATTER (RECOMMENDED METHOD - FOR SENTINEL-2 ONLY)

This is hybrid method of the Chavez (1996) Base of the Histogram Reflectance - .01 and GIS Ag Maps Bin 5 - .008. For Sentinel-2, the Base of the Histogram, Bin 5, and Lowest Connected Value methods (all described below) can all be the same value (which is the case for the tutorial data), but usually are not the same values (Sentinel-2 has modest statistical tails compared to Landsat 8).

Ultimately, a value similar to the ArcGIS Frequency 50 - .008 is sought. 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. The .008 deduction value is based on research here that showed the average and median difference between Frequency 50 reflectance and GIS Ag Maps Lowest Valid Value (Attribute Table Method) reflectance 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. GIS Ag Maps Lowest Valid Value Attribute Table method reflectance method is cited in Remote Sensing of Environment as producing quite accurate results for NDVI calculation (though NDVI was based on a scatter table with slightly different NIR values than the current table). It is a good idea to review the scatter method information given in the other tutorials (Landsat 8 and Sentinel-2 for ArcGIS and QGIS) to have a better overall understanding of why certain methods are used.

The ArcGIS attribute table and QGIS histogram shown next (both on the very low end) are not for the tutorial imagery (the tutorial imagery histogram is shown after). The information next shows a histogram example where the starting scatter should be less than the Bin 5 value (which is used to establish starting scatter for Landsat when using QGIS - the starting scatter is this example should be more similar to the Frequency 50 value of 300). In this case, the Bin 5 value of 312 is associated with a frequency that is modestly too high (represented in the histogram). This example shows a statistical tail in a QGIS histogram after the proper processing (histogram processing is described in the Tutorial above; remember the QGIS histogram does not show every low value as the ArcGIS histogram does). In a situation like this when there are very low bins, establish the base of the histogram at the lowest end of the first sub-5 frequency bin (300 in this case). Also keep in mind that QGIS histogram processing results in frequencies that are generalized based on the algorithm - you can see the the bins do not precisely correspond to the actual values in attribute table. THE TUTORIAL HISTOGRAM IS SHOWN AFTER AND HAS A BASE OF THE HISTOGRAM OF 295.

Sentinel-2 ArcGIS attribute table with scatter low values for atmospheric correction and surface reflectance.

QGIS Base of Histogram & Lowest Connected Value


 

 

 

 

 

 

 

 

 

 

 

 

 

 

The QGIS histogram for the tutorial band 4 image is shown next. The base of the histogram is 295. This is a straightforward histogram - more straightforward than usual.

QGIS histogram low value scatter for Sentinel-2 atmospheric correction and surface reflectance

(We will include more base of the histogram examples, to make it more understandable.)

 

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

This is obviously similar to the Base of Histogram -.008 method. Chavez (1988 and 1996) established reflectance scatter value at the Base of Landsat TM histograms (which was not necessary the lowest value) 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; known as 1% dark object reflectance). By completing the vital Step 6, we did not view a histogram that had low values that were too fragmented to establish a base of the histogram value.

It is important to emphasize that the QGIS histogram does not include all values; it generalizes in a way that the very lowest values are commonly not represented, but very low values are represented (as can be seen by comparing the QGIS histogram and ArcGIS attribute table [includes all values] below). You can see for the band 4 tutorial image here, the lowest histogram value is 295, which converts to .0295 in reflectance units (divide by 10,000) and that is higher than some some values that exist as is shown by the attribute table. The 169 value should be considered erroneous and ignored. The frequency 1 values comprise a very small statistical tail (and they start to break apart) and should mostly not used as base; it would reasonable to select the first frequency 1 value, 285, or 286 as the base of the histogram. If the 286 value is used, for example, the difference from the 295 base of the histogram is .0009 in reflectance units (these differences can vary). To calculate 1% dark object reflectance, .01 is deducted from .0295, for a STARTING RED BAND SCATTER of .0195.

It is worth noting that the Sentinel-2 (which is based on 12-bit pixel depth, then converted to 4-digit Integer rasters) histogram does not have the long low-end statistical tail that Landsat 8 (16-bit pixel depth) has, it is more similar to the Landsat TM (8-bit pixel depth) histograms. Click here to view base of TM histogram from Chavez (1988).

QGIS histogram low value scatter for Sentinel-2 atmospheric correction and surface reflectance

ArcGIS attribute table low value scatter for Sentinel-2 atmospheric correction and surface reflectance

 

 

 

 

 

 

 

 

 

 

 

 

   

 

 

18) The Base of the Histogram - .008 is shown in the steps below. Access the Sentinel-2 Relative Scatter Calculator and enter the scatter amount of .02150. After entering the starting band 4 scatter value (may need to press enter), the relative values will appear. 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. Scatter does not need to be deducted for bands larger than 8a (scatter amounts are irrelevant and negligible at best) - TOA can be used for SR for bands larger than 8a. 

Current band center wavelength values on the Relative Calculator are slightly different than original values which were solely based on Sentinel-2A. Current values (shown below) are based on an the average of Sentinel 2A & 2B values listed in actual metadata files (which are slightly different from each other and, therefore, slightly different than solely Sentinel-2A, as is detailed on the Relative Calculator page. Correlation values remain at least R² = 1.000). If using the Base of the Histogram minus .01 method (based on Chavez [1988 and 1996]), scatter values for bands 4, 2, 3, 5, 6, 7, 8, and 8a would be .01860, .05952, .03569, .01483, .01233, .00993, .00751, and .00677, respectively.


  Using Relative Scatter Calculator Provides a Proper Perfect Power Correlation

 

19) After establishing starting scatter reflectance, convert the TOA raster to reflectance units by dividing it by 10,000 with the raster calculator. Click Raster, the Raster Calculator.

Sentinel-2 imagery in QGIS

 

20) The Raster Calculator will open. Paste the formula below to convert the integer TOA raster to reflectance units. 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.

 QGIS raster calculator

 

21) The TOA layer in reflectance units will appear in the window. Once you deduct band 4 scatter from this band 4 TOA 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 19).

Sentinel-2 satellite imagery in QGIS software

 

22) To calculate surface reflectance for any band, deduct the scatter reflectance from the band TOA raster using the Raster Calculator. For band 2 (blue), 3 (green), 4 (red), 5 (red edge), 6 (red edge), 7 (red edge), 8 (NIR), and 8A (NIR), the values to deduct for this example are (as shown in Step 19) .07082, .04671, .02750, .02282, .01982 .01645, .01311, and .01205, respectively.

If calculating scatter with the the Lowest Connected Value - .01 or Base of the Histogram - .01 Methods, surface reflectance will be slightly higher because scatter is slightly lower (scatter is listed in Step 18).

______________________________________________________________-

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.

 

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

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 that there would be a high correlation between total crop and grass chlorophyll and N contentwith 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.

ESA. 2017. European Space Agency; Sentinel-2 Spatial Resolution. Cited at: https://earth.esa.int/web/sentinel/user-guides/sentinel-2-msi/resolutions/spatial

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.