[GIS] bad results NDVI in rasterio with jp2 sentinel-2 bands

jpeg 2000rasteriosentinel-2

I am using 2 Sentinel-2A bands to compute the NDVI. Here is my code:

import os, rasterio

outfile = r'some\path\ndvi.tif'
#url to the bands
b4 = 'http://sentinel-s2-l1c.s3.amazonaws.com/tiles/30/T/TK/2017/4/12/0/B04.jp2'
b8 = 'http://sentinel-s2-l1c.s3.amazonaws.com/tiles/30/T/TK/2017/4/12/0/B08.jp2'

#open the bands (I can't believe how easy is this with rasterio!)
with rasterio.open(b4) as red:
    RED = red.read()
with rasterio.open(b8) as nir:
    NIR = nir.read()

#compute the ndvi
ndvi = (NIR-RED)/(NIR+RED)
#print(ndvi.min(), ndvi.max()) The problem is alredy here

profile = red.meta
profile.update(driver='GTiff')
profile.update(dtype=rasterio.float32)

with rasterio.open(outfile, 'w', **profile) as dst:
    dst.write(ndvi.astype(rasterio.float32))

The problem is that, this is working fine, except for negative values. It seems that the num (NIR-RED) gets positives values no matter what. So vegetation and all the ndvi values between 0 and 1 keeps good values, but ndvi values between 0 and -1 gets weird values.
I can solve this saving RED and NIR to 2 GTiff rasters in local, but I guess that there must be one solution in order to avoid to do that.

Best Answer

That's because these bands come as unsigned integer 16 so the numpy division returns only positive integers.

You can replace ndvi = (NIR-RED)/(NIR+RED) by ndvi = (NIR.astype(float) - RED.astype(float)) / (NIR+RED)

Related Question