# Different NDVI in Python vs QGIS

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:

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`

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

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

``````red = rasterio.open(S_files[2])
nir = rasterio.open(S_files[3])
ndvi = (nirc - redc) / (nirc + redc)
ndvi.max()
Out[118]: 63.337864077669906

ndvi.min()
Out[119]: 0.0
``````

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.

• 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
– 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
• @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

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.

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

• 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

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