I try to subtract two .tif images. For that I use the following code:
import numpy as np
import os
from osgeo import gdal,ogr
import copy
from PIL import Image
import time
def loadTifAsArray(image, filepath):
print("Loading "+image)
path = os.path.join(filepath,image)
tif = gdal.Open(path)
tifArray = tif.ReadAsArray()
return tifArray
def subtract(image1, image2):
# Copy of one of the images is used for saving calculated values
print("Subtracting ...")
sub = copy.deepcopy(image1)
rows = len(sub)
cols = len(sub[0])
for px in range(cols):
for py in range(rows):
sub[px][py] = image1[px][py] - image2[px][py]
return sub
start_time = time.time()
cwd = os.getcwd()
filepath = os.path.join(cwd,'tifs')
arr = os.listdir(filepath)
tifList = []
for image in arr:
tifList.append(loadTifAsArray(image, filepath))
print("--- %s seconds for loading the images ---" % (time.time() - start_time))
sub = subtract(tifList[0], tifList[1])
print("--- %s seconds for loading and subtracting ---" % (time.time() - start_time))
I subtract the images, loaded as rasterdata, and simply make a deep copy of one of the images to store the calculated values in.
The problem is the calculated value. When I look at the values of both images at index [0][0], I get the following values:
print(image1[0][0])
print(image2[0][0])
505
549
When I try to subtract them I get this:
print(image1[0][0] - image2[0][0])
65492
I don't understand why that is and would appreciate any help!
That smells like overflow! Basically, I'm assuming your images are uint16 images, and negative numbers "wrap around" back to the maximal value.
Note that you expected 44, but got 2^16 - 44.
The work-around is pretty simple; cast your images to float32 for instance, by adding at the beginning of your subtract function:
image1 = np.array(image1).astype(np.float32)
image2 = np.array(image2).astype(np.float32)
Good luck!
P.S: Of course np is just my own personal preference. You can cast in other ways that suit you best.