Перекрытия всех объектов в одном наборе данных в arcpy (st_overlaps)

Question

Пакет R sf имеет удивительный набор функций под названием «Геометрические двоичные предикаты», которые подробно описаны здесь .

Как указано в ссылке,функции рекурсивно применяются ко всем геометриям в одном и том же наборе данных , если предоставляется только один sf -объект (см. пример ниже)

Если y отсутствует, st_predicate(x, x)фактически вызывается, и квадратная матрица возвращается с диагональными элементами st_predicate(x[i], x[i]).

Однако сейчас я создаю некоторые инструменты, где я привязан к arcpy из ArcGIS.Какой быстрый способ получить квадратную матрицу всех объектов в одном наборе данных, указывающую, перекрываются ли соответствующие объекты или нет?

arcpy.SpatialJoin_analysis() сравнивает только два набора данных и arcpy.GenerateNearTable_analysis(), а также arcpy.Near_analysis()вычислять только расстояния между объектами.

Вот как st_overlaps() работает в R:

library(sf)
#> Warning: Paket 'sf' wurde unter R Version 3.5.2 erstellt
#> Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3

b0 = st_polygon(list(rbind(c(-1,-1), c(1,-1), c(1,1), c(-1,1), c(-1,-1))))
a0 = b0 * 0.8
a1 = a0 * 0.5 + c(2, 0.7)
a2 = a0 + 1
a3 = b0 * 0.5 + c(2, -0.5)
y = st_sfc(a0,a1,a2,a3)

plot(y)

st_overlaps(y,sparse = F)
#>       [,1]  [,2]  [,3]  [,4]
#> [1,] FALSE FALSE  TRUE FALSE
#> [2,] FALSE FALSE  TRUE FALSE
#> [3,]  TRUE  TRUE FALSE FALSE
#> [4,] FALSE FALSE FALSE FALSE

^{Создано в2019-04-16 по представлению пакета (v0.2.1)}

Answer 1

Один из способов сделать это:

1. с помощью инструмента «Пересечь», чтобы разбить все полигоны на части, где они перекрываются.
2. С помощью инструмента «Найти дубликаты» сгенерироватьсписок полигонов, имеющих одинаковую форму.
3. Присоедините результаты этих процессов к исходной таблице, используя идентификаторы объектов.Как только объединения будут созданы, вы увидите перекрывающиеся полигоны, отмеченные одинаковым значением в поле «Feat_Seq».

Пример Python:

arcpy.analysis.Intersect("test.shp", "test_Intersect", "ONLY_FID", None, "INPUT")
arcpy.management.FindIdentical("test_Intersect", r"test_Intersect_FindIdentical", "Shape", None, 0, "ONLY_DUPLICATES")
arcpy.management.AddJoin("test", "FID", "test_Intersect", "FID_test", "KEEP_ALL")
arcpy.management.AddJoin("test", "test_Intersect.OBJECTID", "test_Intersect_FindIdentical", "IN_FID", "KEEP_ALL")