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

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.

-
One way to implement raster_diff(old_img, new_img) == "Identical" would be to check that the zonal max of the absolute value of the difference equals 0, where the zone is taken over the entire grid extent. Is this the sort of solution you're looking for? (If so, it would need to be refined to check that any NoData values are consistent, too.) –  whuber Aug 27 '14 at 19:54
@whuber thanks for ensuring proper NoData handling stays in the conversation. –  matt wilkie Aug 27 '14 at 21:59
if you can check that len(numpy.unique(yourraster)) == 2, then you know that it has 2 unique values and you can safely do this. –  RemcoGerlich Aug 28 '14 at 12:48
@Remco The algorithm underlying numpy.unique is going to be more computationally expensive (both in terms of time and space) than most other ways to check that the difference is a constant. When confronted with a difference between two very large floating point rasters that exhibit many differences (such as comparing an original to a lossy compressed version) it would likely bog down forever or fail completely. –  whuber Aug 28 '14 at 14:21
@Aaron, I got pulled off the project to do other things. Part of that was because the development time kept growing: too many edge cases to handle automatically, so the decision was made to throw the problem back at the people generating the images rather than fix them. (e.g. "Your disk quota is X. You learn how to work inside it.") However gdalcompare.py showed great promise (see answer) –  matt wilkie Mar 2 at 22:20

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)

d = numpy.array_equal(r1,r2)

if d == False:
print "They differ"

else:
print "They are the same"

-
That looks sweet and simple. I am curious about two details (which, technical though they are, could be crucial). First, does this solution handle NoData values correctly? Second, how does its speed compare with using built-in functions intended for grid comparisons, such as zonal summaries? –  whuber Aug 27 '14 at 19:56
Good points @whuber. I made a quick adjustment to the script that should take into account the shape and elements. I will check the points you brought up and report the findings. –  Aaron Aug 27 '14 at 20:05
@whuber Regarding the NoData handling, RasterToNumPyArray assigns by default the input raster's NoData value to the array. The user can specify a different value, although that would not apply in Matt's case. Regarding the speed, It took 4.5 seconds for the script to compare 2 4-bit rasters with 6210 columns and 7650 rows (DOQQ extent). I have not compared the method to any zonal summaries. –  Aaron Aug 27 '14 at 21:59
I folded in the gdal equivalent, adapted from gis.stackexchange.com/questions/32995/… –  matt wilkie May 4 at 16:24

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.

-
intriguing, I love gdal, didn't know about this script. Docs for interpreting the results are sparse to non-existent though ;-). In my initial testing it reports differences in colour interpretation and pallettes, meaning it might too specific for my current needs. I'm still exploring it though. (note: this answer too short to be a good fit here, link only answers are discouraged, please consider fleshing it out). –  matt wilkie Aug 27 '14 at 21:42

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.

-
Would this work with 32-bit floats? Would building and comparing two tables actually be any faster (or easier) than examining the values of the difference of the two rasters (which in principle should be just zero and NoData)? –  whuber Aug 27 '14 at 20:08
You are right that it would not work with 32-bit float and I did not check for the speed. However, building the attribute table needs to read the data only once and can help to avoid the 1-bit compression when you know that it will fail. Also I don't know the size of the images, but sometimes you cannot store them in memory. –  radouxju Aug 27 '14 at 20:18
@radouxju the images range up to 100,000px to a side, so ~30GB uncompressed. We don't have machine with that much ram (though perhaps with virtual...) –  matt wilkie Aug 27 '14 at 21:49
It seems likely RAM would not be an issue provided you stick with ArcGIS native operations. It's pretty good with RAM usage when processing grids: internally it can do the processing row-by-row, by groups of rows, and by rectangular windows. Local operations like subtracting one grid from another can operate essentially at the speed of input and output, requiring only one (relatively tiny) buffer for each input dataset. Constructing an attribute table requires an additional hash table--which would be minuscule when only one or two values show up, but could be enormous for arbitrary grids. –  whuber Aug 27 '14 at 22:39
numpy will do a lot of swapping with 2*30Go arrays, this is not ArcGIS anymore. I assumed based on the printscreen that the images are classified images (most with only afew values), so you don't expect so many classes. –  radouxju Aug 28 '14 at 7:09

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

-
One huge way to fool these statistics would be to permute the cell contents (which can happen, and does, when image dimensions are not quite right). On very large rasters neither the SD nor the mean would reliably detect a few small changes scattered about (especially if a few pixels were just dropped). Conceivably they would not detect a wholesale resampling of the grid, either, provided cubic convolution were used (which is intended to preserve the mean and SD). It would seem prudent instead to compare the SD of the difference of the grids to zero. –  whuber Aug 27 '14 at 21:25

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.

-
Subtraction seems like a good way to conduct a comparison. However, I believe the histogram would not be very useful in detecting problems with NoData values. Suppose, for instance, that the compression procedure eliminated a one-pixel border around the grid (this can happen!) but otherwise was accurate: all differences would still be zero. Also, did you notice that the OP needs to do this with 7000 raster data sets? I'm not sure he would relish examining 7000 plots. –  whuber Aug 27 '14 at 21:30