Dif raster: como verificar se as imagens têm valores idênticos?

10

Existe um meio de verificar se duas camadas de varredura fornecidas têm conteúdo idêntico ?

Temos um problema no nosso volume de armazenamento compartilhado corporativo: agora é tão grande que leva mais de três dias para realizar um backup completo. A investigação preliminar revela que um dos maiores culpados por consumir espaço são os controladores on / off que realmente devem ser armazenados como camadas de 1 bit com compressão CCITT.

uma varredura típica presente / não presente

Atualmente, esta imagem de amostra é de 2 bits (portanto, 3 valores possíveis) e salva como tiff compactado LZW, 11 MB no sistema de arquivos. Depois de converter para 1bit (portanto, 2 valores possíveis) e aplicar a compactação do CCITT Group 4, reduzimos para 1,3 MB, quase uma ordem completa de magnitude de economia.

(Na verdade, este é um cidadão muito bem-comportado, existem outros armazenados como bóia de 32 bits!)

Esta é uma notícia fantástica! No entanto, existem quase 7.000 imagens para aplicar isso também. Seria simples escrever um script para compactá-los:

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

... mas está faltando um teste vital: a versão recém-compactada é idêntica ao conteúdo?

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

Existe uma ferramenta ou método que pode (des) provar automaticamente que o conteúdo da imagem A é idêntico ao valor da imagem b?

Eu tenho acesso ao ArcGIS 10.2 e QGIS, mas também estou aberto a quase tudo o que pode evitar a necessidade de inspecionar todas essas imagens manualmente para garantir a correção antes da substituição. Seria horrível converter e substituir por engano uma imagem que realmente tinha mais do que valores de ativação / desativação. A maioria custa milhares de dólares para reunir e gerar.

um resultado muito ruim

atualização: Os maiores infratores são carros alegóricos de 32 bits que variam de até 100.000px para um lado, portanto, ~ 30 GB não compactados.

python raster compression differences Matt Wilson
fonte
1
Uma maneira de implementar raster_diff(old_img, new_img) == "Identical"seria verificar se o máximo zonal do valor absoluto da diferença é igual a 0, onde a zona é ocupada em toda a extensão da grade. Esse é o tipo de solução que você está procurando? (Em caso afirmativo, seria necessário refinar para verificar se quaisquer valores de NoData também são consistentes.)
whuber
1
A @whuber agradece por garantir o NoDatamanuseio adequado durante a conversa.
precisa saber é o seguinte
se você pode verificar isso len(numpy.unique(yourraster)) == 2, sabe que ele possui 2 valores exclusivos e pode fazer isso com segurança.
RemcoGerlich 28/08
@Remco O algoritmo subjacente numpy.uniqueserá mais caro em termos de computação (tanto em termos de tempo quanto de espaço) do que a maioria das outras maneiras de verificar se a diferença é constante. Quando confrontado com uma diferença entre dois rasters de ponto flutuante muito grandes que exibem muitas diferenças (como comparar um original a uma versão compactada com perda), ele provavelmente fica parado para sempre ou falha completamente.
whuber
1
@ Aaron, fui retirado do projeto para fazer outras coisas. Parte disso ocorreu porque o tempo de desenvolvimento continuou crescendo: muitos casos extremos para lidar automaticamente, por isso foi tomada a decisão de devolver o problema às pessoas que geravam as imagens, em vez de corrigi-las. (por exemplo, "Sua cota de disco é X. Você aprender a trabalhar dentro dela.") No entanto gdalcompare.pymostrou uma grande promessa ( ver resposta )
Matt Wilkie

Respostas:

8

Tente converter suas rasters em matrizes numpy e verifique se elas têm a mesma forma e elementos com array_equal . Se forem iguais, o resultado deve 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
fonte
Parece doce e simples. Estou curioso sobre dois detalhes (que, por mais técnicos que sejam, podem ser cruciais). Primeiro, esta solução lida com os valores NoData corretamente? Segundo, como sua velocidade se compara ao uso de funções internas destinadas a comparações de grade, como resumos de zonas?
whuber
1
Bons pontos @whuber. Fiz um rápido ajuste no script que deveria levar em conta a forma e os elementos. Vou verificar os pontos que você levantou e relatar as descobertas.
Aaron
1
@whuber Em relação à NoDatamanipulação, RasterToNumPyArrayatribui por padrão o valor NoData da varredura de entrada à matriz. O usuário pode especificar um valor diferente, embora isso não se aplique no caso de Matt. Em relação à velocidade, o script levou 4,5 segundos para comparar 2 rasters de 4 bits com 6210 colunas e 7650 linhas (extensão DOQQ). Não comparei o método a nenhum resumo de zonas.
Aaron
1
I dobrado no equivalente gdal, adaptado de gis.stackexchange.com/questions/32995/...
Wilkie mate
4

Você pode tentar o script gdalcompare.py http://www.gdal.org/gdalcompare.html . O código fonte do script está em http://trac.osgeo.org/gdal/browser/trunk/gdal/swig/python/scripts/gdalcompare.py e, por ser um script python, deve ser fácil remover os itens desnecessários. testes e adicione novos para atender às suas necessidades atuais. O script parece fazer comparação pixel por pixel lendo dados de imagem das duas imagens, banda por banda, e esse é provavelmente um método rápido e reutilizável.

user30184
fonte
1
intrigante, eu amo gdal, não sabia sobre esse script. Documentos para interpretar os resultados são escassos ou inexistentes ;-). Nos meus testes iniciais, ele relata diferenças na interpretação das cores e paletas, o que significa que pode ser específico demais para minhas necessidades atuais. Ainda estou explorando isso. (nota: esta resposta é muito curta para ser um bom ajuste aqui; as respostas somente para links são desencorajadas; considere aprofundar).
Matt Wilkie
1

Sugiro que você construa sua tabela de atributos raster para cada imagem e compare as tabelas. Essa não é uma verificação completa (como calcular a diferença entre as duas), mas a probabilidade de suas imagens serem diferentes com os mesmos valores de histograma é muito pequena. Além disso, fornece o número de valores exclusivos sem NoData (do número de linhas na tabela). Se sua contagem total for menor que o tamanho da imagem, você saberá que possui pixels NoData.

radouxju
fonte
Isso funcionaria com flutuadores de 32 bits? Criar e comparar duas tabelas seria realmente mais rápido (ou mais fácil) do que examinar os valores da diferença das duas rasters (que em princípio deveriam ser apenas zero e NoData)?
whuber
Você está certo de que não funcionaria com flutuador de 32 bits e eu não verifiquei a velocidade. No entanto, a construção da tabela de atributos precisa ler os dados apenas uma vez e pode ajudar a evitar a compactação de 1 bit quando você souber que ela falhará. Também não sei o tamanho das imagens, mas às vezes você não pode armazená-las na memória.
Radouxju
@radouxju as imagens variam até 100.000px para um lado, então ~ 30GB não são compactados. Não temos máquina com tanta RAM (embora talvez com Virtual ...)
Matt Wilkie
Parece que a RAM não seria um problema, desde que você se atenha às operações nativas do ArcGIS. É muito bom com o uso de RAM ao processar grades: internamente, ele pode fazer o processamento linha por linha, por grupos de linhas e por janelas retangulares. Operações locais, como subtrair uma grade de outra, podem operar essencialmente na velocidade de entrada e saída, exigindo apenas um buffer (relativamente pequeno) para cada conjunto de dados de entrada. A construção de uma tabela de atributos requer uma tabela de hash adicional - que seria minúscula quando apenas um ou dois valores aparecessem, mas poderia ser enorme para grades arbitrárias.
whuber
O numpy fará muitas trocas com matrizes 2 * 30Go, isso não é mais o ArcGIS. Eu assumi com base na tela de impressão que as imagens são classificadas (a maioria com apenas alguns valores), então você não espera tantas classes.
Radouxju 28/08
0

A solução mais simples que encontrei é calcular algumas estatísticas resumidas sobre os rasters e compará-las. Eu geralmente uso desvio e média padrão, que são robustos para a maioria das alterações, embora seja possível enganá-las manipulando intencionalmente os dados.

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
fonte
2
Uma maneira imensa de enganar essas estatísticas seria permutar o conteúdo da célula (o que pode acontecer, e ocorre, quando as dimensões da imagem não são muito corretas). Em rasters muito grandes, nem o SD nem a média detectariam com segurança algumas pequenas alterações espalhadas (especialmente se apenas alguns pixels fossem descartados). É concebível que eles também não detectariam uma reamostragem por atacado da grade, desde que fosse usada convolução cúbica (que visa preservar a média e o desvio padrão). Parece prudente comparar o desvio padrão da diferença das grades com zero.
whuber
0

A maneira mais fácil é subtrair uma varredura da outra, se o resultado for 0, as duas imagens serão iguais. Além disso, você pode ver o histograma ou plotar pela cor do resultado.

Pau
fonte
Subtração parece ser uma boa maneira de realizar uma comparação. No entanto, acredito que o histograma não seria muito útil na detecção de problemas com os valores NoData. Suponha, por exemplo, que o procedimento de compactação tenha eliminado uma borda de um pixel ao redor da grade (isso pode acontecer!), Mas, caso contrário, seria preciso: todas as diferenças ainda seriam nulas. Além disso, você notou que o OP precisa fazer isso com 7000 conjuntos de dados rasterizados? Não sei se ele gostaria de examinar 7.000 parcelas.
whuber