GIS Ag Maps.org

______________________________________________________________

SENTINEL-2 LEVEL 1C SURFACE REFLECTANCE TUTORIAL

(For all bands except band 1 [60 meter Coastal Aerosol])

To convert to surface reflectance, follow the steps below.

1) Sentinel-2 Level 1C imagery upon download from Copernicus is in Top of Atmosphere (TOA) reflectance integer units, which makes conversion to surface reflectance quite simple. Divide any band by 10,000 with the Raster Calculator to convert the TOA four-digit integer raster to TOA reflectance units (reflectance units x 100 = reflectance percent). You may need to first convert the raster to a floating point prior to this step, or the division process may automatically convert the new raster to floating point (depends on the GIS software). A simple way to convert to a floating point raster in ArcGIS is to multiply it by 1.0 (though there is a integer to floating point tool). Surface reflectance is calculated by deducting scatter reflectance (described in later steps) from TOA reflectance.

STEPS 2 & 3 (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 STEP 4.

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

2) Step 2 involves establishing the starting scatter amount (for relative scatter) so it can be deducted from TOA reflectance (which was calculated in Step 1) 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). 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 (Frequency 50 - .008 is recommended and is described in detail below). Differences between the scatter methods below typically result in minor to modest amounts.

FREQUENCY 50 - .008 SCATTER REFLECTANCE (Recommended Method)

For Sentinel-2 (as with Landsat 8), if using ArcGIS (or other GIS software that can generate a raster attribute table), we recommend Frequency 50 - .008 scatter reflectance (which equals 0.8 % reflectance). This is a quite systematic method; the Frequency 50 - .008 method and will result in somewhat more scatter being deducted than the Base of the Histogram - .01 method because the Frequency 50 value is slightly more than the base of the histogram value (described in previous link) and the .008 deduction is slightly more than the .01 deduction. The .008 deduction amount is used, in part, because it is the same amount as the Landsat 8 Frequency 50 - .008 deduction. Also, 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). Also, the Frequency 50 value is similar to the Bin 5 value in QGIS, so similar surface reflectance can be calculated in either software. A comparison here shows that Sentinel-2 and Landsat 8 both converted to surface reflectance with the Frequency 50 method for imagery on the same date only 15 minutes apart, have similar values to each other and USGS Landsat 8 Algorithm Surface Reflectance.

Overall, Frequency 50 - .008 is recommended because it 1) is quite systematic, 2) produces a value very similar to the Bin 5 method with free QGIS, 3) is the same method used for Landsat 8, and 4) derives SR values that seem accurate for well-documented surfaces. 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.

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

For the tutorial band 4 image, from the attribute table below, the band 4 Frequency 50 value is 295, which is .0295 in reflectance units (divide 295 by 10,000 to convert to reflectance units; see the previous link for a detailed description of the Frequency 50 method). After deducting .008, the starting red band scatter is .0215. (If using the QGIS histogram, the Bin 5 value is also 295.) The .0215 value should be entered into the Relative Scatter Calculator to calculate scatter for the other bands as shown in the graphic in Step 3.

FREQUENCY 50 - .01 SCATTER REFLECTANCE

This is a hybrid of the first two methods described. It offers the systematic Frequency 50 rule to establish the scatter value from the histogram and the .01 deduction from Chavez (1996).

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

This method is directly from Chavez (1988 & 1996). The base of the histogram reflectance is used as scatter (1998) and .01 is deducted from that value so that the base of the histogram scatter is .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 Sentinel-2 (originally 12-bit, then converted to a 4-integer TOA raster) histogram does not have a long low-end statistical tail like  Landsat 8 does (16-bit pixel depth), which it makes it easier to establish a base of the histogram. The Sentinel-2 histogram is more similar to the Landsat 5 and 7 histograms. Click here to view base of the histogram pixel selection from Chavez (1988).

If using ArcGIS, you can use the histogram itself (right-click on the raster layer, then select Properties > Classified > Classify) to select a base of the histogram value, but using the attribute table is a simpler and probably better way. The low end of the band 4 (red) attribute table for the tutorial data is shown next. A reasonable base of the histogram value is 286. 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), though it would reasonable to select the first frequency 1 value, 285. For purposes here, the 286 value will be used as the base of the histogram. (See the previous link for guidence selecting the base of the histogram value.) This value represents .0286 in reflectance units (divide 286 by 10,000). Deduct .01 from the value for a starting scatter value of .0186 to establish 1% dark object reflectance (where the base of the histogram value has a surface reflectance of .01; Chavez, 1996). The .0186 value is used to establish relative scatter for other bands (described more below).

3) For purposes here, Frequency 50 - .008 scatter will be applied. Access the Relative Scatter Calculator (has an R² = 1.000 power correlation between band center wavelength and scatter reflectance [as there should be] and has proper power exponents), then deduct the scatter values corresponding to the band 4 value of 0.0215. 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

4) Divide any TOA raster by the corresponding scatter reflectance to calculate surface reflectance.

Tutorial Results Table (F50 -.008  is Frequency 50 reflectance minus .008 and is listed below)

If calculating with the Base of the Histogram minus .01 method (based on Chavez [1988 and 1996]) or the Frequency 50 minus .01 method, surface reflectance will be slightly higher because scatter is slightly lower.           

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.

Additional Sentinel-2 Band Information

Band Center Wavelengths are listed in the Sentinel-2 User Handbook (PDF; opens in new tab) starting on page 53, and are the same as in the graphics below except for band 9 which is shown as 940 nm below but is listed as 945 in the User Handbook.

Sentinel-2 10 meters Resolution Bands

Sentinel-2 20 meters Resolution Bands

Sentinel-2 60 meters Resolution Bands

(Band 9 Central Wavelength is Listed as 945 nm in User Guide)

(ESA, 2017)

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

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