GIS Ag Maps.org

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

 

SENTINEL-2 QGIS CONVERSION TO SURFACE REFLECTANCE STEPS

1) Download the Sentinel-2 B4 (red) image from above by clicking on it. 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) The download window will appear. As stated in previous courses, make sure your internet browser setting asks you where to save downloads. Extract the compressed download by right-clicking on it (if needed, free extraction software can be downloaded from QGIS download page).

 

3) Open QGIS (install Free QGIS, as mentioned above, if necessary).

4) Add the Sentinel-2 band 4 image to QGIS (Layer > Add Layer > Add Raster Layer; see Course 1B for information about adding imagery to QGIS, if necessary).

 

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.

 

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.

 

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.

 

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.

 

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.

 

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.

 

STEPS 11 THROUGH 18 (SHOWN AFTER THE ATMOSPHERIC SCATTER BACKGROUND) SHOW HOW TO CALCULATE ATMOSPHERIC SCATTER REFLECTANCE, WHICH WILL BE DEDUCTED FROM TOA REFLECTANCE TO DERIVE SURFACE REFLECTANCE IN LATER STEPS

ATMOSPHERIC SCATTER REFLECTANCE BACKGROUND

The path of solar radiation to Earth then to a satellite sensor involves many interactions. The following graphic (Jenson, 2007) shows various paths of incoming solar radiation (see descriptions below graphic).

Above: Ls is the total radiance at the sensor, which includes erroneous additive path radiance, Lp, of  Path 2 and Path 4. Path 2 is diffuse irradiance that can include Rayleigh scattering and Mie scattering. Path 4 is reflectance from neighboring areas which is adjacency radiation. Lt is the total radiance from the actual surface/target of interest. Path 1 is irradiance that was reduce very little before striking the surface. Path 3 is irradiance that has been scattered downward onto the study area. Path 5 is radiation that has been reflected from a neighboring area, then scattered down onto the study area.  

Overall, the result of the incoming solar radiation interactions is to erroneously add to Top of Atmosphere (TOA) reflectance, predominantly due to Path 2 radiance (as shown in the simplified diagram below; TOA reflectance is the ratio of total radiance at the sensor [Ls above] to total incoming irradiance at the top of the atmosphere where a satellite sensor is located). Path 2 radiance is the predominant path radiance, and it erroneously increases visible and near infrared (NIR) reflectance (though it increases NIR much less than visible reflectance). IMAGE-BASED ATMOSPHERIC CORRECTION DEDUCTS PATH 2 RADIANCE (WHICH IS THE PREDOMINANT PATH RADIANCE AND IS TYPE OF ATMOSPHERIC SCATTER) TO DERIVE SURFACE REFLECTANCE.

Path 2 can include Rayleigh scatter which is the scattering of light by molecules (much smaller than wavelengths of light) and results in smaller wavelengths scattering much more than larger wavelengths, which is why the sky is blue on a clear day (blue light has smaller wavelengths than green or red). Path 2 can also include Mie scatter (which can includes haze) and is the scattering of light by particles with a similar size to the light wavelengths; Mie results in more larger wavelength scatter than Rayleigh. Clearer days have less overall scatter but a higher ratio of blue light to red light scatter because clearer days have more Rayleigh scatter; while, conversely, hazier days have more overall scatter but a smaller ratio of blue light to red light scatter because hazier days have more Mie scatter. Near infrared radiation (not shown below) scatters a very small amount, though is should still be deducted (it is the longest wavelengths that should be deducted; see Relative Scatter Calculator page for specific values).

Fortunately, even with all the atmospheric radiation paths and interactions, a relatively simple image-based method to calculate surface reflectance quite accurately was developed, the Dark Object Subtraction (DOS) COST method, (Chavez, 1996; opens to tutorial on this website). The DOS COST method was developed with Landsat 5, and can also be applied to Landsat 7. GIS AG MAPS DEVELOPED CUSTOM DOS METHOD FOR LANDSAT 8 AND SENTINEL-2 (Chavez did not research Landsat 8 or Sentinel-2 imagery) that were shown to be accurate when applied to NDVI and compared well to USGS Landsat 8 Surface Reflectance and to each other.

The theory behind DOS atmospheric correction is that in a satellite image scene with tens of millions of pixel (which satellite imagery commonly has; Landsat has about 40 million), there should be some pixels that have zero reflectance (in visible and NIR bands); the reason there is not, is due to atmospheric scattering erroneously increasing reflectance values (mainly caused Path 2 radiance above). (DOS does not consider secondary scattering into shadowed areas [Chavez, 1996]). Chavez (1988) calculated the scatter amount by establishing a base of the Landsat TM histogram  (which is considered the Dark Object and is not necessarily the lowest value); the reflectance corresponding to the base was established as the scatter reflectance amount that should be subtracted (Dark Object Subtraction [DOS]). Chavez (1996) then deducted .01 from the established scatter reflectance value (so that the established dark object had a surface reflectance of .01) because of "the fact that very few targets on the Earth's surface are absolute black, so an assumed one-percent minimum reflectance is better than zero percent" (Chavez, 1996).

Chavez (1988) suggested that continual relative scatter could be applied. GIS Ag Maps followed that suggestion and developed custom continual relative scatter calculators and lookup tables for Landsat 8 and Sentinel-2 based on principles in Chavez (1988). Per Chavez (1988), scatter values should have a power relationship (shown below) between band center wavelength and amount of scatter (as the calculators and lookup tables on this website do). The point of applying relative scatter is to establish scatter for one band (we recommend the red band), then using the relative scatter for the other bands (explained more below). We highly recommend the red band as the starting band for relative scatter (especially if you are using the red band in an index with NIR).

 

IMPORTANT: We recommend using QGIS Version 3 (can be downloaded from this website here) because it can produce a raster attribute table for atmospheric correction, then applying the  Frequency 50 - .008 Attribute Table Method for starting scatter in Step 18 below. To produce a raster attribute table in QGIS Version 3, open the Processing Toolbox, then open Raster Analysis, then open and run Raster Layer Unique Value Report (the Unique Value Report is and .html page, while the Unique Value Table is a .shp file that will open as a table in the QGIS Layers panel). IF USING AN ATTRIBUTE TABLE METHOD TO ESTABLISH A SCATTER DN, YOU CAN SKIP TO STEP 18.

 

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

 

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.

 

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.

 

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.

 

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.

Important: Sentinel-2 pixel values represent Top of Atmosphere (TOA) reflectance units x 10,000 (for example, a pixel value of 5000 = 0.5000 TOA reflectance units [5000/10000], a pixel value of 250 = .0250 reflectance, a pixel value of 3482 = 0.3482 reflectance, and so on).

 

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

 

 

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

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.  

 

* RECOMMENDED SCATTER METHOD FOR SENTINEL-2 (WHEN USING QGIS ONLY)

SENTINEL-2 QGIS BASE OF QGIS HISTOGRAM - .008 (includes base of histogram rules) *

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 correspond to the same value (which is the case for the tutorial data), but are usually 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.

Important: Sentinel-2 pixel values represent Top of Atmosphere (TOA) reflectance units x 10,000 (for example, a pixel value of 5000 = 0.5000 TOA reflectance units [5000/10000], a pixel value of 250 = .0250 reflectance, a pixel value of 3482 = 0.3482 reflectance, and so on).

 

 

 

 

 

 

 

 

 

 

 

 

 

 

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.

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

 

 

 

 

 

 

 

 

 

 

 

 

   

 

 

 

 

 

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.

 

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.

 

 

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

 

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.

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

Rouse, J. W., R. H. Haas, J. A. Schell, and D. W. Deering (1973). Monitoring vegetation systems in the Great Plains with ERTS, Third ERTS Symposium, NASA SP-351 I, 309-317.