Reputation: 515
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!
Upvotes: 0
Views: 439
Reputation: 2905
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.
Upvotes: 2