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

10

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.

typowy obecny / nieobecny raster

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.

bardzo zły wynik

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

matowe wilkie
źródło
1
Jednym ze sposobów realizacji raster_diff(old_img, new_img) == "Identical"byłoby sprawdzenie, czy strefowe maksimum wartości bezwzględnej różnicy wynosi 0, gdzie strefa jest przejmowana przez cały zasięg siatki. Czy tego rodzaju rozwiązania szukasz? (Jeśli tak, należy dopracować, aby sprawdzić, czy wartości NoData również są spójne.)
whuber
1
@ Whuber dziękuję za zapewnienie właściwej NoDataobsługi pozostaje w rozmowie.
matt wilkie
jeśli możesz to sprawdzić len(numpy.unique(yourraster)) == 2, to wiesz, że ma 2 unikalne wartości i możesz to bezpiecznie zrobić.
RemcoGerlich
@Remco Podstawowy algorytm numpy.uniquebędzie droższy pod względem obliczeniowym (zarówno pod względem czasu, jak i przestrzeni) niż większość innych sposobów sprawdzania, czy różnica jest stała. W obliczu różnicy między dwoma bardzo dużymi zmiennoprzecinkowymi rastrami, które wykazują wiele różnic (takich jak porównanie oryginału do wersji ze stratną kompresją), prawdopodobnie zapadłby się na zawsze lub całkowicie zawiódł.
whuber
1
@Aaron, zostałem wyłączony z projektu, aby robić inne rzeczy. Częściowo było to spowodowane tym, że czas projektowania stale się wydłużał: zbyt wiele przypadków brzegowych, aby można było je obsługiwać automatycznie, dlatego postanowiono rzucić problem z powrotem osobom generującym obrazy zamiast je naprawić. (np. „Twój przydział dysku to X. Ty uczysz się, jak w nim pracować.”) gdalcompare.py
Okazało się

Odpowiedzi:

8

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"
Aaron
źródło
To wygląda słodko i prosto. Jestem ciekawy dwóch szczegółów (które, choć techniczne, mogą być kluczowe). Po pierwsze, czy to rozwiązanie poprawnie obsługuje wartości NoData? Po drugie, jak wypada jego prędkość w porównaniu z użyciem wbudowanych funkcji przeznaczonych do porównań siatki, takich jak streszczenia strefowe?
whuber
1
Dobre punkty @ whuber. Szybko dostosowałem scenariusz, który powinien uwzględniać kształt i elementy. Sprawdzę punkty, które poruszyłeś, i przedstawię wyniki.
Aaron
1
@ whuber Jeśli chodzi o NoDataobsługę, RasterToNumPyArraydomyślnie przypisuje tablicę wartość NoData wejściowego rastra. Użytkownik może podać inną wartość, chociaż nie miałoby to zastosowania w przypadku Matta. Jeśli chodzi o szybkość, skrypt zajął 4,5 sekundy, aby porównać 2 4-bitowe rastry z 6210 kolumnami i 7650 wierszami (zakres DOQQ). Nie porównałem tej metody do żadnych streszczeń strefowych.
Aaron
1
Złożyłem ekwiwalent gdal, zaadaptowany z gis.stackexchange.com/questions/32995/...
mat wilkie
4

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.

użytkownik30184
źródło
1
intrygujące, kocham gdal, nie wiedziałem o tym skrypcie. Dokumenty do interpretacji wyników są jednak rzadkie do nieistniejących ;-). Podczas moich pierwszych testów zgłasza różnice w interpretacji kolorów i paletach, co oznacza, że ​​może być zbyt specyficzne dla moich bieżących potrzeb. Nadal jednak to badam. (uwaga: ta odpowiedź jest zbyt krótka, aby można ją było tutaj dobrze dopasować, odradza się tylko odpowiedzi na linki, proszę rozważyć jej uzupełnienie).
matt wilkie
1

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.

radouxju
źródło
Czy działałoby to z 32-bitowymi zmiennoprzecinkowymi? Czy budowanie i porównywanie dwóch tabel byłoby w rzeczywistości szybsze (lub łatwiejsze) niż badanie wartości różnicy między dwoma rastrami (które w zasadzie powinny wynosić tylko zero i NoData)?
whuber
Masz rację, że nie będzie działać z 32-bitową liczbą zmiennoprzecinkową i nie sprawdziłem prędkości. Jednak zbudowanie tabeli atrybutów musi odczytać dane tylko raz i może pomóc uniknąć kompresji 1-bitowej, gdy wiadomo, że się nie powiedzie. Nie znam też rozmiaru obrazów, ale czasami nie można ich zapisać w pamięci.
radouxju
@radouxju obrazy mają rozmiar do 100 000 pikseli na bok, więc ~ 30 GB jest nieskompresowane. Nie mamy maszyny z tak dużą
ilością pamięci
Wygląda na to, że pamięć RAM nie byłaby problemem, pod warunkiem, że będziesz trzymać się natywnych operacji ArcGIS. Jest całkiem dobry w użyciu pamięci RAM podczas przetwarzania siatek: wewnętrznie może przetwarzać przetwarzanie rząd po rzędzie, według grup rzędów i prostokątnych okien. Lokalne operacje, takie jak odejmowanie jednej siatki od drugiej, mogą działać zasadniczo z prędkością wejścia i wyjścia, wymagając tylko jednego (względnie małego) bufora dla każdego wejściowego zestawu danych. Konstruowanie tabeli atrybutów wymaga dodatkowej tabeli skrótów - która byłaby niewielka, gdy pojawi się tylko jedna lub dwie wartości, ale może być ogromna dla dowolnych siatek.
whuber
numpy będzie dużo zamieniać z tablicami 2 * 30Go, to już nie jest ArcGIS. W oparciu o ekran drukowania założyłem, że obrazy są obrazami klasyfikowanymi (większość z tylko kilkoma wartościami), więc nie oczekujesz tylu klas.
radouxju
0

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.")
scw
źródło
2
Jednym wielkim sposobem na oszukanie tych statystyk byłoby permutowanie zawartości komórki (co może się zdarzyć i ma miejsce, gdy wymiary obrazu nie są całkiem właściwe). Na bardzo dużych rastrach ani SD, ani średnia nie wykryłyby pewnie kilku małych zmian rozproszonych wokół (szczególnie, gdyby kilka pikseli zostało po prostu upuszczonych). Możliwe, że nie wykryliby również hurtowego ponownego próbkowania siatki, pod warunkiem, że zastosowano splot sześcienny (który ma na celu zachowanie średniej i SD). Rozsądnym wydaje się zamiast tego porównanie SD różnicy siatek do zera.
whuber
0

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.

Pau
źródło
Odejmowanie wydaje się dobrym sposobem na porównanie. Uważam jednak, że histogram nie byłby bardzo przydatny w wykrywaniu problemów z wartościami NoData. Załóżmy na przykład, że procedura kompresji wyeliminowała jednopikselową ramkę wokół siatki (może się to zdarzyć!), Ale poza tym była dokładna: wszystkie różnice byłyby nadal zerowe. Czy zauważyłeś również, że PO musi to zrobić z 7000 zestawami danych rastrowych? Nie jestem pewien, czy podoba mu się zbadanie 7000 fabuł.
whuber