Rasterio result different than gdal_calc


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

Join main@rasterio.groups.io to automatically receive all group messages.