What could go wrong using gdal_calc to sum pixel values of a multiband raster?

Question

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?

Answer 1

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