Rasterio result different than gdal_calc
anand@...
I am trying to recreate gdal_calc for simple calculation using below code.
import numpy
import rasterio
def merge_bands(calculationType, band1Image, band2Image, outputFilePath):
# Read raster bands directly to Numpy arrays.
#
# We handle the connections with "with"
with rasterio.open(band1Image) as dsRed:
bandRed = dsRed.read(1)
with rasterio.open(band2Image) as dsNIR:
bandNIR = dsNIR.read(1)
# Allow division by zero
numpy.seterr(divide='ignore', invalid='ignore')
ndvi = (bandNIR.astype(float)bandRed.astype(float))/(bandNIR.astype(float) + bandRed.astype(float))
kwargs = dsRed.meta
kwargs.update(dtype=rasterio.float32, count=1, compress='lzw')
with rasterio.open(outputFilePath, 'w', **kwargs) as dst:
dst.write_band(1, ndvi.astype(rasterio.float32))
But the result of this script is different than gdal_calc: gdal_cal result range 0.1  +0.09 rasterio result range 0.02  0.1


Sean Gillies
Hi, Can you also include the expressions that you're using with gdal_calc? Without them we can only guess at the cause of the difference in the results. I don't see any big problems with the way you're using Rasterio. I have a couple comments below: On Wed, Nov 21, 2018 at 3:35 PM anand via Groups.Io <anand=geospoc.com@groups.io> wrote: I am trying to recreate gdal_calc for simple calculation using below code. You could convert this to a float32 now instead of doing it twice in a statement below.
Same comment.
Unless you're doing this in gdal_calc too, this could be a cause of differences. I think it would be better to add keyword argument `masked=True` when you read the source bands. This will give you masked arrays that are safer to work when calculating the NDVI.
If you didn't specify `astype(float)` here, Numpy would convert to float64 as soon as it performed the division operation, and would return a float64 ndarray.
Like I said, we need to see the gdal_calc expression. I suspect nodata pixels are involved. Sean Gillies


anand@...
Hi Sean,
Sure, I will try those changes you have mentioned. here is the gdal_calc expression I used : gdal_calc.py A T1_B5.TIF A_band=1 B T1_B4.TIF B_band=1 outfile=ndvi_gal_1.TIF calc='(A.astype(float)B)/(A.astype(float)+B)' type='Float32' co='COMPRESS=LZW'

