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

10

¿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.

un típico raster presente / no presente

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.

un muy mal resultado

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.

wilkie mate
fuente
1
Una forma de implementar raster_diff(old_img, new_img) == "Identical"sería verificar que el máximo zonal del valor absoluto de la diferencia sea igual a 0, donde la zona se toma en toda la extensión de la cuadrícula. ¿Es este el tipo de solución que estás buscando? (Si es así, tendría que ser refinado para comprobar que todos los valores NoData son consistentes, también.)
whuber
1
@whuber gracias por garantizar que el NoDatamanejo adecuado se mantenga en la conversación.
Matt Wilkie
si puede verificar eso len(numpy.unique(yourraster)) == 2, entonces sabe que tiene 2 valores únicos y puede hacerlo de manera segura.
RemcoGerlich
@Remco El algoritmo subyacente numpy.uniqueserá más costoso desde el punto de vista computacional (tanto en términos de tiempo como de espacio) que la mayoría de las otras formas de verificar que la diferencia es una constante. Cuando se enfrenta a una diferencia entre dos rásteres de punto flotante muy grandes que exhiben muchas diferencias (como comparar un original con una versión comprimida con pérdida), es probable que se empantane para siempre o falle por completo.
whuber
1
@ Aaron, me sacaron del proyecto para hacer otras cosas. Parte de eso se debió a que el tiempo de desarrollo siguió creciendo: demasiados casos extremos para manejarlos automáticamente, por lo que se tomó la decisión de devolver el problema a las personas que generan las imágenes en lugar de solucionarlos. (por ejemplo, "Su cuota de disco es X. Aprende a trabajar dentro de ella"). Sin embargo, gdalcompare.pymostró una gran promesa ( ver respuesta )
matt wilkie

Respuestas:

8

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"
Aaron
fuente
Eso se ve dulce y simple. Tengo curiosidad acerca de dos detalles (que, por técnicos que sean, podrían ser cruciales). Primero, ¿esta solución maneja los valores NoData correctamente? En segundo lugar, ¿cómo se compara su velocidad con el uso de funciones integradas destinadas a las comparaciones de cuadrícula, como los resúmenes zonales?
whuber
1
Buenos puntos @whuber. Hice un ajuste rápido al guión que debería tener en cuenta la forma y los elementos. Revisaré los puntos que mencionaste e informaré los resultados.
Aaron
1
@whuber En cuanto al NoDatamanejo, RasterToNumPyArrayasigna por defecto el valor NoData del ráster de entrada a la matriz. El usuario puede especificar un valor diferente, aunque eso no se aplicaría en el caso de Matt. En cuanto a la velocidad, el script tardó 4,5 segundos en comparar 2 rásteres de 4 bits con 6210 columnas y 7650 filas (extensión DOQQ). No he comparado el método con ningún resumen zonal.
Aaron
1
Doblé el equivalente en gdal, adaptado de gis.stackexchange.com/questions/32995/…
matt wilkie
4

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.

usuario30184
fuente
1
intrigante, me encanta gdal, no sabía sobre este guión. Sin embargo, los documentos para interpretar los resultados son escasos o inexistentes ;-). En mi prueba inicial, informa diferencias en la interpretación del color y las paletas, lo que significa que podría ser demasiado específico para mis necesidades actuales. Aunque todavía lo estoy explorando. (nota: esta respuesta es demasiado corta para encajar bien aquí, se desaconsejan las respuestas de solo enlace, considere la posibilidad de completarla).
Matt Wilkie
1

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.

radouxju
fuente
¿Funcionaría esto con flotadores de 32 bits? ¿Construir y comparar dos tablas en realidad sería más rápido (o más fácil) que examinar los valores de la diferencia de los dos rásteres (que en principio deberían ser solo cero y NoData)?
whuber
Tienes razón en que no funcionaría con flotante de 32 bits y no verifiqué la velocidad. Sin embargo, la creación de la tabla de atributos debe leer los datos solo una vez y puede ayudar a evitar la compresión de 1 bit cuando se sabe que fallará. Además, no sé el tamaño de las imágenes, pero a veces no puedes almacenarlas en la memoria.
radouxju
@radouxju las imágenes varían hasta 100,000 px a un lado, por lo que ~ 30GB sin comprimir. No tenemos máquina con tanto carnero (aunque quizás con virtual ...)
matt wilkie
Parece probable que la RAM no sea un problema siempre que se mantenga con las operaciones nativas de ArcGIS. Es bastante bueno con el uso de RAM al procesar cuadrículas: internamente puede hacer el procesamiento fila por fila, por grupos de filas y por ventanas rectangulares. Las operaciones locales, como restar una cuadrícula de otra, pueden operar esencialmente a la velocidad de entrada y salida, y requieren solo un búfer (relativamente pequeño) para cada conjunto de datos de entrada. La construcción de una tabla de atributos requiere una tabla hash adicional, que sería minúscula cuando solo se muestran uno o dos valores, pero podría ser enorme para cuadrículas arbitrarias.
whuber
numpy hará muchos intercambios con matrices 2 * 30Go, esto ya no es ArcGIS. Supuse, en base a la pantalla de impresión, que las imágenes son imágenes clasificadas (la mayoría con solo unos pocos valores), por lo que no espera tantas clases.
radouxju
0

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.")
scw
fuente
2
Una gran manera de engañar a estas estadísticas sería permutar el contenido de la celda (lo que puede suceder, y ocurre, cuando las dimensiones de la imagen no son del todo correctas). En rásteres muy grandes, ni la SD ni la media detectarían de manera confiable algunos pequeños cambios dispersos (especialmente si solo se eliminaran algunos píxeles). Posiblemente tampoco detectarían un remuestreo mayorista de la cuadrícula, siempre que se utilizara una convolución cúbica (que tiene como objetivo preservar la media y la DE). Parecería prudente en cambio comparar el SD de la diferencia de las cuadrículas a cero.
whuber
0

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.

Pau
fuente
La resta parece ser una buena forma de realizar una comparación. Sin embargo, creo que el histograma no sería muy útil para detectar problemas con los valores NoData. Supongamos, por ejemplo, que el procedimiento de compresión eliminó un borde de un píxel alrededor de la cuadrícula (¡esto puede suceder!) Pero por lo demás fue exacto: todas las diferencias seguirían siendo cero. Además, ¿notó que el OP necesita hacer esto con 7000 conjuntos de datos ráster? No estoy seguro de que le gustaría examinar 7000 parcelas.
whuber