减去 L8 条带时 gdal_calc.py 中的负值饱和

Question

我正在尝试用 gdal_calc.py 减去 Landsat8 的两个波段，如下所示：

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

很明显输出是一个带符号的数字（即Int16），但是我从来没有得到负值，当一个值在减法后应该是负数时我得到值 32767。是否有在 --calc 表达式中解决这个问题的方法？这是计算 L8 的 NDVI 的第一步。

Answer 1

可以用numpy.ndarray的astype()方法解决，也就是calc表达式中A和B表示的数据类型.

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

这些int16实际上是numpy.int16但是计算表达式发生在numpy名字space.

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

如果您事先为输入栅格设置 NoDataValue 会更好。这将确保栅格大部分内的 calculated-zero NDVI 值不会被标记为 NoData，因为它们与 NoData space 中的 0 - 0 "NDVI" 具有相同的值.

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

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