Calculando a distância máxima dentro do polígono na direção x (direção leste-oeste) no PostGIS?

Estou interessado na largura máxima de um polígono, por exemplo, um lago, na direção leste-oeste. Caixas delimitadoras ajudarão apenas em polígonos simples, mas não em polígonos côncavos complexos.

Provavelmente, isso requer algum script em qualquer plataforma GIS.

O método mais eficiente (assintoticamente) é uma varredura de linha vertical: requer a classificação das arestas pelas coordenadas y mínimas e o processamento das arestas de baixo (y mínimo) para cima (y máximo), para um O (e * log ( e)) algoritmo quando e arestas estão envolvidas.

O procedimento, embora simples, é surpreendentemente complicado de acertar em todos os casos. Os polígonos podem ser desagradáveis: podem ter oscilações, lascas, orifícios, serem desconectados, ter vértices duplicados, trechos de vértices ao longo de linhas retas e ter limites não dissolvidos entre dois componentes adjacentes. Aqui está um exemplo que exibe muitas dessas características (e mais):

Buscaremos especificamente o (s) segmento (s) horizontal (is) de comprimento máximo (s) inteiramente dentro do fechamento do polígono. Por exemplo, isso elimina a oscilação entre x = 20 ex = 40 emanando do furo entre x = 10 ex = 25. É então simples mostrar que pelo menos um dos segmentos horizontais de comprimento máximo cruza pelo menos um vértice. (Se não houver soluções cruzam há vértices, eles vão deitar no interior de alguns paralelogramo delimitada na parte superior e inferior por soluções que fazer se cruzam pelo menos um vértice. Isso nos dá um meio de encontrar todas as soluções.)

Assim, a varredura de linha deve começar com os vértices mais baixos e depois mover para cima (ou seja, em direção a valores mais altos de y) para parar em cada vértice. Em cada parada, encontramos novas arestas que emanam para cima dessa elevação; elimine quaisquer arestas que terminem de baixo nessa elevação (esta é uma das idéias principais: simplifica o algoritmo e elimina metade do processamento potencial); e processe cuidadosamente todas as arestas totalmente em elevação constante (as arestas horizontais).

Por exemplo, considere o estado quando um nível de y = 10 for atingido. Da esquerda para a direita, encontramos as seguintes arestas:

      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

Nesta tabela, (x.min, y.min) são coordenadas do ponto final inferior da aresta e (x.max, y.max) são coordenadas do seu ponto final superior. Nesse nível (y = 10), a primeira aresta é interceptada em seu interior, a segunda é interceptada em sua base e assim por diante. Algumas arestas que terminam nesse nível, como (10,0) a (10,10), não estão incluídas na lista.

Para determinar onde estão os pontos internos e externos, imagine começar da extrema esquerda - que está fora do polígono, é claro - e se mover horizontalmente para a direita. Cada vez que cruzamos uma borda que não é horizontal , alternamos alternadamente do exterior para o interior e para trás. (Essa é outra ideia importante.) No entanto, todos os pontos em qualquer aresta horizontal são determinados como estando dentro do polígono, não importa o quê. (O fechamento de um polígono sempre inclui suas bordas.)

Continuando o exemplo, aqui está a lista classificada de coordenadas x em que arestas não horizontais começam na linha y = 10 ou cruzam:

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

(Observe que x = 40 não está nesta lista.) Os valores da interiormatriz marcam os pontos finais esquerdos dos segmentos internos: 1 designa um intervalo interno, 0 um intervalo externo. Assim, o primeiro 1 indica que o intervalo de x = 6,7 ex = 10 está dentro do polígono. O próximo 0 indica que o intervalo de x = 10 ex = 20 está fora do polígono. E assim prossegue: a matriz identifica quatro intervalos separados como dentro do polígono.

Alguns desses intervalos, como o de x = 60 ex = 63,3, não cruzam nenhum vértice: uma verificação rápida contra as coordenadas x de todos os vértices com y = 10 elimina esses intervalos.

Durante a varredura, podemos monitorar os comprimentos desses intervalos, mantendo os dados referentes aos intervalos de comprimento máximo encontrados até o momento.

Observe algumas das implicações dessa abordagem. Um vértice em forma de "v", quando encontrado, é a origem de duas arestas. Portanto, dois comutadores ocorrem ao cruzá-lo. Esses interruptores são cancelados. Qualquer "v" invertido nem é processado, porque as duas arestas são eliminadas antes de iniciar a digitalização da esquerda para a direita. Nos dois casos, esse vértice não bloqueia um segmento horizontal.

Mais de duas arestas podem compartilhar um vértice: isso é ilustrado em (10,0), (60,5), (25, 20) e - embora seja difícil dizer - em (20,10) e (40 10). (Isso ocorre porque o dangle vai (20,10) -> (40,10) -> (40,0) -> (40, -50) -> (40, 10) -> (20, 10) Observe como o vértice em (40,0) também está no interior de outra aresta ... isso é desagradável.) Esse algoritmo lida bem com essas situações.

Uma situação complicada é ilustrada bem no fundo: as coordenadas x de segmentos não horizontais existem

30, 50

Isso faz com que tudo à esquerda de x = 30 seja considerado exterior, tudo entre 30 e 50 seja interior e tudo depois de 50 seja exterior novamente. O vértice em x = 40 nunca é considerado neste algoritmo.

Aqui está a aparência do polígono no final da verificação. Mostro todos os intervalos interiores contendo vértices em cinza escuro, todos os intervalos de comprimento máximo em vermelho e cor os vértices de acordo com as coordenadas y. O intervalo máximo é de 64 unidades.

Os únicos cálculos geométricos envolvidos são calcular onde as bordas cruzam as linhas horizontais: isso é uma interpolação linear simples. Também são necessários cálculos para determinar quais segmentos interiores contêm vértices: essas são determinações de intermediação , facilmente calculadas com algumas desigualdades. Essa simplicidade torna o algoritmo robusto e apropriado para representações de coordenadas de números inteiros e de ponto flutuante.

Se as coordenadas são geográficas , as linhas horizontais estão realmente em círculos de latitude. Seus comprimentos não são difíceis de calcular: basta multiplicar seus comprimentos euclidianos pelo cosseno de sua latitude (em um modelo esférico). Portanto, esse algoritmo se adapta bem às coordenadas geográficas. (Para lidar com o envolvimento do poço do meridiano + -180, talvez seja necessário primeiro encontrar uma curva do polo sul ao polo norte que não passe pelo polígono. Depois de re-expressar todas as coordenadas x como deslocamentos horizontais em relação àquele curva, este algoritmo encontrará corretamente o segmento horizontal máximo.)

A seguir, o Rcódigo implementado para executar os cálculos e criar as ilustrações.

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

Aqui está uma solução baseada em varredura. É rápido (fiz todo o trabalho do início ao fim em 14 minutos), não requer scripts, leva apenas algumas operações e é razoavelmente preciso.

Comece com uma representação raster do polígono. Este usa uma grade de 550 linhas e 1200 colunas:

Nesta representação, as células cinza (internas) têm o valor 1 e todas as outras células são NoData.

Calcule a acumulação de fluxo na direção oeste-leste usando os valores das células unitárias para a grade de peso (quantidade de "precipitação"):

A baixa acumulação é escura, aumentando para as acumulações mais altas no amarelo brilhante.

Um máximo zonal (usando o polígono para a grade e a acumulação de fluxo para os valores) identifica a (s) célula (s) em que o fluxo é maior. Para mostrar isso, tive que aumentar o zoom no canto inferior direito:

As células vermelhas marcam as extremidades das acumulações mais altas de fluxo: elas são as extremidades mais à direita dos segmentos interiores de comprimento máximo do polígono.

Para encontrar esses segmentos, coloque todo o peso nas células vermelhas e execute o fluxo para trás!

A faixa vermelha próxima à parte inferior marca duas linhas de células: dentro delas está o segmento horizontal de comprimento máximo. Use essa representação como está para análises adicionais ou converta-a em uma forma de polilinha (ou polígono).

Há algum erro de discretização feito com uma representação raster. Pode ser reduzido aumentando a resolução, a algum custo no tempo de computação.

Um aspecto realmente interessante dessa abordagem é que geralmente encontramos valores extremos das coisas como parte de um fluxo de trabalho maior, no qual algum objetivo precisa ser alcançado: localização de um oleoduto ou campo de futebol, criação de amortecedores ecológicos e assim por diante. O processo envolve trocas. Portanto, a linha horizontal mais longa pode não fazer parte de uma solução ideal. Podemos cuidar em vez de saber onde quase linhas mais longas iria mentir. Isso é simples: em vez de selecionar o fluxo máximo zonal, selecione todas as células próximas a um máximo zonal. Neste exemplo, o zonal max é igual a 744 (o número de colunas estendidas pelo segmento interior mais longo). Em vez disso, vamos selecionar todas as células dentro de 5% desse máximo:

A execução do fluxo de leste a oeste produz esta coleção de segmentos horizontais:

Este é um mapa de locais onde a extensão leste-oeste ininterrupta é 95% ou maior que a extensão leste-oeste máxima em qualquer lugar dentro do polígono.

Está bem. Eu tenho outra (melhor) ideia ( idéia №2 ). Mas suponho que seja melhor ser realizado como um script python, não como uma consulta SQL. Novamente, aqui está o caso comum, não apenas o EW.

Você precisará de uma caixa delimitadora para o polígono e um azimute (A) como sua direção de medição. Suponha que o comprimento das arestas da BBox seja LA e LB. A distância máxima possível (MD) dentro de um polígono é:, MB = (LA^2 * LB^2)^(1/2)portanto, a busca do valor (V) não é maior que MB:V <= MB .

1. A partir de qualquer vértice do BBox, crie uma linha (LL) com o comprimento MB e o azimute A.
2. Interseção da linha LL com o polígono para obter a linha de interseção (IL)
3. Verifique a geometria da IL - se houver apenas dois pontos na linha IL, calcule seu comprimento. Se 4 ou mais - calcule os segmentos e escolha o comprimento do mais longo. Nulo (sem interseção) - ignore.
4. Continue criando outras linhas LL movendo-se do contador do ponto inicial ou no sentido horário em direção às bordas do BBox até que você não termine no ponto inicial (você fará todo o loop no BBox).
5. Escolha o maior valor de comprimento de IL (na verdade você não precisa armazenar todos os comprimentos, você pode manter o valor máximo 'até agora' durante o loop) - será o que você procura.

Não tenho certeza de que a resposta de Fetzer seja o que você deseja fazer, mas o st_box2d pode fazer o trabalho.

A idéia de SS_Rebelious N ° 1 funcionará em muitos casos, mas não em alguns polígonos côncavos.

Penso que você deve criar linhas lw artificiais cujos pontos seguem arestas quando linhas feitas por vértices cruzam as bordas do polígono se houver uma possibilidade de linha leste-oeste.

Para isso, você pode tentar fazer um polígono de 4 nós em que o comprimento da linha é alto, criar o polígono P *, que é a sobreposição anterior do polígono original e ver se os min (y1) e max (y2) deixam alguma linha x possibilidade. (onde y1 é o conjunto de pontos entre o canto superior esquerdo e o canto superior direito e y2 o conjunto de y entre os cantos inferior esquerdo e inferior direito do polígono de 4 nós). Não é tão fácil assim, espero que você encontre ferramentas psql para ajudá-lo!

Eu tenho uma idéia №1 ( Edit: para casos comuns, não apenas a direção EW, e com algumas limitações descritas nos comentários). Não fornecerei o código, apenas um conceito. A "direção x" é na verdade um azimute, calculado por ST_Azimuth. As etapas propostas são:

1. Extraia todos os vértices do polígono como pontos.
2. Crie linhas entre cada par de pontos.
3. Selecione linhas (vamos chamá-las de linhas lw) que estão dentro do polígono original (não precisamos de linhas que cruzem as bordas do polígono).
4. Encontre distâncias e azimutes para cada linha lw.
5. Selecione a distância mais longa das linhas lw em que o azimute é igual ao procurado ou fica em algum intervalo (pode ser que nenhum azimute seja exatamente igual ao procurado).

Dê uma olhada na minha pergunta e na resposta do Evil Genius.

Esperamos que o polígono do seu lago tenha vários pontos; você pode criar linhas nesses pontos com um azimute (aspecto, direção geográfica). Escolha o comprimento das linhas grandes o suficiente (parte ST_MakePoint), para poder calcular a linha mais curta entre as duas linhas mais distantes.

Aqui está um exemplo:

O exemplo mostra a largura máxima do polígono. Eu escolho ST_ShortestLine (linha vermelha) para esta abordagem. ST_MakeLine aumentaria o valor (linha azul) e o ponto final da linha (canto inferior esquerdo) atingiria a linha azul do polígono. Você deve calcular a distância com os centróides das linhas criadas (de ajuda).

Uma idéia para polígonos irregulares ou côncavos para essa abordagem. Pode ser que você tenha que cruzar o polígono com uma varredura.