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

13

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.

cardiopteris
la source
3
Je ne connais pas la fonctionnalité postgis. Cependant, il peut y avoir un outil de boîte englobante. La largeur du cadre de délimitation serait la distance maximale dans la direction EW.
Fezter
4
@Fetzter qui n'est pas correct: un contre-exemple, même pour un simple polygone complexe, est un losange mince s'étendant de SW à NE. Sa largeur maximale est-ouest peut être une fraction arbitrairement petite de la largeur de sa boîte englobante.
whuber
1
J'ai créé un utilitaire pour cette tâche basé sur ceci et ces propositions. Il est capable de calculer la largeur max ou min du polygone. Actuellement, cela fonctionne avec un fichier shp, mais vous pouvez le réécrire pour travailler avec un PostGIS ou simplement attendre un certain temps jusqu'à ce qu'il évolue vers un plugin QGIS qui fonctionnera également avec PostGIS. Description détaillée et lien de téléchargement ici .
SS_Rebelious

Réponses:

16

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

Un polygone

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.

Après l'analyse

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)
whuber
la source
Existe-t-il un théorème qui prouve que la ligne de longueur maximale à l'intérieur d'un polygone non convexe dans une direction donnée coupe au moins un sommet de ce polygone?
SS_Rebelious
@SS Oui, il y en a. Voici une esquisse d'une preuve: s'il n'y a pas d'intersection, les extrémités du segment se trouvent toutes les deux à l'intérieur des bords et le segment peut être déplacé, au moins un peu, de haut en bas. Sa longueur est une fonction linéaire de la quantité de déplacement. Ainsi, il ne peut avoir une longueur maximale que si la longueur ne change pas lors du déplacement. Cela implique à la fois (a) qu'il fait partie d'un parallélogramme formé de segments de longueur maximale et (b) que les bords supérieur et inférieur de ce parallélogramme doivent rencontrer un sommet, QED.
whuber
Et quel est le nom de ce théorème? J'ai du mal à le trouver. BTW, qu'en est-il des bords courbes qui n'ont pas de sommet (je veux dire une approche théorique)? Un croquis d'exemple de la figure que je veux dire (un polygone en forme de haltère): "C = D".
SS_Rebelious
@SS Lorsque les arêtes sont courbes, le théorème ne tient plus. Des techniques de géométrie différentielle peuvent être appliquées pour obtenir des résultats utiles. J'ai appris ces méthodes dans le livre de Cheeger & Ebin, Comparison Theorems in Riemannian Geometry . Cependant, la plupart des SIG rapprocheront les courbes par des polylignes détaillées, de sorte que la question (en pratique) est théorique.
whuber
pourriez-vous spécifier le nom du théorème (et la page si possible)? J'ai le livre et je n'ai pas pu localiser le théorème nécessaire.
SS_Rebelious
9

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:

Polygone

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"):

Accumulation de flux

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:

Maximum

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!

Résultat

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:

Cellules quasi optimales sélectionnées

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

Solutions quasi optimales

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.

Whuber
la source
3

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.
SS_Rebelious
la source
Cela ressemble à une double boucle sur les sommets: c'est suffisamment inefficace pour qu'il soit évité (sauf pour les polygones très simplifiés).
whuber
@whuber, je ne vois pas de boucles supplémentaires ici. Il y a juste un traitement insignifiant de 2 côtés de BB qui ne donnera que des valeurs nulles. Mais ce traitement a été exclu dans le script que j'ai fourni dans la réponse qui a été supprimée ici (il semble que ce soit un commentaire maintenant, mais je ne le vois pas comme un commentaire - seulement comme une réponse supprimée)
SS_Rebelious
(1) Il s'agit du troisième commentaire à la question. (2) Vous avez raison: en lisant très attentivement votre description, il me semble maintenant que vous trouvez le segment le plus long entre (les quatre) sommets du cadre de délimitation et les sommets du polygone. Je ne vois pas comment cela répond à la question, cependant: le résultat n'est certainement pas ce que le PO recherchait.
whuber
@whuber, l'algorithme proposé trouve l'intersection la plus longue du polygone avec la ligne représentant la direction donnée. Apparemment, le résultat EST ce qui a été demandé si la distance entre les lignes d'intersection -> 0 ou si elle passe tous les sommets (pour les figures non courbes).
SS_Rebelious
3

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. un exemple où cela ne fonctionnera pas

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!

Un nom
la source
C'est sur la bonne voie. Le segment EW le plus long se trouvera parmi les intersections avec l'intérieur du polygone avec des lignes horizontales passant par les sommets du polygone. Cela nécessite que le code boucle sur les sommets. Il existe une méthode alternative (mais équivalente) disponible en suivant un flux artificiel est-ouest sur une représentation raster du polygone: la longueur de flux maximale trouvée dans le polygone (qui est l'une de ses "statistiques zonales") est la largeur souhaitée. La solution raster est obtenue en seulement 3 ou 4 étapes et ne nécessite ni boucle ni script.
whuber
@Aname, veuillez ajouter "№1" à "l'idée de SS_Rebelious" pour éviter tout malentendu: j'ai ajouté une autre proposition. Je ne peux pas modifier moi-même votre réponse car cette modification comporte moins de 6 caractères.
SS_Rebelious
1

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é).
SS_Rebelious
la source
Cela ne fonctionnera même pas pour certains triangles , comme celui aux sommets (0,0), (1000, 1000) et (501, 499). Sa largeur maximale est-ouest est d'environ 2; les azimuts sont tous autour de 45 degrés; et peu importe, le segment de ligne le plus court entre les sommets est plus de 350 fois plus long que la largeur est-ouest.
whuber
@whuber, vous avez raison, cela échouera pour les triangles, mais pour les polygones, représentant certaines caractéristiques de la nature, cela peut être utile.
SS_Rebelious
1
Il est difficile de recommander une procédure qui échoue de façon spectaculaire, même pour les cas simples, dans l'espoir qu'elle obtienne parfois une réponse correcte!
whuber
@whuber, alors ne le recommande pas! ;-) J'ai proposé cette solution de contournement car il n'y avait pas de réponse à cette question. Notez que vous pouvez publier votre propre meilleure réponse. BTW, si vous placez quelques points sur les bords du triangle, ma proposition fonctionnera ;-)
SS_Rebelious
J'ai proposé plusieurs approches. Une trame est disponible sur gis.stackexchange.com/questions/32552/… et élaborée sur forums.esri.com/Thread.asp?c=93&f=982&t=107703&mc=3 . Un autre - pas tout à fait comme applicable, mais avec des utilisations intéressantes - est à gis.stackexchange.com/questions/23664/… (la transformation du radon). Ceci est illustré sur stats.stackexchange.com/a/33102 .
whuber
1

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:

entrez la description de l'image ici

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.

Stefan
la source