Растровые различия: как проверить, имеют ли изображения одинаковые значения?

10

Есть ли способ проверить, имеют ли два заданных растровых слоя одинаковое содержимое ?

У нас есть проблема с нашим корпоративным общим хранилищем: теперь он настолько большой, что для полного резервного копирования требуется более 3 дней. Предварительное расследование показало, что одним из самых больших занимающих пространство виновников являются растры включения / выключения, которые действительно должны храниться в виде 1-битных слоев со сжатием CCITT.

типичный настоящий / не присутствующий растр

Это примерное изображение в настоящее время 2-битное (так 3 возможных значения) и сохранено как сжатый формат LZW, 11 МБ в файловой системе. После преобразования в 1 бит (таким образом, 2 возможных значения) и применения сжатия CCITT Group 4 мы сократили его до 1,3 МБ, что является почти полным порядком экономии.

(Это на самом деле очень хорошо ведущий себя гражданин, другие хранятся как 32-битные числа с плавающей точкой!)

Это фантастические новости! Тем не менее, есть почти 7000 изображений, чтобы применить это тоже. Было бы просто написать скрипт для их сжатия:

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

... но в нем отсутствует жизненно важный тест: идентична ли вновь сжатая версия содержимому?

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

Существует ли инструмент или метод, который может автоматически (не) доказать, что содержимое Image-A идентично содержанию Image-B?

У меня есть доступ к ArcGIS 10.2 и QGIS, но я также открыт для всего остального, кроме того, что я могу избежать необходимости проверять все эти изображения вручную, чтобы убедиться в их корректности перед перезаписью. Было бы ужасно ошибочно преобразовать и перезаписать изображение , которое действительно было иметь больше , чем на / от значений в нем. Большинство собирают и генерируют тысячи долларов.

очень плохой результат

обновление: крупнейшие нарушители - 32-битные числа с плавающей запятой в диапазоне до 100 000 пикселей в стороне, то есть ~ 30 ГБ без сжатия.

Мэтт Уилки
источник
1
Одним из способов реализации raster_diff(old_img, new_img) == "Identical"было бы проверить, что зональный максимум абсолютного значения разности равен 0, где зона берется по всему экстенту сетки. Это то решение, которое вы ищете? (Если это так, его необходимо будет уточнить, чтобы убедиться, что любые значения NoData также являются согласованными.)
whuber
1
@whuber спасибо за обеспечение правильной NoDataобработки остается в разговоре.
Мэтт Уилки
если вы можете это проверить len(numpy.unique(yourraster)) == 2, то вы знаете, что он имеет 2 уникальных значения, и вы можете безопасно это сделать.
RemcoGerlich
@Remco Алгоритм, лежащий в основе алгоритма, numpy.uniqueбудет в вычислительном отношении более дорогим (как с точки зрения времени, так и пространства), чем большинство других способов проверить, является ли разница постоянной. Столкнувшись с разницей между двумя очень большими растрами с плавающей запятой, которые демонстрируют много различий (например, сравнение оригинала с сжатой версией с потерями), он, вероятно, навсегда потерпит неудачу или потерпит неудачу полностью.
whuber
1
@ Аарон, меня сняли с проекта, чтобы заняться другими делами. Частично это объяснялось тем, что время разработки продолжало расти: слишком много крайних случаев обрабатывалось автоматически, поэтому было принято решение бросить проблему людям, генерирующим изображения, а не исправлять их. (Например, «Ваша дисковая квота - X. Вы учитесь, как работать внутри нее».) Тем не менее, он gdalcompare.pyпоказал большое обещание ( см. ответ )
Matt Wilkie

Ответы:

8

Попробуйте преобразовать ваши растры в numy массивы, а затем проверьте, имеют ли они одинаковую форму и элементы с array_equal . Если они одинаковы, результат должен быть 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"
Аарон
источник
Это выглядит мило и просто. Мне любопытно две детали (которые, хотя они и технические, могут иметь решающее значение). Во-первых, правильно ли это решение обрабатывает значения NoData? Во-вторых, как его скорость сравнивается с использованием встроенных функций, предназначенных для сравнения сеток, таких как зональные сводки?
whuber
1
Хорошие очки @whuber. Я быстро настроил скрипт, который должен учитывать форму и элементы. Я проверю приведенные вами вопросы и сообщу о результатах.
Аарон
1
@whuber Что касается NoDataобработки, RasterToNumPyArrayпо умолчанию присваивает массиву значение NoData входного растра. Пользователь может указать другое значение, хотя это не относится к случаю Мэтта. Что касается скорости, сценарию потребовалось 4,5 секунды, чтобы сравнить 2 4-битных растра с 6210 столбцами и 7650 строками (экстент DOQQ). Я не сравнивал метод с какими-либо зональными аннотациями.
Аарон
1
Я сложил в эквивалент gdal, адаптированный с gis.stackexchange.com/questions/32995/…
Мэтт Вилки
4

Вы можете попробовать скрипт gdalcompare.py http://www.gdal.org/gdalcompare.html . Исходный код скрипта находится по адресу http://trac.osgeo.org/gdal/browser/trunk/gdal/swig/python/scripts/gdalcompare.py, и, поскольку это скрипт на python, он должен легко удалить ненужные тестируйте и добавляйте новые в соответствии с вашими текущими потребностями. Похоже, что скрипт выполняет попиксельное сравнение, считывая данные изображения из двух изображений по полосам, и это, вероятно, довольно быстрый и многократно используемый метод.

user30184
источник
1
Интересно, я люблю Гдал, не знал об этом сценарии. Документы для интерпретации результатов хотя и редки, либо вообще отсутствуют ;-). В моем первоначальном тестировании он сообщает о различиях в интерпретации цвета и палитрах, что означает, что он может быть слишком специфичным для моих текущих потребностей. Я все еще исследую это все же. (примечание: этот ответ слишком короткий, чтобы быть здесь подходящим, ответы только на ссылки не приветствуются, пожалуйста, подумайте об этом).
Мэтт Вилки
1

Я бы посоветовал вам создать свою таблицу атрибутов растра для каждого изображения, а затем сравнить их. Это не полная проверка (например, вычисление разницы между ними), но вероятность того, что ваши изображения отличаются с одинаковыми значениями гистограммы, очень и очень мала. Также он дает вам количество уникальных значений без NoData (из числа строк в таблице). Если ваш общий счет меньше размера изображения, вы знаете, что у вас есть пиксели NoData.

radouxju
источник
Будет ли это работать с 32-разрядными числами? Будет ли построение и сравнение двух таблиц на самом деле быстрее (или проще), чем проверка значений разности двух растров (которые в принципе должны равняться нулю и NoData)?
whuber
Вы правы, что он не будет работать с 32-битным поплавком, и я не проверял скорость. Однако построение таблицы атрибутов должно считывать данные только один раз и может помочь избежать 1-битного сжатия, когда вы знаете, что оно не удастся. Также я не знаю размер изображений, но иногда вы не можете сохранить их в памяти.
Radouxju
@radouxju изображения в диапазоне до 100 000 пикселей в сторону, поэтому ~ 30 ГБ без сжатия. У нас нет машины с таким большим количеством оперативной
памяти
Похоже, что оперативная память не будет проблемой, если вы будете придерживаться нативных операций ArcGIS. Это довольно хорошо с использованием оперативной памяти при обработке сеток: внутренне он может выполнять обработку построчно, по группам строк и по прямоугольным окнам. Локальные операции, такие как вычитание одной сетки из другой, могут работать по существу со скоростью ввода и вывода, для чего требуется только один (относительно небольшой) буфер для каждого входного набора данных. Для построения таблицы атрибутов требуется дополнительная хеш-таблица, которая может быть незначительной, когда отображаются только одно или два значения, но может быть огромной для произвольных сеток.
whuber
Numpy будет много менять с массивами 2 * 30Go, это больше не ArcGIS. Исходя из экрана печати, я предположил, что изображения являются классифицированными изображениями (большинство из которых имеют только несколько значений), поэтому вы не ожидаете так много классов.
Radouxju
0

Самое простое решение, которое я нашел, - это вычислить некоторую сводную статистику по растрам и сравнить ее. Я обычно использую стандартное отклонение и среднее значение, которые устойчивы к большинству изменений, хотя их можно обмануть, преднамеренно манипулируя данными.

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
источник
2
Один из огромных способов обмануть эту статистику - перестановка содержимого ячейки (что может случиться и происходит, когда размеры изображения не совсем верны). На очень больших растрах ни SD, ни среднее не могли бы надежно обнаружить несколько небольших изменений, разбросанных по всему (особенно, если несколько пикселей были просто отброшены). Вероятно, они также не будут обнаруживать оптовую повторную выборку сетки при условии использования кубической свертки (которая предназначена для сохранения среднего значения и SD). Вместо этого было бы разумно сравнить SD разницы сеток с нулем.
whuber
0

Самый простой способ - вычесть один растр из другого, если результат равен 0, то оба изображения одинаковы. Также вы можете увидеть гистограмму или график по цвету результата.

Pau
источник
Вычитание кажется хорошим способом провести сравнение. Однако я считаю, что гистограмма не будет очень полезна при обнаружении проблем со значениями NoData. Предположим, например, что процедура сжатия исключила однопиксельную границу вокруг сетки (это может произойти!), Но в противном случае была точна: все различия все равно были бы равны нулю. Кроме того, вы заметили, что OP должен делать это с 7000 наборов растровых данных? Я не уверен, что он с удовольствием изучит 7000 участков.
whuber