Search code examples
pythongdalsubtraction

saturation of negative values in gdal_calc.py when subtracting L8 bands

I am trying to subtract two bands of Landsat8 with gdal_calc.py as following:

gdal_calc.py -A LC8_B4.TIF -B LC8_B5.TIF --outfile="$pathout"/B5minusB4.tif --type='Int16' --calc="B-A"

It is clear that the output is a signed number (i.e. Int16), however I never get negative values, when a value is supposed to be negative after the subtraction I get the value 32767. Is there a way to fix this in the --calc expression? This is the first step in calculating the NDVI for L8.

Solution

  • You can solve this by using the astype() method of numpy.ndarray, which is the datatype represented by A and B in the calc expression.

    --calc="B.astype(int16) - A.astype(int16)"

    These int16 are actually numpy.int16 but the calculation expression takes place within the numpy namespace.

    gdal_calc.py -A LC8_B4.TIF -B LC8_B5.TIF \
  --outfile="$pathout"/B5minusB4.tif --type='Int16' \ 
  --calc="B.astype(int16) - A.astype(int16)"

    What's even better is if you set the NoDataValue for the input rasters beforehand. This will ensure the calculated-zero NDVI values within the good part of the raster aren't flagged as NoData because they have the same value as the 0 - 0 "NDVI" from the NoData space.

    gdal_edit.py -a_nodata 0 LC8_B4.TIF 
gdal_edit.py -a_nodata 0 LC8_B5.TIF 
# gdal_calc.py ...