Comment calculer la distance maximale à l'intérieur d'un polygone dans la direction x (direction est-ouest) dans PostGIS?

Je m'intéresse à la largeur maximale d'un polygone, par exemple un lac, dans la direction est-ouest. Les boîtes de délimitation ne seront utiles que dans les polygones simples mais pas dans les polygones concaves complexes.

Cela nécessite probablement des scripts dans n'importe quelle plate-forme SIG.

La méthode la plus efficace (asymptotiquement) est un balayage de ligne verticale: elle nécessite de trier les bords par leurs coordonnées y minimales puis de traiter les bords de bas (minimum y) à haut (maximum y), pour un O (e * log ( e)) algorithme lorsque les bords e sont impliqués.

La procédure, bien que simple, est étonnamment difficile à obtenir dans tous les cas. Les polygones peuvent être désagréables: ils peuvent avoir des pendants, des éclats, des trous, être déconnectés, avoir des sommets dupliqués, des séries de sommets le long de lignes droites et avoir des limites non dissoutes entre deux composants adjacents. Voici un exemple présentant plusieurs de ces caractéristiques (et plus):

Nous rechercherons spécifiquement le ou les segments horizontaux de longueur maximale entièrement situés à l'intérieur de la fermeture du polygone. Par exemple, cela élimine le balancement entre x = 20 et x = 40 émanant du trou entre x = 10 et x = 25. Il est alors simple de montrer qu'au moins un des segments horizontaux de longueur maximale coupe au moins un sommet. (S'il existe des solutions qui se croisent pas les sommets, ils se trouvent à l'intérieur de quelque parallélogramme délimité en haut et en bas par des solutions qui font Intersection au moins un sommet. Cela nous donne un moyen de trouver toutes les solutions.)

Par conséquent, le balayage de ligne doit commencer par les sommets les plus bas, puis se déplacer vers le haut (c'est-à-dire vers des valeurs y plus élevées) pour s'arrêter à chaque sommet. À chaque arrêt, nous trouvons de nouvelles arêtes émanant vers le haut de cette élévation; éliminer tous les bords se terminant par le bas à cette élévation (c'est l'une des idées clés: cela simplifie l'algorithme et élimine la moitié du traitement potentiel); et traitez soigneusement tous les bords se trouvant entièrement à une altitude constante (les bords horizontaux).

Par exemple, considérons l'état lorsqu'un niveau de y = 10 est atteint. De gauche à droite, on retrouve les bords suivants:

      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

Dans ce tableau, (x.min, y.min) sont les coordonnées de l'extrémité inférieure de l'arête et (x.max, y.max) sont les coordonnées de son extrémité supérieure. A ce niveau (y = 10), le premier bord est intercepté à l'intérieur, le second est intercepté en bas, etc. Certaines arêtes se terminant à ce niveau, comme de (10,0) à (10,10), ne sont pas incluses dans la liste.

Pour déterminer où se trouvent les points intérieurs et ceux extérieurs, imaginez en partant de l'extrême gauche - qui est en dehors du polygone, bien sûr - et en vous déplaçant horizontalement vers la droite. Chaque fois que nous traversons un bord qui n'est pas horizontal , nous basculons alternativement de l'extérieur vers l'intérieur et le dos. (C'est une autre idée clé.) Cependant, tous les points à l'intérieur d'une arête horizontale sont déterminés comme étant à l'intérieur du polygone, quoi qu'il arrive. (La fermeture d'un polygone inclut toujours ses bords.)

Poursuivant l'exemple, voici la liste triée des coordonnées x où les bords non horizontaux commencent à ou traversent la ligne y = 10:

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

(Notez que x = 40 ne figure pas dans cette liste.) Les valeurs du interiortableau marquent les extrémités gauche des segments intérieurs: 1 désigne un intervalle intérieur, 0 un intervalle extérieur. Ainsi, le premier 1 indique que l'intervalle de x = 6,7 à x = 10 est à l'intérieur du polygone. Le 0 suivant indique que l'intervalle de x = 10 à x = 20 est en dehors du polygone. Et donc ça continue: le tableau identifie quatre intervalles distincts comme à l'intérieur du polygone.

Certains de ces intervalles, tels que celui de x = 60 à x = 63,3, ne coupent aucun sommet: une vérification rapide des coordonnées x de tous les sommets avec y = 10 élimine ces intervalles.

Pendant le balayage, nous pouvons surveiller la longueur de ces intervalles, en conservant les données concernant le ou les intervalles de longueur maximale trouvés jusqu'à présent.

Notez certaines des implications de cette approche. Un sommet en "V", lorsqu'il est rencontré, est à l'origine de deux arêtes. Par conséquent, deux commutations se produisent lors de la traversée. Ces commutateurs s'annulent. Tout «v» inversé n'est même pas traité, car ses deux bords sont éliminés avant de lancer le balayage de gauche à droite. Dans les deux cas, un tel sommet ne bloque pas un segment horizontal.

Plus de deux arêtes peuvent partager un sommet: cela est illustré en (10,0), (60,5), (25, 20) et - bien que ce soit difficile à dire - en (20,10) et (40 ,dix). (C'est parce que le balancement va (20,10) -> (40,10) -> (40,0) -> (40, -50) -> (40, 10) -> (20, Remarquez comment le sommet en (40,0) est également à l'intérieur d'un autre bord ... c'est méchant.) Cet algorithme gère très bien ces situations.

Une situation délicate est illustrée tout en bas: les coordonnées x des segments non horizontaux sont

30, 50

Ainsi, tout ce qui se trouve à gauche de x = 30 est considéré comme extérieur, tout ce qui se situe entre 30 et 50 est intérieur et tout ce qui se trouve après 50 est de nouveau extérieur. Le sommet à x = 40 n'est même jamais considéré dans cet algorithme.

Voici à quoi ressemble le polygone à la fin du scan. Je montre tous les intervalles intérieurs contenant des sommets en gris foncé, tous les intervalles de longueur maximale en rouge et colorie les sommets en fonction de leurs coordonnées y. L'intervalle maximum est de 64 unités.

Les seuls calculs géométriques impliqués sont de calculer où les arêtes coupent les lignes horizontales: c'est une simple interpolation linéaire. Des calculs sont également nécessaires pour déterminer quels segments intérieurs contiennent des sommets: ce sont des déterminations d' interdépendance , facilement calculées avec quelques inégalités. Cette simplicité rend l'algorithme robuste et approprié à la fois pour les représentations de coordonnées entières et à virgule flottante.

Si les coordonnées sont géographiques , alors les lignes horizontales sont vraiment sur des cercles de latitude. Leurs longueurs ne sont pas difficiles à calculer: il suffit de multiplier leurs longueurs euclidiennes par le cosinus de leur latitude (dans un modèle sphérique). Par conséquent, cet algorithme s'adapte bien aux coordonnées géographiques. (Pour gérer l'habillage autour du puits méridien + -180, il peut être nécessaire de trouver d'abord une courbe du pôle sud au pôle nord qui ne traverse pas le polygone. Après avoir ré-exprimé toutes les coordonnées x en déplacements horizontaux par rapport à celui courbe, cet algorithme trouvera correctement le segment horizontal maximum.)

Voici le Rcode implémenté pour effectuer les calculs et créer les illustrations.

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

Voici une solution basée sur un raster. C'est rapide (j'ai fait tout le travail du début à la fin en 14 minutes), ne nécessite aucun script, ne prend que quelques opérations et est assez précis.

Commencez avec une représentation raster du polygone. Celui-ci utilise une grille de 550 lignes et 1200 colonnes:

Dans cette représentation, les cellules grises (à l'intérieur) ont la valeur 1 et toutes les autres cellules sont NoData.

Calculez l'accumulation de flux dans le sens ouest-est en utilisant les valeurs des cellules unitaires pour la grille de poids (quantité de "précipitations"):

Une faible accumulation est sombre, augmentant aux accumulations les plus élevées dans le jaune vif.

Un maximum zonal (en utilisant le polygone pour la grille et l'accumulation de flux pour les valeurs) identifie la ou les cellules où le flux est le plus important. Pour les montrer, j'ai dû zoomer en bas à droite:

Les globules rouges marquent les extrémités des accumulations de flux les plus élevées: ce sont les extrémités les plus à droite des segments intérieurs de longueur maximale du polygone.

Pour trouver ces segments, placez tout le poids sur les globules rouges et faites couler le flux vers l'arrière!

La bande rouge près du bas marque deux rangées de cellules: à l'intérieur de celles-ci se trouve le segment horizontal de longueur maximale. Utilisez cette représentation telle quelle pour une analyse plus approfondie ou convertissez-la en forme de polyligne (ou de polygone).

Il y a une erreur de discrétisation faite avec une représentation raster. Il peut être réduit en augmentant la résolution, à un certain coût en temps de calcul.

Un aspect vraiment agréable de cette approche est que, généralement, nous trouvons des valeurs extrêmes des choses dans le cadre d'un flux de travail plus vaste dans lequel un objectif doit être atteint: l'emplacement d'un pipeline ou d'un terrain de football, la création de tampons écologiques, etc. Le processus implique des compromis. Ainsi, la ligne horizontale la plus longue pourrait ne pas faire partie d'une solution optimale. On pourrait plutôt se soucier de savoir où presque les lignes les plus longues se trouveraient. C'est simple: au lieu de sélectionner le débit maximum zonal, sélectionnez toutes les cellules proches d'un maximum zonal. Dans cet exemple, le zonal max est égal à 744 (le nombre de colonnes enjambées par le segment intérieur le plus long). Au lieu de cela, sélectionnons toutes les cellules à moins de 5% de ce maximum:

L'exécution du flux d'est en ouest produit cette collection de segments horizontaux:

Il s'agit d'une carte des emplacements où l'étendue ininterrompue est-ouest est de 95% ou plus que l'étendue maximale est-ouest n'importe où dans le polygone.

D'accord. J'ai une (meilleure) idée ( idée-№2 ). Mais je suppose qu'il vaut mieux être réalisé comme un script python, pas comme un SQL-querry. Encore une fois, voici le cas commun, pas seulement EW.

Vous aurez besoin d'un cadre de délimitation pour le polygone et d'un azimut (A) comme direction de mesure. Supposons que la longueur des bords BBox soit LA et LB. La distance maximale possible (MD) dans un polygone est: MB = (LA^2 * LB^2)^(1/2), pour la recherche de valeur (V) ne soit pas plus grand que MB: V <= MB.

1. À partir de n'importe quel sommet de la BBox, créez une ligne (LL) avec la longueur MB et l'azimut A.
2. Intersection de la ligne LL avec le polygone pour obtenir la ligne d'intersection (IL)
3. Vérifiez la géométrie de l'IL - s'il n'y a que deux points dans la ligne IL, calculez sa longueur. Si 4 ou plus - calculez les segments et choisissez la longueur du plus long. Null (pas d'intersection du tout) - ignorer.
4. Continuez à créer d'autres lignes LL en vous déplaçant du compteur du point de départ ou dans le sens des aiguilles d'une montre vers les bords de la BBox jusqu'à ce que vous ne vous retrouviez pas au point de départ (vous ferez toute la boucle à travers la BBox).
5. Choisissez la plus grande valeur de longueur IL (en fait, vous n'avez pas besoin de stocker toutes les longueurs, vous pouvez conserver la valeur maximale `` jusqu'à présent '' pendant la boucle) - ce sera ce que vous recherchez.

Je ne suis pas sûr que la réponse de Fetzer soit ce que vous voulez faire, mais c'est ainsi que la st_box2d peut faire le travail.

L'idée N ° 1 de SS_Rebelious fonctionnera dans de nombreux cas mais pas pour certains polygones concaves.

Je pense que vous devez créer des lignes lw artificielles dont les points suivent les bords lorsque des lignes faites par des sommets traversent les frontières du polygone s'il existe une possibilité de ligne est-ouest.

Pour cela, vous pouvez essayer de créer un polygone à 4 nœuds où la longueur de la ligne est élevée, créer le polygone P * qui est le précédent chevauchant avec vous le polygone d'origine et voir si les min (y1) et max (y2) laissent une ligne x possibilité. (où y1 est l'ensemble de points entre le cornet supérieur gauche et le coin supérieur droit et y2 l'ensemble de y entre les coins inférieur gauche et inférieur droit de votre polygone à 4 nœuds). Ce n'est pas si simple j'espère que vous trouverez des outils psql pour vous aider!

J'ai une idée-№1 ( Edit: pour le cas commun, pas seulement la direction EW, et avec quelques limitations qui sont décrites dans les commentaires). Je ne fournirai pas le code, juste un concept. La "direction x" est en fait un azimut, qui est calculé par ST_Azimuth. Les étapes proposées sont les suivantes:

1. Extrayez tous les sommets du polygone en tant que points.
2. Créez des lignes entre chaque paire de points.
3. Sélectionnez des lignes (appelons-les lw-lines) qui se trouvent dans le polygone d'origine (nous n'avons pas besoin de lignes qui traverseront les frontières du polygone).
4. Trouvez les distances et les azimuts pour chaque ligne lw.
5. Sélectionnez la plus longue distance des lignes lw où l'azimut est égal à l'azimut recherché ou se situe dans un certain intervalle (il se peut qu'aucun azimut ne soit exactement égal à l'azimut recherché).

Jetez un œil à ma question et à la réponse de Evil Genius.

Si tout va bien votre polygone de lac a un certain nombre de points, vous pouvez créer des lignes sur ces points avec un azimut (aspect, direction géographique). Choisissez la longueur des lignes suffisamment grande (partie ST_MakePoint), afin de pouvoir calculer la ligne la plus courte entre les deux lignes les plus éloignées.

Voici un exemple:

L'exemple montre la largeur maximale du polygone. Je choisis ST_ShortestLine (ligne rouge) pour cette approche. ST_MakeLine augmenterait la valeur (ligne bleue) et l'extrémité de la ligne (en bas à gauche) atteindrait la ligne bleue du polygone. Vous devez calculer la distance avec les centroïdes des lignes (d'aide) créées.

Une idée de polygones irréguliers ou concaves pour cette approche. Il se peut que vous deviez couper le polygone avec un raster.