Определение, окружена ли точка, используя растровую обработку

Я пытаюсь улучшить в настоящее время чрезвычайно громоздкий процесс вектор / питон для модели естественной опасности. На данный момент у нас есть длинный скрипт, который генерирует линии расстояния / пеленга из заданной точки, чтобы определить:

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

Там гораздо больше, но это суть. Я пытаюсь найти способ улучшить это, и в настоящее время нахожусь в тупике в части 3. Идея состоит в том, чтобы определить, полностью ли точка окружена полигонами, скажем, в 200 м.

Поэтому на моем прикрепленном изображении я бы хотел, чтобы точка A была отмечена как находящаяся под более высоким риском, чем точка B, поскольку она полностью окружена моими полигонами. Это повторяется примерно для 13 миллионов точек, так что это не маленькая задача, и я бы предпочел иметь поверхность для получения значений, а не запускать наш скрипт. Я думаю, что для этого должен быть вариант гидрологических инструментов или затрат, но я никак не могу обойти это.

Как я мог пойти об этом?

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

Черные клетки являются частями окружающих многоугольников. Цвета, от светло-оранжевого (короткий) до синего (длинный), показывают максимальное расстояние (максимум до 50 ячеек), которое может быть достигнуто путем прохождения линии прямой видимости без перехвата многоугольных ячеек. (Любая ячейка за пределами этого изображения рассматривается как часть многоугольников.)

Давайте обсудим эффективный способ сделать это, используя растровое представление данных. В этом представлении все «окружающие» многоугольные ячейки будут иметь, скажем, ненулевые значения, а любая ячейка, которая может быть «просвечена», будет иметь нулевое значение.

Шаг 1: Предварительный расчет структуры данных окрестности

Сначала вы должны решить, что значит для одной клетки блокировать другую. Одно из самых справедливых правил, которые я могу найти, заключается в следующем: используя интегральные координаты для строк и столбцов (и предполагая, что квадратные ячейки), давайте рассмотрим, какие ячейки могут блокировать ячейку (i, j) из представления в начале координат (0,0). Я назначаю ячейку (i ', j'), которая находится ближе всего к отрезку, соединяющему (i, j) с (0,0) среди всех ячеек, координаты которых отличаются от i и j не более чем на 1. Поскольку это не всегда дать уникальное решение (например, при (i, j) = (1,2) оба (0,1) и (1,1) будут работать одинаково хорошо), необходимы некоторые средства для устранения связей. Было бы хорошо, если бы это разрешение связей учитывало симметрию круговых окрестностей в сетках: отрицание либо координаты, либо переключение координат сохраняют эти окрестности. Поэтому мы можем решить, какие клетки блокировать (я,

Это правило иллюстрирует следующий код прототипа, написанный на R. Этот код возвращает структуру данных, которая будет удобна для определения «окруженности» произвольных ячеек в сетке.

screen <- function(k=1) {
  #
  # Returns a data structure:
  #   $offset is an array of offsets
  #   $screened is a parallel array of screened offset indexes.
  #   $distance is a parallel array of distances.
  # The first index always corresponds to (0,0).
  #
  screened.by <- function(xy) {
    uv <- abs(xy)
    if (reversed <- uv[2] > uv[1]) {
      uv <- rev(uv)
    }
    i <- which.min(c(uv[1], abs(uv[1]-uv[2]), uv[2]))
    ij <- uv + c(floor((1-i)/3), floor(i/3)-1)
    if (reversed) ij <- rev(ij)
    return(ij * sign(xy))
  }
  #
  # For each lattice point within the circular neighborhood,
  # find the unique lattice point that screens it from the origin.
  #
  xy <- subset(expand.grid(x=(-k:k), y=(-k:k)), 
               subset=(x^2+y^2 <= k^2) & (x != 0 | y != 0))
  g <- t(apply(xy, 1, function(z) c(screened.by(z), z)))
  #
  # Sort by distance from the origin.
  #
  colnames(g) <- c("x", "y", "x.to", "y.to")
  ij <- unique(rbind(g[, 1:2], g[, 3:4]))
  i <- order(abs(ij[,1]), abs(ij[,2])); ij <- ij[i, , drop=FALSE]
  rownames(ij) <- 1:length(i)
  #
  # Invert the "screened by" relation to produce the "screened" relation.
  #
  # (Row, column) offsets.
  ij.df <- data.frame(ij, i=1:length(i))
  #
  # Distances from the origin (in cells).
  distance <- apply(ij, 1, function(u) sqrt(sum(u*u)))
  #
  # "Screens" relation (represented by indexes into ij).
  g <- merge(merge(g, ij.df), ij.df, 
             by.x=c("x.to", "y.to"), by.y=c("x","y"))
  g <- subset(g, select=c(i.x, i.y))
  h <- by(g$i.y, g$i.x, identity)

  return( list(offset=ij, screened=h, distance=distance) )
}

Значение screen(12)было использовано для создания этого отношения скрининга: стрелки указывают от ячеек к тем, которые немедленно их экранируют. Оттенки пропорциональны расстоянию до начала координат, которое находится в середине этого района:

Это вычисление быстрое и должно быть сделано только один раз для данной окрестности. Например, при поиске 200 м на сетке с 5-метровыми ячейками размер окрестности будет 200/5 = 40 единиц.

Шаг 2. Применение вычислений к выбранным точкам

Все остальное просто: чтобы определить, окружена ли ячейка, расположенная в точке (x, y) (в координатах строки и столбца) относительно этой структуры данных окрестности, выполнить тест рекурсивно, начиная со смещения (i, j) = (0,0) (окрестность происхождения). Если значение в многоугольной сетке в точке (x, y) + (i, j) отлично от нуля, тогда видимость там блокируется. В противном случае нам нужно будет рассмотреть все смещения, которые могли быть заблокированы по смещению (i, j) (которое найдено за время O (1) с использованием возвращаемой структуры данных screen). Если нет ни одного заблокированного объекта, мы достигли периметра и пришли к выводу, что (x, y) не окружен, поэтому мы прекращаем вычисление (и не пытаемся осмотреть оставшиеся точки в окрестности).

Мы можем собрать еще больше полезной информации, отслеживая самое дальнее расстояние прямой видимости, достигнутое в течение алгоритма. Если это меньше, чем желаемый радиус, ячейка окружена; в противном случае это не так.

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

#
# Test a grid point `ij` for a line-of-sight connection to the perimeter
# of a circular neighborhood.  
#   `xy` is the grid.
#   `counting` determines whether to return max distance or count of stack ops.
#   `perimeter` is the assumed values beyond the extent of `xy`.
#
# Grid values of zero admit light; all others block visibility
# Returns maximum line-of-sight distance found within `nbr`.
#
panvisibility <- function(ij, xy, nbr=screen(), counting=FALSE, perimeter=1) {
  #
  # Implement a stack for the algorithm.
  #
  count <- 0 # Stack count
  stack <- list(ptr=0, s=rep(NA, dim(nbr$offset)[1]))
  push <- function(x) {
    n <- length(x)
    count <<- count+n         # For timing
    stack$s[1:n + stack$ptr] <<- x
    stack$ptr <<- stack$ptr+n
  }
  pop <- function() {
    count <<- count+1         # For timing
    if (stack$ptr <= 0) return(NULL)
    y <- stack$s[stack$ptr]
    #stack$s[stack$ptr] <<- NA # For debugging
    stack$ptr <<- stack$ptr - 1
    return(y)
  }
  #
  # Initialization.
  #
  m <- dim(xy)[1]; n <- dim(xy)[2]
  push(1) # Stack the *indexes* of nbr$offset and nbr$screened.
  dist.max <- -1
  #
  # The algorithm.
  #
  while (!is.null(i <- pop())) {
    cell <- nbr$offset[i, ] + ij
    if (cell[1] <= 0 || cell[1] > m || cell[2] <= 0 || cell[2] > n) {
      value <- perimeter
    } else {  
      value <- xy[cell[1], cell[2]]
    }
    if (value==0) {
      if (nbr$distance[i] > dist.max) dist.max <- nbr$distance[i]
      s <- nbr$screened[[paste(i)]]
      if (is.null(s)) {
        #exited = TRUE
        break
      }
      push(s)
    }
  }
  if (counting) return ( count )
  return(dist.max)
}

В этом примере многоугольные ячейки черные. Цвета дают максимальное расстояние прямой видимости (до 50 ячеек) для неполигональных ячеек: от светло-оранжевого для коротких расстояний до темно-синего для самых длинных расстояний. (Ячейки имеют ширину и высоту в одну единицу.) Видимые полосы создаются маленькими многоугольными «островками» в середине «реки»: каждая из них блокирует длинную линию других ячеек.

Анализ алгоритма

Структура стека реализует поиск в графе видимости окрестности в первую очередь для подтверждения того, что ячейка не окружена. В тех случаях, когда ячейки находятся далеко от любого многоугольника, этот поиск потребует проверки только O (k) ячеек для круговой окрестности радиуса-k. Наихудшие случаи возникают, когда в окрестности имеется небольшое количество рассеянных многоугольных ячеек, но даже при этом граница соседства не может быть полностью достигнута: для этого требуется осмотреть почти все ячейки в каждой окрестности, что является O (k ^ 2) операция.

Следующее поведение типично для того, что будет встречаться. Для малых значений k, если многоугольники не заполняют большую часть сетки, большинство неполигональных ячеек будут явно необоснованными, и алгоритм масштабируется как O (k). Для промежуточных значений масштабирование начинает выглядеть как O (k ^ 2). Поскольку k становится действительно большим, большинство ячеек будут окружены, и этот факт можно будет определить задолго до проверки всей окрестности: тем самым вычислительные усилия алгоритма достигают практического предела. Этот предел достигается, когда радиус окрестности приближается к диаметру самых больших соединенных неполигональных областей в сетке.

В качестве примера я использовал countingопцию, закодированную в прототипе screenдля возврата количества операций стека, используемых в каждом вызове. Это измеряет вычислительные усилия. На следующем графике показано среднее число операций стека в зависимости от радиуса окрестности. Это демонстрирует предсказанное поведение.

Мы можем использовать это для оценки вычислений, необходимых для оценки 13 миллионов точек на сетке. Предположим, что используется окрестность k = 200/5 = 40. Тогда в среднем потребуется несколько сотен операций стека (в зависимости от сложности сетки многоугольников и от того, где расположены 13 миллионов точек относительно многоугольников), что означает, что в эффективном скомпилированном языке не более нескольких тысяч простых числовых операций потребуется (сложение, умножение, чтение, запись, смещение и т. д.). Большинство компьютеров смогут оценить окружение примерно в миллион очков с такой скоростью. (TheRреализация намного, намного медленнее, чем это, потому что это плохо в алгоритме такого рода, поэтому его можно рассматривать только как прототип.) Соответственно, мы могли бы надеяться, что эффективная реализация на достаточно эффективном и соответствующем языке - C ++ и Python приходит на ум - может завершить оценку 13 миллионов точек в минуту или меньше, при условии, что вся многоугольная сетка находится в оперативной памяти.

Когда сетка слишком велика для размещения в ОЗУ, эту процедуру можно применить к мозаичным частям сетки. Они должны перекрываться только kстроками и столбцами; возьмите максимумы в перекрытиях при составлении мозаики результатов.

Другие приложения

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

Я определенно вижу, как можно сделать это с помощью растрового решения, но, учитывая даже уменьшенное количество точек, я ожидал бы очень большое / высокое разрешение и, следовательно, сложную обработку сетки или набора сеток. Учитывая это, мне интересно, может ли использование топологии в GDB быть более эффективным. Вы можете найти все внутренние пустоты с чем-то вроде:

arcpy.env.workspace = 'myGDB'
arcpy.CreateTopology_management('myGDB', 'myTopology', '')    
arcpy.AddFeatureClassToTopology_management('myTopology', 'myFeatures', '1','1')    
arcpy.AddRuleToTopology_management ('myToplogy', 'Must Not Have Gaps (Area)', 'myFeatures', '', '', '')    
arcpy.ValidateTopology_management('myTopology', 'Full_Extent')
arcpy.ExportTopologyErrors_management('myTopology', 'myGDB', 'topoErrors')
arcpy.FeatureToPolygon_management('topoErrors_line','topoErrorsVoidPolys', '0.1')`

затем вы можете работать topoErrorsVoidPolysв обычном порядке Intersect_analysis()или как угодно. Возможно, вам придется возиться с извлечением интересующих вас политик topoErrorsVoidPolys. У @whuber есть много отличных постов о подобных вещах в других местах здесь, на gis.stackexchange.com.

Если вы действительно хотите сделать растр ... я бы сделал что-то в духе этого псевдокода (не смущайтесь, потому что очевидно, что я возвращаю AML!: P)

1. растеризация точек ("pts_g") и polys ("polys_g" (
2. пустоты = региональная группа (con (isnull (polys_g), 1))
3. возможно, потребуется что-то сделать, чтобы очистить пустоты, чтобы устранить нежелательные внешние полигоны / открытую область вселенной
4. pts_surounded = con (пустоты, pts_g)

Просто придумал, так что, возможно, потребуется доработка.

Если вы используете пороговое значение расстояния (здесь вы говорите о 200 м), то лучшим решением будет использовать векторный анализ:

1) создайте 200-метровый буфер вокруг каждой точки (черным на иллюстрации)

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

3) использовать функцию полигона (управление) для создания полигонов, где ваши точки полностью окружены (красным на иллюстрации)

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

Варианты 2b) В зависимости от ваших потребностей, вы можете выбрать расположение окружающих полигонов на расстоянии 200 м, затем вы можете определить некоторые виды «вложения», но не так строго, как в 2).

Рассматривая «случай лабиринта», это могло бы помочь: оценить, сколько времени нужно «убежать» из локации.

• Вы уже можете исключить из анализа баллы, которые полностью включены или полностью свободны

• затем вы преобразуете свои препятствия в растр и устанавливаете значения для NoData, где у вас есть многоугольник, и для размера ячейки в метрах, где вы этого не делаете (это сделает ваш растр стоимости).

• в-третьих, вы можете вычислить расстояние затрат, используя только что созданный растр затрат

• наконец, вы используете зональную статистику в качестве таблицы, основанной на границах буфера, преобразованных в растр (образуя кольцо). Если вы можете убежать во всех направлениях, минимум должен быть примерно 200 (в зависимости от размера ячейки вашего анализа). Но если вы находитесь в лабиринте, максимум будет больше 200. Таким образом, максимум зональной статистики минус 200 будет непрерывным значением, указывающим, насколько трудно «сбежать».