Cálculo de la distancia máxima dentro del polígono en dirección x (dirección este-oeste) en PostGIS?

Estoy interesado en el ancho máximo de un polígono, por ejemplo, un lago, en dirección este-oeste. Los cuadros delimitadores ayudarán solo en polígonos simples pero no en polígonos cóncavos complejos.

Esto probablemente requiere algunas secuencias de comandos en cualquier plataforma SIG.

El método más eficiente (asintóticamente) es un barrido de línea vertical: requiere ordenar los bordes por sus coordenadas mínimas y y luego procesar los bordes de abajo (mínimo y) a arriba (máximo y), para un O (e * log ( e)) algoritmo cuando están involucrados los bordes e .

El procedimiento, aunque simple, es sorprendentemente difícil de acertar en todos los casos. Los polígonos pueden ser desagradables: pueden tener colgantes, astillas, agujeros, estar desconectados, tener vértices duplicados, tramos de vértices a lo largo de líneas rectas y tener límites no disueltos entre dos componentes adyacentes. Aquí hay un ejemplo que exhibe muchas de estas características (y más):

Buscaremos específicamente los segmentos horizontales de longitud máxima que se encuentran completamente dentro del cierre del polígono. Por ejemplo, esto elimina la oscilación entre x = 20 yx = 40 que emana del agujero entre x = 10 yx = 25. Entonces es sencillo mostrar que al menos uno de los segmentos horizontales de longitud máxima se cruza con al menos un vértice. (Si no hay soluciones que se cruzan no hay vértices, se encuentran en el interior de algunos de paralelogramo delimitado en la parte superior e inferior mediante soluciones que hacen intersecan al menos un vértice. Esto nos da un medio para encontrar todas las soluciones.)

En consecuencia, el barrido de línea debe comenzar con los vértices más bajos y luego moverse hacia arriba (es decir, hacia valores de y más altos) para detenerse en cada vértice. En cada parada, encontramos nuevos bordes que emanan hacia arriba desde esa elevación; elimine cualquier borde que termine desde abajo en esa elevación (esta es una de las ideas clave: simplifica el algoritmo y elimina la mitad del procesamiento potencial); y procese cuidadosamente cualquier borde que se encuentre totalmente a una elevación constante (los bordes horizontales).

Por ejemplo, considere el estado cuando se alcanza un nivel de y = 10. De izquierda a derecha, encontramos los siguientes bordes:

      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

En esta tabla, (x.min, y.min) son coordenadas del punto final inferior del borde y (x.max, y.max) son coordenadas de su punto final superior. En este nivel (y = 10), el primer borde se intercepta dentro de su interior, el segundo se intercepta en su parte inferior, y así sucesivamente. Algunos bordes que terminan en este nivel, como de (10,0) a (10,10), no están incluidos en la lista.

Para determinar dónde están los puntos interiores y los exteriores, imagine comenzar desde el extremo izquierdo, que está fuera del polígono, por supuesto, y moverse horizontalmente hacia la derecha. Cada vez que cruzamos un borde que no es horizontal , alternadamente cambiamos de exterior a interior y viceversa. (Esta es otra idea clave). Sin embargo, se determina que todos los puntos dentro de cualquier borde horizontal están dentro del polígono, pase lo que pase. (El cierre de un polígono siempre incluye sus bordes).

Continuando con el ejemplo, aquí está la lista ordenada de coordenadas x donde los bordes no horizontales comienzan o cruzan la línea y = 10:

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 no está en esta lista). Los valores de la interiormatriz marcan los puntos finales izquierdos de los segmentos interiores: 1 designa un intervalo interior, 0 un intervalo exterior. Por lo tanto, el primer 1 indica que el intervalo de x = 6.7 a x = 10 está dentro del polígono. El siguiente 0 indica que el intervalo de x = 10 a x = 20 está fuera del polígono. Y así continúa: la matriz identifica cuatro intervalos separados como dentro del polígono.

Algunos de estos intervalos, como el de x = 60 a x = 63.3, no se cruzan con ningún vértice: una verificación rápida contra las coordenadas x de todos los vértices con y = 10 elimina dichos intervalos.

Durante el escaneo podemos monitorear las longitudes de estos intervalos, reteniendo los datos relativos a los intervalos de longitud máxima encontrados hasta ahora.

Observe algunas de las implicaciones de este enfoque. Un vértice en forma de "v", cuando se encuentra, es el origen de dos bordes. Por lo tanto, se producen dos interruptores al cruzarlo. Esos interruptores se cancelan. Cualquier "v" invertida ni siquiera se procesa, porque sus dos bordes se eliminan antes de comenzar el escaneo de izquierda a derecha. En ambos casos, dicho vértice no bloquea un segmento horizontal.

Más de dos bordes pueden compartir un vértice: esto se ilustra en (10,0), (60,5), (25, 20) y, aunque es difícil de distinguir, en (20,10) y (40 , 10). (Eso se debe a que el cuelga va (20,10) -> (40,10) -> (40,0) -> (40, -50) -> (40, 10) -> (20, 10). Observe cómo el vértice en (40,0) también está en el interior de otro borde ... eso es desagradable.) Este algoritmo maneja esas situaciones muy bien.

Una situación difícil se ilustra en la parte inferior: las coordenadas x de los segmentos no horizontales que hay

30, 50

Esto hace que todo a la izquierda de x = 30 se considere exterior, que todo entre 30 y 50 sea interior, y que todo después de 50 sea exterior nuevamente. El vértice en x = 40 nunca se considera en este algoritmo.

Así es como se ve el polígono al final del escaneo. Muestro todos los intervalos interiores que contienen vértices en gris oscuro, cualquier intervalo de longitud máxima en rojo, y coloreo los vértices de acuerdo con sus coordenadas y. El intervalo máximo es de 64 unidades de largo.

Los únicos cálculos geométricos involucrados son calcular dónde los bordes se cruzan con las líneas horizontales: es una simple interpolación lineal. También se necesitan cálculos para determinar qué segmentos interiores contienen vértices: estas son determinaciones de intermediación , que se calculan fácilmente con un par de desigualdades. Esta simplicidad hace que el algoritmo sea robusto y apropiado tanto para representaciones de coordenadas de punto entero como flotante.

Si las coordenadas son geográficas , entonces las líneas horizontales están realmente en círculos de latitud. Sus longitudes no son difíciles de calcular: simplemente multiplique sus longitudes euclidianas por el coseno de su latitud (en un modelo esférico). Por lo tanto, este algoritmo se adapta muy bien a las coordenadas geográficas. (Para manejar la envoltura alrededor del pozo del meridiano + -180, primero podría ser necesario encontrar una curva desde el polo sur hasta el polo norte que no pase por el polígono. Después de volver a expresar todas las coordenadas x como desplazamientos horizontales en relación con ese curva, este algoritmo encontrará correctamente el segmento horizontal máximo).

El siguiente es el Rcódigo implementado para realizar los cálculos y crear las ilustraciones.

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

Aquí hay una solución basada en ráster. Es rápido (hice todo el trabajo de principio a fin en 14 minutos), no requiere secuencias de comandos, solo requiere algunas operaciones y es razonablemente preciso.

Comience con una representación ráster del polígono. Éste usa una cuadrícula de 550 filas y 1200 columnas:

En esta representación, las celdas grises (dentro) tienen el valor 1 y todas las demás celdas son NoData.

Calcule la acumulación de flujo en dirección oeste a este utilizando los valores de celda unitaria para la cuadrícula de peso (cantidad de "lluvia"):

La acumulación baja es oscura, aumentando a las acumulaciones más altas en el amarillo brillante.

Un máximo zonal (que usa el polígono para la cuadrícula y la acumulación de flujo para los valores) identifica las celdas donde el flujo es mayor. Para mostrar esto, tuve que acercarme a la esquina inferior derecha:

Las células rojas marcan los extremos de las acumulaciones más altas de flujo: son los puntos finales más a la derecha de los segmentos interiores de longitud máxima del polígono.

Para encontrar estos segmentos, coloque todo el peso en los glóbulos rojos y ejecute el flujo hacia atrás.

La franja roja cerca de la parte inferior marca dos filas de celdas: dentro de ellas se encuentra el segmento horizontal de longitud máxima. Utilice esta representación tal cual para su posterior análisis o conviértala en una forma de polilínea (o polígono).

Se produce un error de discretización con una representación ráster. Se puede reducir aumentando la resolución, a un costo en tiempo de cálculo.

Un aspecto realmente agradable de este enfoque es que normalmente encontramos valores extremos de las cosas como parte de un flujo de trabajo más amplio en el que se debe alcanzar algún objetivo: ubicar una tubería o campo de fútbol, ​​crear amortiguadores ecológicos, etc. El proceso implica compensaciones. Por lo tanto, la línea horizontal más larga podría no ser parte de una solución óptima. Podríamos cuidar su lugar para saber dónde casi se encontrarían líneas más largas. Esto es simple: en lugar de seleccionar el flujo máximo zonal, seleccione todas las celdas cercanas a un máximo zonal. En este ejemplo, el máximo zonal es igual a 744 (el número de columnas abarcadas por el segmento interior más largo). En su lugar, seleccionemos todas las celdas dentro del 5% de ese máximo:

Ejecutar el flujo de este a oeste produce esta colección de segmentos horizontales:

Este es un mapa de ubicaciones donde la extensión ininterrumpida este-oeste es 95% o mayor que la extensión máxima este-oeste en cualquier lugar dentro del polígono.

Okay. Tengo otra (mejor) idea ( idea-№2 ). Pero supongo que es mejor realizarse como un script de Python, no como una consulta SQL. Una vez más, este es el caso común, no solo EW.

Necesitará un cuadro delimitador para el polígono y un acimut (A) como dirección de medición. Suponga que la longitud de los bordes de BBox son LA y LB. La máxima distancia posible (MD) dentro de un polígono es: MB = (LA^2 * LB^2)^(1/2), por lo que la búsqueda de valor (V) no es más grande que MB: V <= MB.

1. Comenzando desde cualquier vértice del BBox, cree una línea (LL) con la longitud MB y azimut A.
2. Intersecte la línea LL con el polígono para obtener la línea de intersección (IL)
3. Verifique la geometría de IL: si solo hay dos puntos en la línea IL, calcule su longitud. Si tiene 4 o más, calcule segmentos y tome la longitud del más largo. Nulo (sin intersección): ignorar.
4. Siga creando otras líneas LL moviéndose desde el punto de inicio contrario o en el sentido de las agujas del reloj hacia los bordes del BBox hasta que no termine en el punto de inicio (hará todo el bucle a través de BBox).
5. Elija el valor de longitud IL más grande (en realidad no es necesario que almacene todas las longitudes, puede mantener el valor máximo "hasta ahora" durante el bucle): será lo que busca.

No estoy seguro de que la respuesta de Fetzer sea lo que quieres hacer, pero es así que st_box2d puede hacer el trabajo.

La idea N ° 1 de SS_Rebelious funcionará en muchos casos, pero no para algunos polígonos cóncavos.

Creo que debe crear líneas lw artificiales cuyos puntos sigan a los bordes cuando las líneas de vértice cruzan los bordes del polígono si existe una posibilidad de línea este-oeste.

Para esto, puede intentar hacer un polígono de 4 nodos donde la longitud de la línea sea alta, hacer que el polígono P *, que es el anterior, se superponga con su polígono original, y ver si min (y1) y max (y2) dejan algo de línea x posibilidad. (donde y1 es el conjunto de puntos entre la esquina superior izquierda y la esquina superior derecha y y2 el conjunto de y entre las esquinas inferior izquierda e inferior derecha de su polígono de 4 nodos). ¡No es tan fácil, espero que encuentres herramientas psql para ayudarte!

Tengo una idea-№1 ( Editar: para el caso común, no solo la dirección EW, y con algunas limitaciones que se describen en los comentarios). No proporcionaré el código, solo un concepto. La "dirección x" es en realidad un acimut, que se calcula mediante ST_Azimuth. Los pasos propuestos son:

1. Extrae todos los vértices del polígono como puntos.
2. Crea líneas entre cada par de puntos.
3. Seleccione líneas (llamémoslas lw-lines) que están dentro del polígono original (no necesitamos líneas que crucen los bordes del polígono).
4. Encuentra distancias y acimutes para cada línea lw.
5. Seleccione la distancia más larga desde las líneas lw donde el acimut es igual al acimut buscado o se encuentra en algún intervalo (puede ser que ningún acimut sea exactamente igual al acimut buscado).

Echa un vistazo a mi pregunta y la respuesta de Evil Genius.

Esperemos que su polígono de lago tenga varios puntos, puede crear líneas en estos puntos con un acimut (aspecto, dirección geográfica). Elija la longitud de las líneas lo suficientemente grande (parte ST_MakePoint), para que pueda calcular la línea más corta entre las dos líneas más distantes.

Aquí hay un ejemplo:

El ejemplo muestra el ancho máximo del polígono. Elijo ST_ShortestLine (línea roja) para este enfoque. ST_MakeLine aumentaría el valor (línea azul) y el punto final de la línea (abajo a la izquierda) golpearía la línea azul del polígono. Debe calcular la distancia con los centroides de las líneas creadas (ayuda).

Una idea para polígonos irregulares o cóncavos para este enfoque. Es posible que tenga que intersecar el polígono con un ráster.