Как создать ориентированный буфер, используя arcpy?

Я хотел бы создать ориентированный буфер для каждого многоугольника в моем шейп-файле, используя arcpy. Под ориентацией я подразумеваю, что у меня есть два угла a1 и a2, которые ограничивают направление буфера. Это представлено на графике ниже:

Любые идеи?

Резюме

Этот ответ помещает вопрос в более широкий контекст, описывает эффективный алгоритм, применимый к представлению объектов в шейп-файле (в виде «векторов» или «строк линий» точек), показывает некоторые примеры его применения и дает рабочий код для использования или переноса в среда ГИС.

Фон

Это пример морфологического расширения. В полной общности дилатация «распространяет» точки региона в их окрестности; коллекция точек, где они заканчиваются, является «расширением». Приложения в ГИС многочисленны: моделирование распространения огня, движение цивилизаций, распространение растений и многое другое.

Математически и в очень большой (но полезной) общности расширение расширяет набор точек риманова многообразия (таких как плоскость, сфера или эллипсоид). Распространение обусловлено подмножеством касательного расслоения в этих точках. Это означает, что в каждой из точек задан набор векторов (направлений и расстояний) (я называю это «окрестностью»); каждый из этих векторов описывает геодезический путь, начинающийся в его базовой точке. Базовая точка «распространяется» на концы всех этих путей. (Гораздо более ограниченное определение «растяжения», которое обычно используется при обработке изображений, см . В статье в Википедии . Функция распространения называется экспоненциальной картой. в дифференциальной геометрии.)

«Буферизация» объекта является одним из простейших примеров такого расширения: диск с постоянным радиусом (буферным радиусом) создается (по крайней мере, концептуально) вокруг каждой точки объекта. Объединение этих дисков является буфером.

Этот вопрос требует вычисления немного более сложного расширения, где распространение может происходить только в пределах заданного диапазона углов (то есть в круглом секторе). Это имеет смысл только для элементов, которые не имеют заметно изогнутой поверхности (например, небольшие элементы на сфере или эллипсоиде или любые элементы на плоскости). Когда мы работаем в самолете, также имеет смысл ориентировать все сектора в одном направлении. (Однако, если бы мы моделировали распространение огня по ветру, мы бы хотели, чтобы сектора были ориентированы по ветру, и их размеры могли бы также изменяться в зависимости от скорости ветра: это одна из причин для общего определения расширения, которое я дал. ) (На изогнутых поверхностях, таких как эллипсоид, вообще невозможно ориентировать все сектора в одном и том же направлении.)

При следующих обстоятельствах дилатацию относительно легко вычислить:

• Объект находится на плоскости (то есть мы расширяем карту объекта и, надеюсь, карта достаточно точная).

• Расширение будет постоянным : распространение в каждой точке объекта будет происходить в конгруэнтных окрестностях с одинаковой ориентацией.

• Эта общая окрестность является выпуклой. Выпуклость значительно упрощает и ускоряет вычисления.

Этот вопрос подходит для таких специализированных обстоятельств: он требует расширения произвольных многоугольников круговыми секторами, чьи источники (центры дисков, из которых они вышли) расположены в базовых точках. Если эти сектора не охватывают более 180 градусов, они будут выпуклыми. (Большие сектора всегда можно разделить пополам на два выпуклых сектора; объединение двух меньших расширений даст желаемый результат.)

Реализация

Поскольку мы проводим Евклида вычисления - делают распространение в плоскости - мы можем растянуть точку только путем перевода дилатации окрестности этой точки. (Чтобы быть в состоянии сделать это, район нуждается в происхожденииэто будет соответствовать базовой точке. Например, источником секторов в этом вопросе является центр круга, из которого они образованы. Это происхождение лежит на границе сектора. В стандартной операции буферизации ГИС окрестность представляет собой круг с началом в центре; Теперь происхождение лежит внутри круга. Выбор источника не имеет большого значения в вычислительном отношении, потому что изменение происхождения просто сдвигает всю дилатацию, но может иметь большое значение с точки зрения моделирования природных явлений. sectorФункции в коде ниже показано , как может быть определено происхождение.)

Расширение отрезка может быть сложным, но для выпуклой окрестности мы можем создать расширение как объединение расширений двух конечных точек вместе с тщательно выбранным параллелограммом. (В интересах пространства я не буду останавливаться, чтобы доказать подобные математические утверждения, но призываю читателей попробовать свои собственные доказательства, потому что это проницательное упражнение.) Вот иллюстрация, использующая три сектора (показаны розовым цветом). Они имеют единичные радиусы и их углы указаны в заголовках. Сам отрезок имеет длину 2, горизонтальный и показан черным цветом:

Параллелограммы находятся путем нахождения розовых точек, которые находятся как можно дальше от сегмента только в вертикальном направлении . Это дает две нижние точки и две верхние точки вдоль линий, параллельных сегменту. Нам просто нужно объединить четыре точки в параллелограмм (показан синим цветом). Обратите внимание, справа, как это имеет смысл, даже если сам сектор является просто отрезком (а не истинным многоугольником): там каждая точка на отрезке была переведена в направлении 171 градуса к востоку от севера на расстояние от 0 до 1. Множество этих конечных точек является показанным параллелограммом. Детали этого вычисления появляются в bufferфункции, определенной dilate.edgesв приведенном ниже коде.

Чтобы расширить полилинию , мы формируем объединение растяжений точек и отрезков, которые ее образуют. Последние две строчки dilate.edgesвыполняют этот цикл.

Расширение многоугольника требует включения его внутренней части вместе с расширением его границы. (Это утверждение делает некоторые предположения о окрестности дилатации. Во-первых, все окрестности содержат точку (0,0), которая гарантирует, что многоугольник включен в его дилатацию. В случае переменных окрестностей также предполагается, что дилатация любой внутренней части точка многоугольника не будет выходить за границы растяжения за пределы границ. Это относится к постоянным окрестностям.)

Давайте рассмотрим некоторые примеры того, как это работает, сначала с неагоном (выбранным для раскрытия деталей), а затем с кружком (выбранным для соответствия иллюстрации в вопросе). В примерах будут по-прежнему использоваться те же три окрестности, но суженные до радиуса 1/3.

На этом рисунке внутренняя часть многоугольника серая, точечные дилатации (сектора) розовые, а краевые дилатации (параллелограммы) синие.

«Круг» на самом деле всего лишь 60 гонов, но он хорошо приближается к кругу.

Представление

Когда базовый элемент представлен N точками, а окрестность дилатации - M точками, этот алгоритм требует усилий O (N M) . Это должно сопровождаться упрощением путаницы вершин и ребер в объединении, что может потребовать усилий O (N M log (N M)): это то, что нужно просить ГИС; мы не должны программировать это.

Вычислительное усилие может быть улучшено до O (M + N) для выпуклых базовых объектов (потому что вы можете решить, как перемещаться вокруг новой границы, соответствующим образом объединяя списки вершин, описывающих границы исходных двух фигур). Это не потребует какой-либо последующей очистки.

Когда окрестность дилатации медленно меняет размер и / или ориентацию по мере продвижения вокруг базового объекта, расширение края может быть близко аппроксимировано от выпуклой оболочки объединения расширений его конечных точек. Если две окрестности дилатации имеют точки M1 и M2, это можно найти с помощью усилия O (M1 + M2), используя алгоритм, описанный в Shamos & Preparata, Computational Geometry . Следовательно, если K = M1 + M2 + ... + M (N) быть общим числом вершин в N окрестностях дилатации, мы можем вычислить дилатацию за O (K * log (K)).

Зачем нам хотеть заняться таким обобщением, если все, что мы хотим, это простой буфер? Для больших объектов на Земле окрестность дилатации (например, диск), которая в действительности имеет постоянный размер, может иметь различный размер на карте, где выполняются эти вычисления. Таким образом, мы получаем способ производить точные вычисления для эллипсоида , продолжая пользоваться всеми преимуществами евклидовой геометрии.

Код

Примеры были созданы с этим Rпрототипом, который можно легко перенести на ваш любимый язык (Python, C ++ и т. Д.). По своей структуре он соответствует анализу, представленному в этом ответе, и поэтому не требует отдельного объяснения. Комментарии уточняют некоторые детали.

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

#
# Dilate the vertices of a polygon/polyline by a shape.
#
dilate.points <- function(p, q) {
  # Translate a copy of `q` to each vertex of `p`, resulting in a list of polygons.
  pieces <- apply(p, 1, function(x) list(t(t(q)+x)))
  lapply(pieces, function(z) z[[1]]) # Convert to a list of matrices
}
#
# Dilate the edges of a polygon/polyline `p` by a shape `q`. 
# `p` must have at least two rows.
#
dilate.edges <- function(p, q) {
  i <- matrix(c(0,-1,1,0), 2, 2)       # 90 degree rotation
  e <- apply(rbind(p, p[1,]), 2, diff) # Direction vectors of the edges
  # Dilate a single edge from `x` to `x+v` into a parallelogram
  # bounded by parts of the dilation shape that are at extreme distances
  # from the edge.
  buffer <- function(x, v) {
    y <- q %*% i %*% v # Signed distances orthogonal to the edge
    k <- which.min(y)  # Find smallest distance, then the largest *after* it
    l <- (which.max(c(y[-(1:k)], y[1:k])) + k-1) %% length(y)[1] + 1
    list(rbind(x+q[k,], x+v+q[k,], x+v+q[l,], x+q[l,])) # A parallelogram
  }
  # Apply `buffer` to every edge.
  quads <- apply(cbind(p, e), 1, function(x) buffer(x[1:2], x[3:4]))
  lapply(quads, function(z) z[[1]]) # Convert to a list of matrices
}
#----------------------- (This ends the dilation code.) --------------------------#
#
# Display a polygon and its point and edge dilations.
# NB: In practice we would submit the polygon, its point dilations, and edge 
#     dilations to the GIS to create and simplify their union, producing a single
#     polygon.  We keep the three parts separate here in order to illustrate how
#     that polygon is constructed.
#
display <- function(p, d.points, d.edges, ...) {
  # Create a plotting region covering the extent of the dilated figure.
  x <- c(p[,1], unlist(lapply(c(d.points, d.edges), function(x) x[,1])))
  y <- c(p[,2], unlist(lapply(c(d.points, d.edges), function(x) x[,2])))
  plot(c(min(x),max(x)), c(min(y),max(y)), type="n", asp=1, xlab="x", ylab="y", ...)
  # The polygon itself.
  polygon(p, density=-1, col="#00000040")
  # The dilated points and edges.
  plot.list <- function(l, c) lapply(l, function(p) 
                  polygon(p, density=-1, col=c, border="#00000040"))
  plot.list(d.points, "#ff000020")
  plot.list(d.edges, "#0000ff20")
  invisible(NULL) # Doesn't return anything
}
#
# Create a sector of a circle.
# `n` is the number of vertices to use for approximating its outer arc.
#
sector <- function(radius, arg1, arg2, n=1, origin=c(0,0)) {
  t(cbind(origin, radius*sapply(seq(arg1, arg2, length.out=n), 
                  function(a) c(cos(a), sin(a)))))
}
#
# Create a polygon represented as an array of rows.
#
n.vertices <- 60 # Inscribes an `n.vertices`-gon in the unit circle.
angles <- seq(2*pi, 0, length.out=n.vertices+1)
angles <- angles[-(n.vertices+1)]
polygon.the <- cbind(cos(angles), sin(angles))
if (n.vertices==1) polygon.the <- rbind(polygon.the, polygon.the)
#
# Dilate the polygon in various ways to illustrate.
#
system.time({
  radius <- 1/3
  par(mfrow=c(1,3))
  q <- sector(radius, pi/12, 2*pi/3, n=120)
  d.points <- dilate.points(polygon.the, q)
  d.edges <- dilate.edges(polygon.the, q)
  display(polygon.the, d.points, d.edges, main="-30 to 75 degrees")

  q <- sector(radius, pi/3, 4*pi/3, n=180)
  d.points <- dilate.points(polygon.the, q)
  d.edges <- dilate.edges(polygon.the, q)
  display(polygon.the, d.points, d.edges, main="-150 to 30 degrees")

  q <- sector(radius, -9/20*pi, -9/20*pi)
  d.points <- dilate.points(polygon.the, q)
  d.edges <- dilate.edges(polygon.the, q)
  display(polygon.the, d.points, d.edges, main="171 degrees")
})

Время вычисления для этого примера (из последнего рисунка) с N = 60 и M = 121 (слева), M = 181 (в середине) и M = 2 (справа) составило одну четверть секунды. Тем не менее, большая часть этого была для дисплея. Как правило, этот Rкод будет обрабатывать около N M = 1,5 миллиона в секунду (что занимает всего 0,002 секунды или около того, чтобы выполнить все показанные в примере вычисления). Тем не менее, появление продукта M N подразумевает дилатацию множества фигур или сложных фигур через детальное соседство, может занять значительное время, так что будьте осторожны! Оцените сроки для небольших проблем, прежде чем заняться большой. В таких обстоятельствах можно обратиться к растровому решению (которое гораздо проще реализовать, требуя по существу только одного вычисления соседства).

Это довольно широко, но вы могли бы:

1. буфер исходного полигона
2. найти точку начала «ориентированных» лучей, которые будут созданы на границе многоугольника (какая-то точка касания?)
3. создайте / вытяните линию из этой точки на расстояние за пределами буфера, используя угол, обсуждаемый в комментариях к вопросу.
4. пересекают эту линию с буфером и исходным многоугольником. Вероятно, это можно сделать одновременно с 3) с соответствующими аргументами Extend.
5. извлечь новый полигон "ориентированный буфер" из результирующего набора полигонов