0

When deriving NDVI out of a Sentinel-2 image - S2A_MSIL2A_20190502T074621_N0211_R135_T38SLF_20190502T102541 - I am getting a strange result when doing it in python vs in QGIS.

In QGIS (using raster calculator) the range of values is as expected, between -1 and +1, which I can confirm from the product metadata:

enter image description here

When doing the NDVI calculation in python usign the following code:

NIR = rasterio.open(S_files[3]).read() # band 8
Red = rasterio.open(S_files[2]).read() # band 4

NDVI_2 = (NIR - Red) / (NIR + Red)

I am getting the following min and max values:

NDVI_2.min()
Out[756]: 0.0

NDVI_2.max()
Out[757]: 21845.0

If you have a look the the histogram of NDVI_2 enter image description here

In addition, looking for the pixels that have value > 1:

print("Values bigger than 1 =", NDVI_2[NDVI_2> 1])

Values bigger than 1 = [ 4.7842736   4.72178932  4.71683968 ... 20.07305095 15.2136194
     17.2828496 ]

print("Their indices are ", numpy.nonzero(NDVI_2 > 1))

Their indices are  (array([0, 0, 0, ..., 0, 0, 0]), array([    0,     0,     0, ..., 10979, 10979, 10979]), array([ 2382,  2383,  2384, ..., 10444, 10445, 10446]))

Any idea what I am doing wrong?

--- EDIT ---

Testing the masked. vs non-masked statistics:

NON MASKED

red = rasterio.open(S_files[2]).read()
nir = rasterio.open(S_files[3]).read()
ndvi = (nir - red) / (nir + red)
ndvi.max()
Out[111]: 21845.0
ndvi.min()
Out[112]: 0.0

MASKED

red = rasterio.open(S_files[2])
nir = rasterio.open(S_files[3])
redc, out_transform = rasterio.mask.mask(red, features, crop=True, all_touched=True)
nirc, out_transform = rasterio.mask.mask(nir, features, crop=True, all_touched=True)
ndvi = (nirc - redc) / (nirc + redc)
ndvi.max()
Out[118]: 63.337864077669906

ndvi.min()
Out[119]: 0.0

3 Answers 3

1

After using Python read the data, you may have ignored the no data value such as 0 or -999 from the numpy array.

When no data value gets involved in the calculation, it could affect the result but I think QGIS will ignore no data value while calculating.

7
  • When checking the min and max values of the original products, the min is 1 while max has 2456 for example. I thouhg that could be an issue but I do not have such values in the products. Any idea of a test code to check for that?
    – GCGM
    Commented Oct 7, 2019 at 7:32
  • @GCGM, Does your data have black areas which means no data value? Can you mask the no data value(using np.ma.masked_equal(array, 0)), then calculating NDVI again. Also, I checked your NDVI formula, I think it should be (NIR - RED) / (NIR + RED) instead of (red - nir ) / (red + nir). Your new edit used a wrong expression. In addition, I think your data is DN value which is all positive, so (NIR - RED) / (NIR + RED) should be range from -1 to 1. So check each band value read from rasterio first, or try gdal.Open().GetRasterBand().ReadAsArray() to read each band again.
    – Zac
    Commented Oct 7, 2019 at 12:16
  • I am using level 2 Sentinel2 products which are already corrected for atmospheric distortions and (not sure if you are familiar with Sentinel products) the pixel values should be already reflectance and not DN. Could it be then a matter of the data type being used (uint16 vs float)
    – GCGM
    Commented Oct 7, 2019 at 12:16
  • Just update the question with the right values and NDVI formula. According to earth.esa.int/web/sentinel/user-guides/sentinel-2-msi/…, level 2A products are atmospherically corrected reflectance images
    – GCGM
    Commented Oct 7, 2019 at 12:27
  • 1
    @GCGM, You need to divide your values by the QUATIFICATION value. Its value is 10000 and can be found in the metadata. Check this link , forum.step.esa.int/t/how-to-convert-dn-into-toa-reflectance/….
    – Zac
    Commented Oct 7, 2019 at 12:30
1

I got a similar weird result when I used complete jp2 images (with DN -digital numbers- values without transforming in reflectances) for calculating NDVI index for a Spanish region.

import rasterio
b8 = '/home/zeito/pyqgis_data/jp2/T30STF_20170422T110651_B08.jp2'
b4 = '/home/zeito/pyqgis_data/jp2/T30STF_20170422T110651_B04.jp2'
NIR = rasterio.open(b8).read() # band 8
Red = rasterio.open(b4).read() # band 4
NDVI_2 = (NIR - Red) / (NIR + Red)
print(NDVI_2.min())
0.0
print(NDVI_2.max())
53.7329498767461

However, when I masked them by using layer of following image, with format jp2 or geotiff, I got values as expected.

enter image description here

By using mask and jp2 images:

b8 = '/home/zeito/pyqgis_data/jp2/b8.jp2'
b4 = '/home/zeito/pyqgis_data/jp2/b4.jp2'
NIR = rasterio.open(b8).read() # band 8
Red = rasterio.open(b4).read() # band 4
NDVI_2 = (NIR - Red) / (NIR + Red)
print(NDVI_2.min())
0.1742668509135575
print(NDVI_2.max())
0.6112084063047285

By using mask and geotiff images:

b8 = '/home/zeito/pyqgis_data/jp2/b8.tif'
b4 = '/home/zeito/pyqgis_data/jp2/b4.tif'
NIR = rasterio.open(b8).read() # band 8
Red = rasterio.open(b4).read() # band 4
NDVI_2 = (NIR - Red) / (NIR + Red)
print(NDVI_2.min())
0.17488100721633656
print(NDVI_2.max())
0.6098130841121495

Concluding, I think it could be a memory issue due handling a big jp2 images or it is necessary to convert previously jp2 images in reflectance values multiplying by 0.0001 factor (because default values in these numpy arrays are uint16 dtype).

1
  • Thanks. I have tested the masked vs. non-masked (see question edit) and still getting more than -1 / +1 values. I think the issue is more with the unit16 dtype of the array. Specifying as.type('float') seems to solve the issue as now the values are wihin the expected NDVI ranges. Any suggestion on why using uint16 is causing the problem and float not?
    – GCGM
    Commented Oct 7, 2019 at 7:40
0

For those of you interested in the future, divinding by the quantification value as suggested by @Zac Wang produces the right result

Another option is to change the data type from uint16 to float. Using the following code, the issue is solved:

red = rasterio.open(S_files[2]).read().astype('float')
nir = rasterio.open(S_files[3]).read().astype('float')
ndvi = (nir - red) / (nir + red)
ndvi.max()
Out[111]: 0.85
ndvi.min()
Out[112]: -0.33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.