Contents

(Search or Access Pages Above and Below)

Crop Agriculture

DOWNLOAD FREE QGIS HERE; ACCESS QGIS USER MANUAL HERE.

COURSE 2B - Image-Based Atmospheric Correction / Conversion of Sentinel-2 Imagery to Surface Reflectance with free QGIS (click here for Atmospheric Correction with ArcGIS).

Sentinel-2 surface reflectance can be calculated with a downloadable toolkit; this course shows how to convert to SR independently - please see the following page for toolkit access information: Background information about USGS & GIS Ag Maps surface reflectance.

Related Pages: ArcGIS Sentinel-2 SR Tutorial    Landsat 8 / Sentinel-2 Rare Comparison & Data          

* IMAGERY FOR THIS TUTORIAL CAN  BE DOWNLOADED HERE *

This course shows how to convert Sentinel-2 imagery to surface reflectance (radiance reflected from a surface / total radiance striking a surface) in free QGIS software. (The Sentinel-2 Surface Reflectance Tutorial shows how to convert if using ArcGIS.) BEING ABLE TO CONVERT IMAGERY TO SURFACE REFLECTANCE WILL ENABLE YOU TO INDEPENDENTLY PROPERLY ASSESS CROPS, OTHER VEGETATION, AND MUCH MORE. Sentinel-2 can be downloaded directly from the USGS GloVis Next website in digital number (DN) integer raster form that represents Top of Atmosphere reflectance.

IMPORTANT: As is the case with Landsat, we have found that visible band, 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 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 (based on Landsat TM) is needed for solar elevations less than 35⁰ (solar zenith angles greater than 55⁰). CLICK HERE FOR LOW SOLAR ELEVATION IMAGERY DOWNLOADS.

For this course, imagery downloaded from this website will be used; however, you can download imagery from GloVis Next (as described in Course 1D) if you prefer.

 

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

 

RECOMMENDED SCATTER CALCULATION METHODS

(Two Recommended Methods - 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)

17)  Chavez (1988 and 1996) established a 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). THE RED BAND IS RECOMMENDED FOR STARTING SCATTER. This Chavez scatter method (which is known as 1% dark object reflectance) is ONE OF THE RECOMMENDED METHODS TO ESTABLISH STARTING SCATTER FOR SENTINEL-2; SCATTER FOR OTHER BANDS SHOULD THEN BE COMPUTED WITH THE RELATIVE SCATTER CALCULATOR OR LOOKUP TABLE, AS DESCRIBED IN NEXT STEP. 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; however, if you feel the band 4 histogram is too incoherent and are not comfortable using it, use another band for the starting scatter (try the blue band next).

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 sen by viewing 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 (below), which converts to .0295 in reflectance units (divide by 10,000) and that is higher than some actual values in 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 (12-bit pixel depth) 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 (12-bit pixel depth) histograms. Click here to view base of TM histogram selection 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

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

        

 

LOWEST CONNECTED VALUE - .008

After processing the histogram properly (shown in Course 2b), the QGIS histogram should not have a break in it - if it does ignore values to the lower side of the break. In this case, select the lowest value in the histogram, which is also 295 (same as the Base of the Histogram in this case). The Lowest Connected Value and Base of Histogram should always be similar, but may not always be interpreted to be the same - see another example here. This value is close to the Frequency 50 value that can be selected with the ArcGIS attribute table (QGIS cannnot genereate an attribute table). Convert the 295 TOA integer value to reflectance by dividing by 10,000. The 295 red band value is converted to .0295. Deduct .008 from .0295 for a starting red band scatter of .0215. The .008 deduction value has been established here based on research with Landsat 8 data that shows the average and median difference between the Frequency 50 and Lowest Valid Value reflectance for 34 images throughout the year (10 of 12 months included; has a significant amount of a variety of surfaces, including agriculture, mountain, desert, and water) were both .008. This value was then applied as an option for Sentinel-2. (See the ArcGIS Landsat 8 SR Tutorial for more details.)

The .0215 value is then entered in the Relative Scatter Calculator to compute scatter for other bands (as shown in next step).

LOWEST CONNECTED VALUE -.01

This method is the same as the previous method, except .01 is deducted. This method should produce a similar or the same value (as it does in this case) as the Base of the Histogram - .01 method.

 

18) The Lowest Connected Value - .008 Method will be shown in the steps below (the same steps are used for the Base of the Histogram - .01 and Lowest Connected Value - .01 Methods). 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. (Band center wavelength values on the current calculator are slightly different than those shown below - current values are based on an the average of Sentinel 2A & 2B and are based on values listed in actual metadata files, as is explained on the calculator page.) If using the Base of the Histogram - .01 or Lowest Connected Value - .01 Methods, scatter values for bands 4, 2, 3, 5, 6, 7, 8, and 8a would both be .01950, .06075, .03685, .01562, .01304, .01056, .00803, and .00726, respectively, because the starting scatter is the same, .01950, for both methods

 Relative scatter from Relative Scatter Calculator for Sentinel-2 atmospheric correction

 

 Using Relative Scatter Calculator Provides a Proper Perfect Power Correlation

Correlation between relative scatter and Sentinel-2 band center wavelength


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. 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, or about how to apply Sentinel-2 Red Edge values.  

 

 

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.