Comparando duas geometrias no ArcPy?

18

Eu estou tentando comparar duas classes de recurso separadas para identificar diferenças entre elas (tipo de uma função diff). Meu fluxo de trabalho básico:

  1. Extraio as geometrias usando um SearchCursor
  2. Salve as geometrias das duas classes de recursos como GeoJSON usando um modificado __geo_interface__( obtenha -o de valveLondonreturn {'type': 'Polygon', 'coordinates': [[((pt.X, pt.Y) if pt else None) for pt in part] for part in self]} ). Isso é para evitar o objeto de geometria compartilhada que a ESRI usa com cursores e a incapacidade de fazer cópias profundas (algumas discussões aqui no gis.stackexchange falam sobre isso).
  3. Verifique as geometrias das duas classes de recursos com base em um identificador exclusivo. Por exemplo, compare a geometria FC1 OID1 com a geometria FC2 OID1. Para obter a geometria como uma instância de objeto ESRI, chame arcpy.AsShape()(modificado para ler polígonos com furos (consulte o ponto 2 acima)) com return cls(Array([map(lambda p: Point(*p) if p is not None else Point(), part) for part in coordinates])). A comparação é simplesmente geom1.equals(geom2)como indicado na classe Geometry .

Espero encontrar ~ 140 mudanças nas geometrias, mas meu script insiste que existem 430. Tentei verificar essas representações do GeoJSON e elas são idênticas, mas a classe Geometry equals () se recusa a dizer isso.

Um exemplo está abaixo:

>>> geom1geoJSON 
{'type': 'Polygon', 'coordinates': [[(-122.8423481559999, 47.060497293000083), (-122.84239755599992, 47.059262423000064), (-122.84416913599989, 47.059309693000046), (-122.84416913599989, 47.060497293000083), (-122.8423481559999, 47.060497293000083)]]}
>>> geom2geoJSON 
{'type': 'Polygon', 'coordinates': [[(-122.8423481559999, 47.060497293000083), (-122.84239755599992, 47.059262423000064), (-122.84416913599989, 47.059309693000046), (-122.84416913599989, 47.060497293000083), (-122.8423481559999, 47.060497293000083)]]}
>>> geom1 = arcpy.AsShape(geom1geoJSON)
>>> geom2 = arcpy.AsShape(geom2geoJSON)
>>> geom1.equals(geom2)
False
>>> geom2.equals(geom1)
False

O comportamento esperado aqui deve ser verdadeiro (não falso).

Alguém tem alguma sugestão antes de eu mudar tudo para geometrias ogr? (Eu hesito em ogr.CreateGeometryFromGeoJSON () espera uma string e o arcpy's __geo_interface__retorna um dicionário e sinto que estou adicionando complexidade extra).

Os seguintes recursos foram úteis, mesmo que eles não respondam à pergunta:

  1. pergunta arcpy.Geometry aqui em gis.stackexchange.com, que foi vinculada acima no meu texto.
  2. Erros na classe Polygon do arcpy nos fóruns do arcgis.com (aparentemente existem muitos erros de precisão no ArcGIS 10.0 que teoricamente foram corrigidos na 10.1, mas não posso verificar se, no 10.0 SP5, você ainda recebe o erro).
arcpy geometry Michalis Avraam
fonte

Respostas:

12

O problema provavelmente é o da precisão do ponto flutuante . No seu caso, você já extraiu as geometrias usando o arcpy e as combinou com o seu RUID.

Felizmente, desde que você instalou o arcpy, ficou numpy, o que facilita a comparação de conjuntos de matrizes numéricas. Nesse caso, sugiro a função numpy.allclose , disponível no numpy 1.3.0 (instalado com o ArcGIS 10).

Das amostras que você deu acima

geom1geoJSON = {'type': 'Polygon', 'coordinates': [[(-122.8423481559999, 47.060497293000083), (-122.84239755599992, 47.059262423000064), (-122.84416913599989, 47.059309693000046), (-122.84416913599989, 47.060497293000083), (-122.8423481559999, 47.060497293000083)]]}
geom2geoJSON = {'type': 'Polygon', 'coordinates': [[(-122.8423481559999, 47.060497293000083), (-122.84239755599992, 47.059262423000064), (-122.84416913599989, 47.059309693000046), (-122.84416913599989, 47.060497293000083), (-122.8423481559999, 47.060497293000083)]]}

import numpy as np

close = np.allclose(np.array(geom1geoJSON["coordinates"]), np.array(geom2geoJSON["coordinates"]), atol=1e-7)
#Returns True

A atolpalavra-chave especifica o valor da tolerância.

Observe que você não deve usar arcpy.AsShapenada. Sempre. Como observei nesta pergunta (/ shameless plug), existe um bug conhecido no ArcGIS que trunca as geometrias quando elas são criadas sem um sistema de coordenadas (mesmo depois de definir a env.XYTolerancevariável de ambiente). Em arcpy.AsShapenão há nenhuma maneira de evitar isso. Felizmente geometry.__geo_interface__, extrai as geometrias corretas das geometrias existentes (embora não lide com polígonos complexos sem a correção de @JasonScheirer).

om_henners
fonte
Obrigado. Eu não pensei em usar numpy para fazer isso. Outra solução parece estar usando o módulo decimal e trabalhar com isso, mas requer muito mais trabalho.
Michalis Avraam
Eu acho que seria importante definir o numpy.allclose() rtolparâmetro como 0. Por padrão, é 1e-05 e pode levar a uma grande tolerância se os valores das matrizes forem grandes, consulte: stackoverflow.com/a/57063678/1914034
Abaixo do radar
11

A precisão de coordenadas será uma consideração importante aqui. Números de ponto flutuante não podem ser armazenados exatamente.

Se você usar a ferramenta Comparação de recursos , ela obtém o resultado esperado usando a tolerância XY padrão?

blah238
fonte
Não verifiquei a ferramenta Comparação de recursos, pois a ferramenta que estou construindo na verdade compara recursos individuais que se movem entre diferentes classes de recursos. Ou seja, um recurso pode passar de CityRoads para CountyRoads, portanto, preciso descobrir se alguma coisa mudou na geometria e nos atributos que não sejam a classe de recurso que o contém. Há um total de 24 classes de recursos, e os recursos podem se mover entre elas. A Comparação de recursos comparará apenas duas classes de recursos, para que ele possa me dizer se não existe mais em um FC. Então eu ainda preciso comparar o recurso para garantir que ele não se alterou
Michalis Avraam
Eu verifiquei a ferramenta de comparação de recursos com a tolerância padrão (8.983e-009, que é bastante pequena, mas este é um GDB de arquivo) e ela relata algumas alterações, mas não as corretas. Especificamente, ele diz que existem 69 alterações na geometria (acho melhor do que antes), mas parece assumir que o OID é o caminho para identificar recursos exclusivos (pesquisas antigas do OID1 e do novo OID1), o que não é necessariamente verdadeiro (eu o configurei para usá-lo) meu RUID como uma espécie, mas não gostou). Então, de volta à prancheta.
Michalis Avraam
4

ao lado da resposta @ blah328, você tem a opção de comparar duas tabelas para relatar diferenças e semelhanças com valores tabulares e definições de campo com Tabela Comparada .

Exemplo:

import arcpy
from arcpy import env
arcpy.TableCompare_management(r'c:\Workspace\wells.dbf', r'c:\Workspace
\wells_new.dbf', 'WELL_ID', 'ALL', 'IGNORE_EXTENSION_PROPERTIES', 'WELL_DEPTH 0.001',
'#','CONTINUE_COMPARE', r'C:\Workspace\well_compare.txt'
Aragão
fonte
Obrigado, analisarei quando tentar comparar dados de atributos. Por enquanto, parece que não consigo comparar geometrias, o que é mais importante.
Michalis Avraam
3 
def truncateCoordinates(myGeometry)
    trucated_coords = []
    partnum = 0

    for part in (myGeometry):
        for pnt in myGeometry.getPart(partnum):
            if pnt:
                trucated_coords.append("{:10.4f}".format(pnt.X))
                trucated_coords.append("{:10.4f}".format(pnt.Y))
             else:
                continue
        partnum += 1     
    return truncated_coords

Se a .equals()função não estiver funcionando conforme o esperado e / ou as coordenadas forem ligeiramente alteradas no ArcGIS, você poderá massagear as coordenadas XY e comparar o equivalente em cadeia da geometria. Observe, truncateCoordinates()corta todos os valores além da quarta casa decimal.

geom1 = truncateCoordinates(feature1.Shape)
geom2 = truncateCoordinates(feature2.Shape)

geom1 == geom2
klewis
fonte
@ klewis- Essa é uma maneira de comparar uma geometria, mas parece que geometry.equals (geometry) deve retornar verdadeiro quando você está comparando exatamente a mesma geometria. Truncar as coordenadas é uma espécie de hack em um sentido. Talvez a ESRI precise começar a usar o tipo decimal () em vez de float se eles não puderem manipular valores de ponto flutuante corretamente internamente, mas puderem representá-los como cadeias iguais.
Michalis Avraam