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