Raster diff: how to check if images have identical values?

Question

Is there a means to check to see if any 2 given raster layers have identical content?

We have a problem on our corporate shared storage volume: it's now so big that it takes over 3 days to conduct a full back up. Preliminary investigation reveals one of the biggest space consuming culprits are on/off rasters that really should be stored as 1-bit layers with CCITT compression.

This sample image is currently 2bit (so 3 possible values) and saved as LZW compressed tiff, 11 MB in file system. After converting to 1bit (so 2 possible values) and applying CCITT Group 4 compression we get it down to 1.3 MB, almost a full order of magnitude of savings.

(This is actually a very well behaved citizen, there are others stored as 32bit float!)

This is fantastic news! However there are almost 7,000 images to apply this too. It would be straightforward to write a script to compress them:

for old_img in [list of images]:
    convert_to_1bit_and_compress(old_img)
    remove(old_img)
    replace_with_new(old_img, new_img)

...but it's missing a vital test: is the newly compressed version content-identical?

  if raster_diff(old_img, new_img) == "Identical":
      remove(old_img)
      rename(new_img, old_img)

Is there a tool or method which can automatically (dis)prove the content of Image-A is value-identical to the content of Image-B?

I have access to ArcGIS 10.2 and QGIS, but am also open to most anything else than can obviate the need to inspect all these images manually to ensure correctness before overwriting. It would be horrible to mistakenly convert and overwrite an image that really did have more than on/off values in it. Most cost thousands of dollars to gather and generate.

update: The biggest offenders are 32bit floats that range up to 100,000px to a side, so ~30GB uncompressed.

Answer 1

Try converting your rasters to numpy arrays and then check to see if they have the same shape and elements with array_equal. If they are the same, the result should be True:

ArcGIS:

import arcpy, numpy

raster1 = r'C:\path\to\raster.tif'
raster2 = r'C:\path\to\raster.tif'

r1 = arcpy.RasterToNumPyArray(raster1)
r2 = arcpy.RasterToNumPyArray(raster2)

d = numpy.array_equal(r1,r2)

if d == False:
    print "They differ"

else:
    print "They are the same"

GDAL:

import numpy
from osgeo import gdal        
    
raster1 = r'C:\path\to\raster.tif'
raster2 = r'C:\path\to\raster.tif'

ds1 = gdal.Open(raster1)
ds2 = gdal.Open(raster2)

r1 = numpy.array(ds1.ReadAsArray())
r2 = numpy.array(ds2.ReadAsArray())

d = numpy.array_equal(r1,r2)

if d == False:
    print "They differ"

else:
    print "They are the same"

Answer 2

You could have a try with gdalcompare.py script http://www.gdal.org/gdalcompare.html. The source code of the script is at http://trac.osgeo.org/gdal/browser/trunk/gdal/swig/python/scripts/gdalcompare.py and because it is a python script it should be easy to remove the unnecessary tests and add new ones to suit your current needs. The script seems to do pixel by pixel comparison by reading image data from the two images band by band and that is probably quite a fast and reusable method.

Answer 3

I would suggest that you build your raster attribute table for each image, then you can compare the tables. This is not a complete check (like computing the difference between the two), but the probability that your images are different with the same histogram values is very very small. Also it gives you the number of unique values without NoData (from the number of rows in the table). If your total count is less than the image size, you know that you have NoData pixels.

Answer 4

The simplest solution I've found is to compute some summary statistics on the rasters, and compare those. I usually use standard deviation and mean, which are robust to most changes, though it is possible to fool them by intentionally manipulating the data.

mean_obj = arcpy.GetRasterProperties(input_raster, 'MEAN')
mean = float(mean_obj.getOutput(0))
if round(mean, 4) != 0.2010:
    print("raster differs from expected mean.")

std_obj = arcpy.GetRasterProperties(input_raster, 'STD')
std = float(std_obj.getOutput(0))
if round(std, 4) != 0.0161:
    print("raster differs from expected standard deviation.")

Answer 5

The easiest way is to subtract one raster from the other, if the result is 0, then both images are the same. Also you can see the histogram or plot by color the result.