Вычислить максимальное расстояние внутри полигона в направлении х (восток-запад) в PostGIS?

Меня интересует максимальная ширина многоугольника, например озера, в направлении восток-запад. Ограничительные рамки помогут только в простых многоугольниках, но не в сложных вогнутых многоугольниках.

Это, вероятно, требует некоторых сценариев на любой платформе ГИС.

Наиболее эффективный метод (асимптотически) - это развертка по вертикальной линии: он требует сортировки ребер по их минимальным y-координатам и последующей обработки ребер снизу (минимум y) до верха (максимум y) для O (e * log ( e)) алгоритм, когда задействованы e ребра.

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

Мы специально будем искать горизонтальный сегмент (ы) максимальной длины, полностью лежащий в замыкании многоугольника. Например, это устраняет колебание между x = 20 и x = 40, исходящее из отверстия между x = 10 и x = 25. Тогда легко показать, что, по крайней мере, один из горизонтальных сегментов максимальной длины пересекает, по крайней мере, одну вершину. (Если есть решения не пересекающиеся вершин , они будут лежать внутри некоторого параллелограмма , ограниченного сверху и снизу решениями , которые делают пересекаются по меньшей мере , одну вершины. Это дает нам средство , чтобы найти все решения.)

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

Например, рассмотрим состояние при достижении уровня y = 10. Слева направо мы находим следующие ребра:

      x.min x.max y.min y.max
 [1,]    10     0     0    30
 [2,]    10    24    10    20
 [3,]    20    24    10    20
 [4,]    20    40    10    10
 [5,]    40    20    10    10
 [6,]    60     0     5    30
 [7,]    60    60     5    30
 [8,]    60    70     5    20
 [9,]    60    70     5    15
[10,]    90   100    10    40

В этой таблице (x.min, y.min) - это координаты нижней конечной точки ребра, а (x.max, y.max) - координаты ее верхней конечной точки. На этом уровне (y = 10) первый ребро пересекается внутри, второй пересекается снизу и так далее. Некоторые ребра, оканчивающиеся на этом уровне, например от (10,0) до (10,10), не включены в список.

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

Продолжая пример, вот отсортированный список x-координат, где негоризонтальные ребра начинаются с линии y = 10 или пересекают ее:

x.array    6.7 10 20 48 60 63.3 65 90
interior     1  0  1  0  1    0  1  0

(Обратите внимание, что x = 40 отсутствует в этом списке.) Значения interiorмассива отмечают левые конечные точки внутренних сегментов: 1 обозначает внутренний интервал, 0 - внешний интервал. Таким образом, первая 1 указывает, что интервал от х = 6,7 до х = 10 находится внутри многоугольника. Следующий 0 указывает, что интервал от х = 10 до х = 20 находится за пределами многоугольника. И так продолжается: массив идентифицирует четыре отдельных интервала как внутри многоугольника.

Некоторые из этих интервалов, например интервал от x = 60 до x = 63,3, не пересекают никакие вершины: быстрая проверка по x-координатам всех вершин с y = 10 устраняет такие интервалы.

Во время сканирования мы можем отслеживать длину этих интервалов, сохраняя данные, относящиеся к максимальному интервалу (длинам), найденному до сих пор.

Обратите внимание на некоторые последствия этого подхода. Вершина в форме буквы «V», если встречается, является источником двух ребер. Поэтому два перехода происходят при его пересечении. Эти переключатели отменяются. Любая перевернутая буква "v" даже не обрабатывается, поскольку оба ее края удаляются перед началом сканирования слева направо. В обоих случаях такая вершина не блокирует горизонтальный сегмент.

Более двух ребер могут иметь общую вершину: это показано в (10,0), (60,5), (25, 20) и - хотя это трудно сказать - в (20,10) и (40). , 10). (Это потому, что болтается (20,10) -> (40,10) -> (40,0) -> (40, -50) -> (40, 10) -> (20, 10). Обратите внимание, что вершина в (40,0) также находится внутри другого ребра ... это противно.) Этот алгоритм прекрасно справляется с этими ситуациями.

Сложная ситуация проиллюстрирована в самом низу: x-координаты не горизонтальных сегментов есть

30, 50

Это приводит к тому, что все слева от x = 30 считается внешним, все между 30 и 50 - внутренним, а все после 50 снова внешним. Вершина в x = 40 никогда не рассматривается даже в этом алгоритме.

Вот как выглядит полигон в конце сканирования. Я показываю все внутренние вершины, содержащие вершины, темно-серым цветом, любые интервалы максимальной длины - красным, и окрашиваю вершины в соответствии с их y-координатами. Максимальный интервал составляет 64 единицы.

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

Если координаты географические , то горизонтальные линии действительно находятся на кругах широты. Их длины не сложно вычислить: просто умножьте их евклидову длину на косинус их широты (в сферической модели). Поэтому этот алгоритм хорошо адаптируется к географическим координатам. (Чтобы справиться с намоткой на скважину с меридианом + -180, сначала может потребоваться найти кривую от южного полюса до северного полюса, которая не проходит через многоугольник. После повторного выражения всех x-координат в виде горизонтальных смещений относительно этого кривая, этот алгоритм будет правильно найти максимальный горизонтальный сегмент.)

Ниже приведен Rкод, реализованный для выполнения расчетов и создания иллюстраций.

#
# Plotting functions.
#
points.polygon <- function(p, ...) {
  points(p$v, ...)
}
plot.polygon <- function(p, ...) {
  apply(p$e, 1, function(e) lines(matrix(e[c("x.min", "x.max", "y.min", "y.max")], ncol=2), ...))
}
expand <- function(bb, e=1) {
  a <- matrix(c(e, 0, 0, e), ncol=2)
  origin <- apply(bb, 2, mean)
  delta <-  origin %*% a - origin
  t(apply(bb %*% a, 1, function(x) x - delta))
}
#
# Convert polygon to a better data structure.
#
# A polygon class has three attributes:
#   v is an array of vertex coordinates "x" and "y" sorted by increasing y;
#   e is an array of edges from (x.min, y.min) to (x.max, y.max) with y.max >= y.min, sorted by y.min;
#   bb is its rectangular extent (x0,y0), (x1,y1).
#
as.polygon <- function(p) {
  #
  # p is a list of linestrings, each represented as a sequence of 2-vectors 
  # with coordinates in columns "x" and "y". 
  #
  f <- function(p) {
    g <- function(i) {
      v <- p[(i-1):i, ]
      v[order(v[, "y"]), ]
    }
    sapply(2:nrow(p), g)
  }
  vertices <- do.call(rbind, p)
  edges <- t(do.call(cbind, lapply(p, f)))
  colnames(edges) <- c("x.min", "x.max", "y.min", "y.max")
  #
  # Sort by y.min.
  #
  vertices <- vertices[order(vertices[, "y"]), ]
  vertices <- vertices[!duplicated(vertices), ]
  edges <- edges[order(edges[, "y.min"]), ]

  # Maintaining an extent is useful.
  bb <- apply(vertices <- vertices[, c("x","y")], 2, function(z) c(min(z), max(z)))

  # Package the output.
  l <- list(v=vertices, e=edges, bb=bb); class(l) <- "polygon"
  l
}
#
# Compute the maximal horizontal interior segments of a polygon.
#
fetch.x <- function(p) {
  #
  # Update moves the line from the previous level to a new, higher level, changing the
  # state to represent all edges originating or strictly passing through level `y`.
  #
  update <- function(y) {
    if (y > state$level) {
      state$level <<- y
      #
      # Remove edges below the new level from state$current.
      #
      current <- state$current
      current <- current[current[, "y.max"] > y, ]
      #
      # Adjoin edges at this level.
      #
      i <- state$i
      while (i <= nrow(p$e) && p$e[i, "y.min"] <= y) {
        current <- rbind(current, p$e[i, ])
        i <- i+1
      }
      state$i <<- i
      #
      # Sort the current edges by x-coordinate.
      #
      x.coord <- function(e, y) {
        if (e["y.max"] > e["y.min"]) {
          ((y - e["y.min"]) * e["x.max"] + (e["y.max"] - y) * e["x.min"]) / (e["y.max"] - e["y.min"])
        } else {
          min(e["x.min"], e["x.max"])
        }
      }
      if (length(current) > 0) {
        x.array <- apply(current, 1, function(e) x.coord(e, y))
        i.x <- order(x.array)
        current <- current[i.x, ]
        x.array <- x.array[i.x]     
        #
        # Scan and mark each interval as interior or exterior.
        #
        status <- FALSE
        interior <- numeric(length(x.array))
        for (i in 1:length(x.array)) {
          if (current[i, "y.max"] == y) {
            interior[i] <- TRUE
          } else {
            status <- !status
            interior[i] <- status
          }
        }
        #
        # Simplify the data structure by retaining the last value of `interior`
        # within each group of common values of `x.array`.
        #
        interior <- sapply(split(interior, x.array), function(i) rev(i)[1])
        x.array <- sapply(split(x.array, x.array), function(i) i[1])

        print(y)
        print(current)
        print(rbind(x.array, interior))


        markers <- c(1, diff(interior))
        intervals <- x.array[markers != 0]
        #
        # Break into a list structure.
        #
        if (length(intervals) > 1) {
          if (length(intervals) %% 2 == 1) 
            intervals <- intervals[-length(intervals)]
          blocks <- 1:length(intervals) - 1
          blocks <- blocks - (blocks %% 2)
          intervals <- split(intervals, blocks)  
        } else {
          intervals <- list()
        }
      } else {
        intervals <- list()
      }
      #
      # Update the state.
      #
      state$current <<- current
    }
    list(y=y, x=intervals)
  } # Update()

  process <- function(intervals, x, y) {
    # intervals is a list of 2-vectors. Each represents the endpoints of
    # an interior interval of a polygon.
    # x is an array of x-coordinates of vertices.
    #
    # Retains only the intervals containing at least one vertex.
    between <- function(i) {
      1 == max(mapply(function(a,b) a && b, i[1] <= x, x <= i[2]))
    }
    is.good <- lapply(intervals$x, between)
    list(y=y, x=intervals$x[unlist(is.good)])
    #intervals
  }
  #
  # Group the vertices by common y-coordinate.
  #
  vertices.x <- split(p$v[, "x"], p$v[, "y"])
  vertices.y <- lapply(split(p$v[, "y"], p$v[, "y"]), max)
  #
  # The "state" is a collection of segments and an index into edges.
  # It will updated during the vertical line sweep.
  #
  state <- list(level=-Inf, current=c(), i=1, x=c(), interior=c())
  #
  # Sweep vertically from bottom to top, processing the intersection
  # as we go.
  #
  mapply(function(x,y) process(update(y), x, y), vertices.x, vertices.y)
}


scale <- 10
p.raw = list(scale * cbind(x=c(0:10,7,6,0), y=c(3,0,0,-1,-1,-1,0,-0.5,0.75,1,4,1.5,0.5,3)),
             scale *cbind(x=c(1,1,2.4,2,4,4,4,4,2,1), y=c(0,1,2,1,1,0,-0.5,1,1,0)),
             scale *cbind(x=c(6,7,6,6), y=c(.5,2,3,.5)))

#p.raw = list(cbind(x=c(0,2,1,1/2,0), y=c(0,0,2,1,0)))
#p.raw = list(cbind(x=c(0, 35, 100, 65, 0), y=c(0, 50, 100, 50, 0)))

p <- as.polygon(p.raw)

results <- fetch.x(p)
#
# Find the longest.
#
dx <- matrix(unlist(results["x", ]), nrow=2)
length.max <- max(dx[2,] - dx[1,])
#
# Draw pictures.
#
segment.plot <- function(s, length.max, colors,  ...) {
  lapply(s$x, function(x) {
      col <- ifelse (diff(x) >= length.max, colors[1], colors[2])
      lines(x, rep(s$y,2), col=col, ...)
    })
}
gray <- "#f0f0f0"
grayer <- "#d0d0d0"
plot(expand(p$bb, 1.1), type="n", xlab="x", ylab="y", main="After the Scan")
sapply(1:length(p.raw), function(i) polygon(p.raw[[i]], col=c(gray, "White", grayer)[i]))
apply(results, 2, function(s) segment.plot(s, length.max, colors=c("Red", "#b8b8a8"), lwd=4))
plot(p, col="Black", lty=3)
points(p, pch=19, col=round(2 + 2*p$v[, "y"]/scale, 0))
points(p, cex=1.25)

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

Начните с растрового представления многоугольника. Этот использует сетку из 550 строк и 1200 столбцов:

В этом представлении серые (внутри) ячейки имеют значение 1, а все остальные ячейки - NoData.

Вычислите накопление потока в направлении запад-восток, используя значения элементарной ячейки для весовой сетки (количество «осадков»):

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

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

Красные клетки отмечают концы самых высоких скоплений потока: они являются самыми правыми конечными точками внутренних сегментов максимальной длины многоугольника.

Чтобы найти эти сегменты, поместите весь вес на красные клетки и запустите поток назад!

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

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

Одним из действительно приятных аспектов этого подхода является то, что обычно мы находим экстремальные ценности вещей как часть более широкого рабочего процесса, в котором необходимо достичь какой-то цели: размещения трубопровода или футбольного поля, создания экологических буферов и так далее. Процесс включает в себя компромиссы. Таким образом, самая длинная горизонтальная линия не может быть частью оптимального решения. Мы могли бы позаботиться о том, чтобы знать, где лежат почти самые длинные строки. Это просто: вместо выбора зонального максимального потока выберите все ячейки, близкие к зональному максимуму. В этом примере зональный максимум равен 744 (количество столбцов, охватываемых самым длинным внутренним сегментом). Вместо этого давайте выделим все ячейки в пределах 5% от этого максимума:

Запуск потока с востока на запад производит эту коллекцию горизонтальных сегментов:

Это карта местоположений, где непрерывная протяженность восток-запад составляет 95% или больше максимальной протяженности восток-запад в любом месте полигона.

Ok. У меня есть другая (лучшая) идея ( идея №2 ). Но я полагаю, что его лучше реализовывать как скрипт на python, а не как SQL-запрос. Опять же, это общий случай, а не только EW.

Вам понадобится ограничительная рамка для многоугольника и азимут (A) в качестве направления измерения. Предположим, что длина ребер BBox равна LA и LB. Максимально возможное расстояние (MD) в пределах многоугольника: MB = (LA^2 * LB^2)^(1/2)так ищет значение (V) не больше , чем МБ: V <= MB.

1. Начиная с любой вершины BBox, создайте линию (LL) с длиной MB и азимутом A.
2. Пересечь линию LL с полигоном, чтобы получить линию пересечения (IL)
3. Проверьте геометрию IL - если в линии IL только две точки, рассчитайте ее длину. Если 4 или более - подсчитать отрезки и подобрать длину самого длинного. Нуль (пересечения нет вообще) - игнорировать.
4. Продолжайте создавать новые линии LL, двигаясь от начальной точки против часовой стрелки или по часовой стрелке к краям BBox, пока вы не окажетесь в начальной точке (вы будете делать весь цикл через BBox).
5. Подберите самое большое значение длины IL (на самом деле вам не нужно хранить все длины, вы можете сохранять максимальное значение «пока» при цикле) - это будет то, что вы ищете.

Я не уверен, что ответ Фетцера - это то, что вы хотите сделать, но именно так st_box2d может сделать эту работу.

Идея SS_Rebelious № 1 будет работать во многих случаях, но не для некоторых вогнутых многоугольников.

Я думаю, что вы должны создать искусственные линии lw, точки которых следуют за ребрами, когда линии вершин пересекают границы многоугольника, если существует возможность линии восток-запад.

Для этого вы можете попытаться создать многоугольник из 4 узлов с большой длиной линии, сделать многоугольник P *, который является предыдущим, перекрывающимся с вашим исходным многоугольником, и посмотреть, оставляют ли min (y1) и max (y2) некоторую x-линию возможность. (где y1 - набор точек между верхним левым корнетом и верхним правым углом, а y2 - набор y между нижним левым и нижним правыми углами полигона из 4 узлов). Это не так просто, я надеюсь, вы найдете psql инструменты, которые помогут вам!

У меня есть идея-№1 ( Правка: для общего случая, не только для направления ВП, но и с некоторыми ограничениями, которые описаны в комментариях). Я не буду предоставлять код, просто концепцию. «Направление x» на самом деле является азимутом, который вычисляется с помощью ST_Azimuth. Предлагаемые шаги:

1. Извлеките все вершины из многоугольника в виде точек.
2. Создайте линии между каждой парой точек.
3. Выделите линии (назовем их lw-линиями), которые находятся внутри исходного многоугольника (нам не нужны линии, которые будут пересекать границы многоугольника).
4. Найти расстояния и азимуты для каждой линии lw.
5. Выберите самое длинное расстояние от линий lw, где азимут равен искуемому азимуту или лежит в некотором интервале (может быть, что ни один азимут не будет точно равен искуемому азимуту).

Взгляните на мой вопрос и ответ от Evil Genius.

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

Вот пример:

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

Идея для неправильных или вогнутых многоугольников для этого подхода. Может быть, вы должны пересечь многоугольник с растром.