Raster diff: ¿cómo verificar si las imágenes tienen valores idénticos?

¿Hay algún medio para verificar si 2 capas ráster dadas tienen contenido idéntico ?

Tenemos un problema con nuestro volumen de almacenamiento compartido corporativo: ahora es tan grande que lleva más de 3 días realizar una copia de seguridad completa. La investigación preliminar revela que uno de los principales culpables que consumen espacio son los rásteres de encendido / apagado que realmente deberían almacenarse como capas de 1 bit con compresión CCITT.

Esta imagen de muestra es actualmente de 2 bits (por lo tanto, 3 valores posibles) y se guarda como tiff comprimido LZW, 11 MB en el sistema de archivos. Después de convertir a 1 bit (es decir, 2 valores posibles) y aplicar la compresión CCITT Group 4, lo reducimos a 1.3 MB, casi un orden completo de magnitud de ahorro.

(Este es en realidad un ciudadano muy bien educado, ¡hay otros almacenados como flotantes de 32 bits!)

¡Estas son noticias fantásticas! Sin embargo, hay casi 7,000 imágenes para aplicar esto también. Sería sencillo escribir un script para comprimirlos:

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

... pero le falta una prueba vital: ¿la versión recién comprimida es idéntica al contenido?

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

¿Existe alguna herramienta o método que pueda (des) probar automáticamente que el contenido de la Imagen-A es de valor idéntico al contenido de la Imagen-B?

Tengo acceso a ArcGIS 10.2 y QGIS, pero también estoy abierto a casi cualquier otra cosa que pueda evitar la necesidad de inspeccionar todas estas imágenes manualmente para garantizar la corrección antes de sobrescribir. Sería horrible para convertir por error y sobrescribir una imagen que realmente hizo tiene más de encendido / apagado valores en ella. La mayoría cuesta miles de dólares para reunir y generar.

Actualización: los mayores infractores son flotadores de 32 bits que varían hasta 100,000 px por lado, por lo que ~ 30 GB sin comprimir.

Intente convertir sus rásteres en matrices numpy y luego verifique si tienen la misma forma y elementos con array_equal . Si son iguales, el resultado debería ser 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"

Puede probar con el script gdalcompare.py http://www.gdal.org/gdalcompare.html . El código fuente del script está en http://trac.osgeo.org/gdal/browser/trunk/gdal/swig/python/scripts/gdalcompare.py y, dado que es un script de Python, debería ser fácil eliminar lo innecesario pruebas y agregue nuevas para satisfacer sus necesidades actuales. El guión parece hacer una comparación píxel por píxel al leer los datos de imagen de las dos imágenes banda por banda y ese es probablemente un método bastante rápido y reutilizable.

Sugeriría que construya su tabla de atributos ráster para cada imagen, luego puede comparar las tablas. Esta no es una verificación completa (como calcular la diferencia entre los dos), pero la probabilidad de que sus imágenes sean diferentes con los mismos valores de histograma es muy muy pequeña. También le da el número de valores únicos sin NoData (del número de filas en la tabla). Si su recuento total es menor que el tamaño de la imagen, sabe que tiene píxeles NoData.

La solución más simple que he encontrado es calcular algunas estadísticas de resumen en los rásteres y compararlos. Usualmente uso la desviación estándar y la media, que son robustas para la mayoría de los cambios, aunque es posible engañarlos manipulando intencionalmente los datos.

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

La forma más fácil es restar un ráster del otro, si el resultado es 0, ambas imágenes son iguales. También puede ver el histograma o trazar por color el resultado.