The QUANTIFICATION_VALUE has a true value (DN to TOA Reflectance)?

[Igor]

Dear experts! How to obtain the values of the TOA Reflectance L1C data when more than QUANTIFICATION_VALUE? Give a reply please.

Hello!
A question about a product Sentinel-2A MSI (L1C), which can be freely accessed through Sentinels Scientific Data Hub.
There is a file from the granules in JPEG2000 format.
Using the program opj_decompress from the package SNAP I converted the data to GeoTIFF format.

For example:
Product:

https://scihub.copernicus.eu/dhus/odata/v1/Products(‘92643e16-9b9e-484c-ab5b-c8a9b4f1c772’)/$value

S2A_OPER_PRD_MSIL1C_PDMC_20160104T193934_R064_V20160103T084718_20160103T084718.SAFE
TILE_ID: S2A_OPER_MSI_L1C_TL_EPA__20160104T162841_A002778_T36UXV_N02.01

./opj_decompress -i ./S2A_OPER_MSI_L1C_TL_EPA__20160104T162841_A002778_T36UXV_B02.jp2 -o ./B2.tif
gdalinfo ./B2.tif

Driver: GTiff/GeoTIFF
Files: ./B2.tif
Size is 10980, 10980
Coordinate System is `’
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( 0.0, 0.0)
Lower Left ( 0.0,10980.0)
Upper Right (10980.0, 0.0)
Lower Right (10980.0,10980.0)
Center ( 5490.0, 5490.0)
Band 1 Block=10980x1 Type=UInt16, ColorInterp=Gray

In the xml-file specified QUANTIFICATION_VALUE = 10000.
According to the documentation, to obtain the
TOA Reflectance (float) = DN / QUANTIFICATION_VALUE.
TOA Reflectance values should range from 0 to 1!

Question: what to do if the statistics of the source data is such that MAX_DN = 28508?

Sincerely,
Your user Igor.

[Igor]

The same result in SNAP 2.0:

[unnic]

Did you check where the concerning pixels are? I sometimes get “impossible” values for pixels outside the image range.

[Igor]

This is one single tile and has no background. Pay attention to the histogram. Most of the values for 10000. It really is such data and it is not only in this file. This is just an example, one of many.

[Igor]

Here’s a mask for values higher than 10000:

[Igor]

Dear forum members and all who joined the discussion.
Please excuse me for my English.
After long correspondences by mail and various discussions and thinking, I decided to summarize.

I thank all participants of this process who was able to help me.
I will be glad to continue the discussion on this issue.
So.

Here is one of my favorite responses from one of the leading experts telespazio.com, expert
Jérôme LOUIS (Earth Observation Missions Expert, Sentinel-2 Expert Support Laboratory L2A Leader):

“… Concerning your particular request, it can happen that the TOA reflectance values are higher than 1.0. It usually happens for bright, high altitude clouds or for some particular solar/viewing geometry of acquisition, e.g. specular reflection on inclined artificial roofs, it could also happen on slopes covered by snow facing the sun and the Sentinel-2 sensor. …”

Yes this is true. This explains the fact that the values of the TOA Reflectance can be greater than 1.
In support of this, you can refer to the few publications on this issue.
For example:

Yale Guide to Landsat 8 Image Processing
http://www.yale.edu/ceo/Documentation/Landsat%208%20image%20processing.pdf

On page 2 we read:
“Reflectance exceeding unity: Landsat 8 also shows a few reflectance values greater than unity; usually in bright areas of snow or clouds. These values are physically acceptable. They usually arise from bright surfaces that are oriented toward the sun: that is, a sloping snow surface or cloud top slanted towards the sun. Such a surface receives more illumination per unit area than a horizontal surface and thus appears to have a higher reflectance.”

Landsat 8: Greater than 1 Reflectance Values
(The pictures show the objects for which Reflectance > 1)
http://surfaceheat.sites.yale.edu/sites/default/files/Greater%20than%201%20Reflectance%20Values.pdf

MODIS Surface Reflectance User’s Guide
(On page 10 we see that the boundary of the original data MOD09: from -100 to 16000, and the coefficient: 0.0001. I.e. it is possible to obtain Reflectance 1.6)
https://lpdaac.usgs.gov/sites/default/files/public/product_documentation/mod09_user_guide_v1.4.pdf

And other sources…

I unfortunately have not found how to consider such values in the calculations.
But it is critical in calculating various vegetation indices!
Here’s what I suggest.

Case option
if(DN>QUANTIFICATION_VALUE) rTOA=1.0
I rejected because when a large number of such pixels, we will get a distorted image. It is very incorrect
I tried to redistribute the values in the required limits and to preserve their distribution.
It is important that these values were close to 1!
The figure shows a histogram of the distribution of values for my test image (Band 2, S2A, MSI, L1C Tile), covered with snow.

And here is the histogram after my conversion – all values are within 0 - 1:

The proposed conversion.

1. Determine the number of pixels with values greater QUANTIFICATION_VALUE (Light Pixels).
2. Get total number of pixels = rows*cols (here taken as a generalized case without considering the background, as in the test).
3. Determine the percentage of Light Pixels in the image:
   double PrBadPixels = (countBadPixels100.0)/(double)(colsrows);
4. We calculate the value of lower bound to QUANTIFICATION_VALUE, which will place the values of pixels with large values of QUANTIFICATION_VALUE:
   double Pr = QUANTIFICATION_VALUE*PrBadPixels/100.0;
   double dQ = (QUANTIFICATION_VALUE - Pr);
5. The histogram computed and determined by the level of 1% of the rightmost value - the Threshold - value exceeding QUANTIFICATION_VALUE,
   the number of which is small and can be neglected.
6. Perform conversion of values according to the following expression (DN, QUANTIFICATION_VALUE, rTOA - float values):
   if( DN < dQ) rTOA = DN / QUANTIFICATION_VALUE;
   else
   {
   if(DN > Threshold) rTOA = 1.0;
   else rTOA = (dQ + ( ( (DN - dQ)*Pr)/(Threshold-dQ) ) ) / QUANTIFICATION_VALUE;
   }

To simplify you can probably reserve for these values already predetermined range from 0.9 to 1, but I have a distribution of values was not so good. Actually it caused a preliminary determination of the percentage ratio of Light Pixels.

Best regards,

Igor Garkusha

Ukraine, Dnipropetrovsk
National Mining University

[jmalik]

Dear @Igor

By doing this you are changing the physical meaning of the reflectance values. You cannot trust the radiometric index you compute after that.

[Igor]

Dear Julien
Yes. Of course, after this transformation can no longer be used reverse the conversion to TOA radiance. For this purpose there is the original dataset and a total conversion of the DN / QUANTIFICATION_VALUE.

[jmalik]

OK so I guess I did not understand what you are trying to achieve by transforming the reflectance values like you describe.

[Igor]

How in the calculation of vegetation indices to take into account the values of the Reflectance of more than 1, or to compare with data obtained from other scanners? What do You suggest.

[jmalik]

Just discard those out-of-range pixels. You cannot interpret them as reflectances.

[Igor]

The presented test image number of pixels is about 60%. You offer them to mask and not use in the calculation?

[Igor]

On the image, the number of such pixels is about 60%. You offer them to mask and not use in the calculation?

[Igor]

Hello everyone,

My thoughts were almost correct. And I can prove it. Slightly modified their formula and brought them to minimum.
If after analyzing the histogram to determine that the Threshold value at 1% <= 1, then save it as is and correction is not needed.
If the Threshold level is 1% greater than 1, then do the conversion based on the following rules:
if( rTOA < 0 ) rTOA := 0
else
if( 0 <= rTOA < ThresholdValueOnLevel1Percent ) rTOA := rTOA/ThresholdValueOnLevel1Percent
else
if( rTOA >= ThresholdValueOnLevel1Percent ) rTOA := 1

It works. Give an example. I was chosen as the test images of Terra MODIS data (product MOD09GQ - Surface Reflectance 250 m) - band 2
and Landsat-8 OLI band 5. Accordingly, data from Sentinel-2A MSI I chose band 8.
Dates: Sentinel-2A – 24.12.2015, 03.01.2016; Terra MODIS – 24.12.2015, 03.01.2016; Landsat-8 – 29.12.2015.
All the data chosen, cover the same region:

All data are calculated values of the TOA Reflectance. All images were converted to geographic projection EPSG:4326.
Over all images was chosen as a small test area that was without clouds. For example:

Area of interest, S2A MSI, Band8 (24.12.2015):

Area of interest, S2A MSI, Band8 (03.01.2016):

Were then constructed maps of the differences of TOA Reflectance values for the three scanners.
Pre-plot with Terra MODIS data was scaled to the size of the Landsat image fragment and then the fragment image S2A.
Using GDAL utilities, own developed programs have been conducted all calculations.
In the pdf file I’m giving the histogram of values on maps the difference in TOA Reflectance.

Graphics.pdf (366.4 KB)

From the graphs it is seen that the values of the TOA Reflectance data from Landsat-8 and MOD09GQ similar in many ways.
This similarity is also manifested in the case of the snow cover. Sentinel-2 – varies considerably!
After the correction are more close to TOA Reflectance data from Landsat and MODIS.

Yes. Of course shooting is done at slightly different times. But even close to the original Sentinel is not comparable with the data MOD09GQ received on the same day.
Whatever blastema snow surface is not, it cannot cover almost 2/3 of one tile. In addition, I know for a fact that there is elevation changes.
All this is another reason to think about QUANTIFICATION_VALUE!

Note – the period of calibrating Sentinel-2A in several modes from every two weeks to once a month!

Best regards,

Igor Garkusha

Ukraine, Dnipropetrovsk
National Mining University

[Igor]

Hello everyone,

The processing algorithm implemented in the program package check_reflectance from my utilities RSUtils:
https://github.com/IgorGarkusha/RSUtils/archive/master.zip
Training materials presented at https://www.youtube.com/channel/UC13LHwOnC2qkuEWCrPinInw
Description training materials: https://github.com/IgorGarkusha/RSUtils/tree/master/tutorial

Best regards,
Igor

[kienbt2410]

Dear Igor,

my question is not relate to your topic. But, can you show me how to use the program opj_decompress from the package SNAP I converted the data to GeoTIFF format.?
I cannot see opj_decompress in SNAP by default.

Best.

[Igor]

Hello,
kienbt2410

If opj_decompress is not in the specified location, it is possible to take, for example, here:

https://github.com/senbox-org/s2tbx/tree/master/lib-openjpeg/src/main/resources/auxdata/openjpeg

or here:

uclouvain/openjpeg

( http://www.openjpeg.org/ )

Best regards,
Igor

[kienbt2410]

Dear Igor,

Thank a lot for your answer.

I am a new in remote sensing image processing. I downloaded sentinel 2 product (1C level) and tried to apply scene classification algorithm (https://sentinel.esa.int/web/sentinel/technical-guides/sentinel-2-msi/level-2a/algorithm) to 13 bands of each tile of a product.
This algorithm needs the reflectance value of each pixel to compute. I saw the method you wrote above: “TOA Reflectance (float) = DN / QUANTIFICATION_VALUE”.
So, as I know, DN is the digital number of each pixel. How is about “QUANTIFICATION_VALUE”? How can I found it in the downloaded product?

Best regards.

[Igor]

In the data directory there are different metadata.
The required value is in the xml file format of root directory of product L1C.
For example file name:
S2A_OPER_MTD_SAFL1C_PDMC_201 . . . . . . .xml

[kienbt2410]

And the Digital Number is the value that I get in pixel info of Snap?