1. main@rasterio.groups.io
Topics

Rasterio result different than gdal_calc

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

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)

You could convert this to a float32 now instead of doing it twice in a statement below.

 
with rasterio.open(band2Image) as dsNIR:
bandNIR = dsNIR.read(1)

Same comment.

 
# Allow division by zero
numpy.seterr(divide='ignore', invalid='ignore')

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.

ndvi = (bandNIR.astype(float)-bandRed.astype(float))/(bandNIR.astype(float) + bandRed.astype(float))

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

Like I said, we need to see the gdal_calc expression. I suspect nodata pixels are involved.

--
Sean Gillies
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
 
1 - 3 of 3