Raster diff: jak sprawdzić, czy obrazy mają identyczne wartości?

Czy istnieje sposób, aby sprawdzić, czy jakieś 2 podane warstwy rastrowe mają identyczną zawartość ?

Mamy problem z wolumenem naszej wspólnej pamięci masowej w firmie: jest on teraz tak duży, że wykonanie pełnej kopii zapasowej zajmuje ponad 3 dni. Wstępne dochodzenie ujawnia, że ​​jednym z największych zajmujących miejsce sprawców są rastry on / off, które naprawdę powinny być przechowywane jako warstwy 1-bitowe z kompresją CCITT.

Ten przykładowy obraz jest obecnie 2-bitowy (czyli 3 możliwe wartości) i zapisany jako skompresowany plik tiff LZW, 11 MB w systemie plików. Po konwersji na 1bit (czyli 2 możliwe wartości) i zastosowaniu kompresji CCITT Group 4, zmniejszamy ją do 1,3 MB, czyli prawie pełnego rzędu wielkości oszczędności.

(W rzeczywistości jest to bardzo dobrze wychowany obywatel, istnieją inne przechowywane jako 32-bitowe float!)

To fantastyczna wiadomość! Istnieje jednak prawie 7 000 zdjęć, które można zastosować. Łatwo byłoby napisać skrypt do ich skompresowania:

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

... ale brakuje istotnego testu: czy nowo skompresowana wersja jest identyczna z treścią?

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

Czy istnieje narzędzie lub metoda, która może (nie) automatycznie udowodnić, że zawartość obrazu A jest identyczna z treścią obrazu B?

Mam dostęp do ArcGIS 10.2 i QGIS, ale jestem również otwarty na większość innych rzeczy, niż jest w stanie uniknąć konieczności ręcznego sprawdzania wszystkich tych obrazów w celu zapewnienia poprawności przed nadpisaniem. Byłoby straszne błędnie konwertować i nadpisać obraz, który naprawdę nie mają więcej niż on / off wartości w nim. Większość kosztuje tysiące dolarów na zebranie i wygenerowanie.

aktualizacja: największymi przestępcami są 32-bitowe zmiennoprzecinkowe o zasięgu do 100 000 pikseli na bok, więc ~ 30 GB bez kompresji.

Spróbuj przekonwertować swoje rastry na tablice numpy, a następnie sprawdź, czy mają one ten sam kształt i elementy z array_equal . Jeśli są takie same, wynikiem powinno być 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"

Możesz spróbować ze skryptem gdalcompare.py http://www.gdal.org/gdalcompare.html . Kod źródłowy skryptu znajduje się na stronie http://trac.osgeo.org/gdal/browser/trunk/gdal/swig/python/scripts/gdalcompare.py, a ponieważ jest to skrypt w języku Python, powinno być łatwo usunąć niepotrzebne testuj i dodawaj nowe, aby dopasować je do bieżących potrzeb. Skrypt wydaje się porównywać piksel po pikselu, odczytując dane obrazu z dwóch obrazów pasmo po paśmie, i jest to prawdopodobnie dość szybka i wielokrotnego użytku metoda.

Sugerowałbym zbudowanie tabeli atrybutów rastrowych dla każdego obrazu, a następnie porównanie tabel. Nie jest to pełna kontrola (jak obliczenie różnicy między nimi), ale prawdopodobieństwo, że twoje obrazy różnią się przy tych samych wartościach histogramu, jest bardzo małe. Daje także liczbę unikalnych wartości bez NoData (od liczby wierszy w tabeli). Jeśli łączna liczba jest mniejsza niż rozmiar obrazu, wiesz, że masz piksele NoData.

Najprostszym rozwiązaniem, jakie znalazłem, jest obliczenie statystyk podsumowujących rastrów i ich porównanie. Zwykle używam odchyleń standardowych i średnich, które są odporne na większość zmian, chociaż można je oszukać, celowo manipulując danymi.

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

Najprostszym sposobem jest odjęcie jednego rastra od drugiego, jeśli wynikiem jest 0, oba obrazy są takie same. Możesz także zobaczyć histogram lub wykres według koloru wyniku.