Сравнивая две геометрии в ArcPy?

18

Я пытаюсь сравнить два отдельных класса пространственных объектов, чтобы определить различия между ними (своего рода функция сравнения). Мой основной рабочий процесс:

  1. Я извлекаю геометрию с помощью SearchCursor
  2. Сохраните геометрию двух классов пространственных объектов как GeoJSON, используя модифицированный __geo_interface__(полученный из valveLondon return {'type': 'Polygon', 'coordinates': [[((pt.X, pt.Y) if pt else None) for pt in part] for part in self]} ). Это делается для того, чтобы избежать использования объекта общей геометрии, который ESRI использует с курсорами, и невозможности делать глубокие копии (об этом говорят некоторые обсуждения здесь на gis.stackexchange).
  3. Проверьте геометрию двух классов объектов на основе уникального идентификатора. Например, сравните геометрию OID1 FC1 с геометрией OID1 FC2. Чтобы получить геометрию как экземпляр объекта ESRI, вызовите arcpy.AsShape()(модифицированный для чтения полигонов с отверстиями (см. Пункт 2 выше) с помощью return cls(Array([map(lambda p: Point(*p) if p is not None else Point(), part) for part in coordinates])). Сравнение просто, geom1.equals(geom2)как указано в классе геометрии .

Я ожидаю найти ~ 140 изменений в геометрии, но мой сценарий настаивает на том, что их 430. Я пытался проверить эти представления GeoJSON, и они идентичны, но класс геометрии equals () отказывается это говорить.

Пример ниже:

>>> 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

Ожидаемое поведение здесь должно быть True (не False).

У кого-нибудь есть какие-либо предложения, прежде чем я перенесу все в ogr геометрии? (Я сомневаюсь, так как ogr.CreateGeometryFromGeoJSON () ожидает строку, а arcpy __geo_interface__возвращает словарь, и я чувствую, что добавляю дополнительную сложность).

Нашли следующие полезные ресурсы, хотя они не отвечают на вопрос:

  1. arcpy.Geometry вопрос здесь на gis.stackexchange.com, который был связан выше в моем тексте.
  2. Ошибки в классе Polygon в Arcpy с форумов arcgis.com (очевидно, существует множество ошибок точности в ArcGIS 10.0, которые теоретически были исправлены в 10.1, но я не могу проверить это, в 10.0 SP5 вы все еще получаете ошибку).
Михалис Авраам
источник

Ответы:

12

Проблема, скорее всего, связана с точностью с плавающей точкой . В вашем случае вы уже извлекли геометрии, используя arcpy, и сопоставили их с вашим RUID.

К счастью, так как у вас установлена ​​arcpy, у вас есть numpy, что упрощает сравнение наборов числовых массивов. В этом случае я бы предложил функцию numpy.allclose , которая доступна в numpy 1.3.0 (установлена ​​с ArcGIS 10).

Из образцов, которые вы дали выше

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

atolКлючевое слово указывает значение допуска.

Обратите внимание, что вы не должны использовать arcpy.AsShapeвообще. Когда-либо. Как я уже отмечал в этом вопросе (/ shameless plug), в ArcGIS есть известная ошибка, которая усекает геометрии, когда они создаются без системы координат (даже после установки env.XYToleranceпеременной окружения). Там arcpy.AsShapeнет никакого способа избежать этого. К счастью, geometry.__geo_interface__он извлекает правильные геометрии из существующих геометрий (хотя он не обрабатывает сложные полигоны без исправления из @JasonScheirer).

om_henners
источник
Спасибо. Я не думал использовать numpy для этого. Другое решение, похоже, использует десятичный модуль и работает через него, но требует гораздо больше работы.
Михалис Авраам
Я думаю, что было бы важно установить numpy.allclose() rtolпараметр на 0. По умолчанию это 1e-05, и это может привести к большому допуску, если значения массивов большие, см. Stackoverflow.com/a/57063678/1914034
Ниже радара
11

Точность координат будет важным фактором здесь. Числа с плавающей точкой не могут быть сохранены точно.

Если вы используете инструмент сравнения функций , он дает ожидаемый результат, используя допуск XY по умолчанию?

blah238
источник
Я не проверял инструмент сравнения объектов, так как инструмент, который я создаю, фактически сравнивает отдельные объекты, которые перемещаются между различными классами объектов. То есть объект может перемещаться из CityRoads в CountyRoads, поэтому мне нужно выяснить, изменилось ли что-либо в геометрии и атрибутах, отличных от класса объектов, в котором она содержится. Всего существует 24 класса пространственных объектов, и объекты могут перемещаться между ними. Сравнение функций будет сравнивать только 2 класса объектов, поэтому он может сказать мне, если он больше не существует в FC. Затем мне все еще нужно сравнить функцию, чтобы убедиться, что она не изменилась
Михалис Авраам
Я проверил инструмент сравнения функций с допуском по умолчанию (8.983e-009, который довольно мал, но это файловая GDB), и он сообщает о некоторых изменениях, но не о правильных. В частности, говорится, что есть 69 изменений геометрии (я думаю, что лучше, чем раньше), но, похоже, предполагается, что OID - это способ идентифицировать уникальные объекты (ищет старый OID1 и новый OID1), что не обязательно верно (я настроил его на использование мой РУИД как род, но это не понравилось). Итак, вернемся к чертежной доске.
Михалис Авраам
4

Помимо ответа @ blah328, у вас есть выбор - сравнить две таблицы для отчета о различиях и сходствах с табличными значениями и определениями полей с помощью функции « Сравнение таблиц» .

Пример:

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' 
Арагон
источник
Спасибо, я посмотрю на это, когда попытаюсь сравнить данные атрибутов. На данный момент, похоже, я не могу сравнивать геометрии, что более важно.
Михалис Авраам
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

Если .equals()функция не работает должным образом и / или координаты немного изменены в ArcGIS, вы можете помассировать координаты XY, а затем сравнить строковый эквивалент геометрии. Обратите внимание, truncateCoordinates()отсекает все значения за пределами 4-го знака после запятой.

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

geom1 == geom2
klewis
источник
@ klewis- Это один из способов сравнения геометрии, но похоже, что geometry.equals (geometry) должен возвращать true при сравнении одной и той же геометрии. Усечение координат является в некотором смысле взломом. Возможно, ESRI нужно начать использовать тип decimal () вместо float, если они не могут правильно обрабатывать значения с плавающей запятой, но могут представлять их как равные строки.
Михалис Авраам
1

Вы можете использовать инструмент Выбрать слой по расположению (Управление данными) с параметром overlap_type «ARE_IDENTICAL_TO», переключить выбор , проверить количество строк и затем выполнить цикл по строкам для сбора объектов или любой другой соответствующей информации.

Брент Эдвардс
источник