2

I have a raster, input.tif. I can open it using rasterio and read 3 bands:

import rasterio
with rasterio.open('input.tif') as src:
    data1 = src.read(1)
    data2 = src.read(2)
    data3 = src.read(3)

Get the different values appearing in each band:

import numpy as np
print(np.unique(data1))
print(np.unique(data2))
print(np.unique(data3))

Outputs:

[181 201 217 222 230 237 255]
[120 156 186 219 245 250 255]
[145 156 176 179 191 196 199 255]

I want to add the values in each band to a new layer. I use gdal_calc to achieve this:

python3 gdal_calc.py -A input.tif -B input.tif -C input.tif  --A_band=1  --B_band=2 --C_band=3 --outfile=result.tif --calc="A+B+C" 

Read in the results:

with rasterio.open('result.tif') as src:
    data = src.read(1)

Print unique values:

print(np.unique(data))

Output:

[ 30 105 133 141 149 171 189 190 231 253]

I find this unexpected. If the values of the original bands were all more than a hundred (see above), how can the sum of bands have a value 30 for example?

What went wrong in the process above?

3
  • 2
    If you add values and fail to increase the bit depth you're likely to get integer overflow.
    – Vince
    Commented Jun 10, 2022 at 11:33
  • Oh ok. I did not find anything about bit depth at gdal_calc docs (gdal.org/programs/gdal_calc.html). Is there a simple way to increase bit depth for the output?
    – zabop
    Commented Jun 10, 2022 at 11:36
  • 2
    gdal.org/programs/gdal_calc.html: --type=<datatype>. Read also the notice Despite the datatype set using --type, when doing intermediate aritmethic operations using operands of the same type, the operation result will honor the original datatype. This may lead into unexpected results in the final result. The program is simple and it is really mostly up to user to avoid things that should not be done, like sum 181+120+145 into byte that has max at 255. I can also see an example about on-the-fly conversion of source data with --calc="(A.astype(numpy.float64).
    – user30184
    Commented Jun 10, 2022 at 11:51

1 Answer 1

3

Integer overflow in Numpy addition is a modulo operation and does not saturate and clip the values at the maximum datatype value, or does not change the datatype to accommodate the overflow.

try:

import numpy as np
a = np.uint8(250)
b = np.uint8(10)
print(a + b) # prints 4

Change the datatype to uint16, float32 or float64, depending on the use case.

Something along the lines of:

import rasterio
import numpy as np

with rasterio.open('input.tif') as src:
    data1 = src.read(1).astype(np.uint16) # float32 or float64 can also be used.
    data2 = src.read(2).astype(np.uint16)
    data3 = src.read(2).astype(np.uint16)

sum_data = data1 + data2 + data3
print(np.unique(sum_data)

for gdal_calc.py use --calc="(A.astype(numpy.float64) + (B.astype(numpy.float64) + (C.astype(numpy.float64)" (as per the comment from user30184).

1
  • The suggested gdal command fails with 0.. evaluation of calculation (A.astype(numpy.float64) + (B.astype(numpy.float64) + (C.astype(numpy.float64) failed unexpected EOF while parsing (<string>, line 1), but otherwise great answer, thanks.
    – zabop
    Commented Jun 10, 2022 at 14:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.