Comment créer un tampon orienté en utilisant arcpy?

Sommaire

Cette réponse place la question dans un contexte plus large, décrit un algorithme efficace applicable à la représentation de fichiers de formes d'entités (comme «vecteurs» ou «chaînes de lignes» de points), montre quelques exemples de son application et donne un code de travail pour l'utilisation ou le portage dans un environnement SIG.

Contexte

Ceci est un exemple de dilatation morphologique. En toute généralité, une dilatation "étale" les points d'une région dans leurs quartiers; la collection de points où ils se terminent est la «dilatation». Les applications dans les SIG sont nombreuses: modélisation de la propagation du feu, du mouvement des civilisations, de la propagation des plantes, et bien plus encore.

Mathématiquement, et dans une très grande généralité (mais utile), une dilatation répartit un ensemble de points dans une variété riemannienne (comme un plan, une sphère ou un ellipsoïde). L'étalement est stipulé par un sous-ensemble du faisceau tangent en ces points. Cela signifie qu'à chacun de ces points un ensemble de vecteurs (directions et distances) est donné (j'appelle cela un "voisinage"); chacun de ces vecteurs décrit un chemin géodésique commençant à son point de base. Le point de base est "étalé" aux extrémités de tous ces chemins. (Pour la définition beaucoup plus limitée de «dilatation» qui est classiquement utilisée dans le traitement d'image, voir l'article Wikipedia . La fonction d'étalement est connue sous le nom de carte exponentielle en géométrie différentielle.)

La "mise en mémoire tampon" d'une entité est l'un des exemples les plus simples d'une telle dilatation: un disque de rayon constant (le rayon de la mémoire tampon) est créé (au moins conceptuellement) autour de chaque point de l'entité. L'union de ces disques est le tampon.

Cette question demande le calcul d'une dilatation légèrement plus compliquée où l'étalement ne peut se produire que dans une plage d'angles donnée (c'est-à-dire dans un secteur circulaire). Cela n'a de sens que pour les entités qui ne renferment aucune surface sensiblement incurvée (telles que les petites entités sur la sphère ou l'ellipsoïde ou toutes les entités dans le plan). Lorsque nous travaillons dans le plan, il est également important d'orienter tous les secteurs dans la même direction. (Cependant, si nous modélisons la propagation du feu par le vent, nous souhaiterions que les secteurs soient orientés avec le vent et que leurs tailles puissent également varier avec la vitesse du vent: c'est une des motivations de la définition générale de la dilatation que j'ai donnée. ) (Sur des surfaces courbes comme un ellipsoïde, il est en général impossible d'orienter tous les secteurs dans la "même" direction.)

Dans les circonstances suivantes, la dilatation est relativement facile à calculer:

L'entité est dans l'avion (c'est-à-dire que nous dilatons une carte de l'entité et, espérons-le, la carte est assez précise).
La dilatation sera constante : l'étalement en chaque point de la caractéristique se produira dans des voisinages congruents de même orientation.
Ce quartier commun est convexe. La convexité simplifie et accélère considérablement les calculs.

Cette question s'inscrit dans de telles circonstances spécialisées: elle demande la dilatation de polygones arbitraires par des secteurs circulaires dont les origines (les centres des disques dont ils sont issus) se situent aux points de base. À condition que ces secteurs ne s'étendent pas sur plus de 180 degrés, ils seront convexes. (Les secteurs plus grands peuvent toujours être divisés en deux en deux secteurs convexes; l'union des deux dilatations plus petites donnera le résultat souhaité.)

la mise en oeuvre

Parce que nous effectuons des calculs euclidiens - faisant l'étalement dans le plan - nous pouvons dilater un point simplement en traduisant le voisinage de dilatation en ce point. (Pour pouvoir le faire, le quartier a besoin d'une originequi correspondra au point de base. Par exemple, l'origine des secteurs dans cette question est le centre du cercle à partir duquel ils sont formés. Cette origine se trouve à la limite du secteur. Dans l'opération de tamponnage SIG standard, le voisinage est un cercle dont l'origine est en son centre; maintenant l'origine se trouve à l'intérieur du cercle. Le choix d'une origine n'est pas un gros problème de calcul, car un changement d'origine ne fait que déplacer la dilatation entière, mais cela peut être un gros problème en termes de modélisation des phénomènes naturels. La sectorfonction dans le code ci-dessous illustre comment une origine peut être spécifiée.)

Dilater un segment de ligne peut être délicat, mais pour un voisinage convexe, nous pouvons créer la dilatation comme l'union des dilatations des deux extrémités avec un parallélogramme soigneusement choisi. (Dans l'intérêt de l'espace, je ne m'arrêterai pas pour prouver des affirmations mathématiques comme celle-ci, mais j'encourage les lecteurs à essayer leurs propres preuves car c'est un exercice perspicace.) Voici une illustration utilisant trois secteurs (montrés en rose). Ils ont des rayons unitaires et leurs angles sont donnés dans les titres. Le segment de ligne lui-même a une longueur de 2, est horizontal et est représenté en noir:

Dilatations de segment

Les parallélogrammes sont trouvés en localisant les points roses qui sont aussi loin que possible du segment dans la direction verticale uniquement . Cela donne deux points inférieurs et deux points supérieurs le long de lignes parallèles au segment. Il suffit de réunir les quatre points en un parallélogramme (représenté en bleu). Remarquez, à droite, comment cela a du sens même lorsque le secteur lui-même n'est qu'un segment de ligne (et non un vrai polygone): là, chaque point du segment a été traduit dans une direction à 171 degrés à l'est du nord sur une distance allant de 0 à 1. L'ensemble de ces points d'extrémité est le parallélogramme illustré. Les détails de ce calcul apparaissent dans la bufferfonction définie dans dilate.edgesle code ci-dessous.

Pour dilater une polyligne , nous formons l'union des dilatations des points et segments qui la forment. Les deux dernières lignes dilate.edgeseffectuent cette boucle.

La dilatation d'un polygone nécessite d'inclure l'intérieur du polygone ainsi que la dilatation de sa frontière. (Cette assertion fait quelques hypothèses sur le quartier de dilatation. L'une est que tous les quartiers contiennent le point (0,0), ce qui garantit que le polygone est inclus dans sa dilatation. Dans le cas des quartiers variables, il suppose également que la dilatation de tout intérieur le point du polygone ne s'étendra pas au-delà de la dilatation des points limites. C'est le cas pour les voisinages constants.)

Examinons quelques exemples de la façon dont cela fonctionne, d'abord avec un nonagone (choisi pour révéler les détails) puis avec un cercle (choisi pour correspondre à l'illustration de la question). Les exemples continueront d'utiliser les trois mêmes quartiers, mais seront réduits à un rayon de 1/3.

Dilatations d'un nonagone

Sur cette figure, l'intérieur du polygone est gris, les dilatations ponctuelles (secteurs) sont roses et les dilatations de bord (parallélogrammes) sont bleues.

Dilatations d'un cercle

Le "cercle" n'est vraiment qu'un 60-gon, mais il se rapproche bien d'un cercle.

Performance

Lorsque l'entité de base est représentée par N points et le voisinage de dilatation par M points, cet algorithme nécessite un effort O (N M) . Cela doit être suivi en simplifiant le désordre des sommets et des arêtes dans l'union, ce qui peut nécessiter un effort O (N M log (N M)): c'est quelque chose à demander au SIG de faire; nous ne devrions pas avoir à programmer cela.

L'effort de calcul pourrait être amélioré à O (M + N) pour les entités de base convexes (car vous pouvez déterminer comment contourner la nouvelle frontière en fusionnant de manière appropriée les listes de sommets décrivant les limites des deux formes d'origine). Cela ne nécessiterait pas non plus de nettoyage ultérieur.

Lorsque le voisinage de dilatation change lentement de taille et / ou d'orientation au fur et à mesure que vous progressez autour de l'entité de base, la dilatation du bord peut être étroitement approchée à partir de la coque convexe de l'union des dilatations de ses extrémités. Si les deux quartiers de dilatation ont des points M1 et M2, cela peut être trouvé avec un effort O (M1 + M2) en utilisant un algorithme décrit dans Shamos & Preparata, Computational Geometry . Par conséquent, en laissant K = M1 + M2 + ... + M (N) le nombre total de sommets dans les N quartiers de dilatation, nous pouvons calculer la dilatation en temps O (K * log (K)).

Pourquoi voudrions-nous aborder une telle généralisation si tout ce que nous voulons, c'est un simple tampon? Pour les grandes entités terrestres, un voisinage de dilatation (tel qu'un disque) qui a en réalité une taille constante peut avoir une taille variable sur la carte où ces calculs sont effectués. Ainsi, nous obtenons un moyen de faire des calculs précis pour l'ellipsoïde tout en continuant à profiter de tous les avantages de la géométrie euclidienne.

Code

Les exemples ont été produits avec ce Rprototype, qui peut facilement être porté dans votre langage préféré (Python, C ++, etc.). Dans sa structure, elle est parallèle à l'analyse rapportée dans cette réponse et n'a donc pas besoin d'explication séparée. Les commentaires clarifient certains détails.

(Il peut être intéressant de noter que les calculs trigonométriques ne sont utilisés que pour créer les exemples d'entités - qui sont des polygones réguliers - et les secteurs. Aucune partie des calculs de dilatation ne nécessite de trigonométrie.)

#
# Dilate the vertices of a polygon/polyline by a shape.
#
dilate.points <- function(p, q) {
  # Translate a copy of `q` to each vertex of `p`, resulting in a list of polygons.
  pieces <- apply(p, 1, function(x) list(t(t(q)+x)))
  lapply(pieces, function(z) z[[1]]) # Convert to a list of matrices
}
#
# Dilate the edges of a polygon/polyline `p` by a shape `q`. 
# `p` must have at least two rows.
#
dilate.edges <- function(p, q) {
  i <- matrix(c(0,-1,1,0), 2, 2)       # 90 degree rotation
  e <- apply(rbind(p, p[1,]), 2, diff) # Direction vectors of the edges
  # Dilate a single edge from `x` to `x+v` into a parallelogram
  # bounded by parts of the dilation shape that are at extreme distances
  # from the edge.
  buffer <- function(x, v) {
    y <- q %*% i %*% v # Signed distances orthogonal to the edge
    k <- which.min(y)  # Find smallest distance, then the largest *after* it
    l <- (which.max(c(y[-(1:k)], y[1:k])) + k-1) %% length(y)[1] + 1
    list(rbind(x+q[k,], x+v+q[k,], x+v+q[l,], x+q[l,])) # A parallelogram
  }
  # Apply `buffer` to every edge.
  quads <- apply(cbind(p, e), 1, function(x) buffer(x[1:2], x[3:4]))
  lapply(quads, function(z) z[[1]]) # Convert to a list of matrices
}
#----------------------- (This ends the dilation code.) --------------------------#
#
# Display a polygon and its point and edge dilations.
# NB: In practice we would submit the polygon, its point dilations, and edge 
#     dilations to the GIS to create and simplify their union, producing a single
#     polygon.  We keep the three parts separate here in order to illustrate how
#     that polygon is constructed.
#
display <- function(p, d.points, d.edges, ...) {
  # Create a plotting region covering the extent of the dilated figure.
  x <- c(p[,1], unlist(lapply(c(d.points, d.edges), function(x) x[,1])))
  y <- c(p[,2], unlist(lapply(c(d.points, d.edges), function(x) x[,2])))
  plot(c(min(x),max(x)), c(min(y),max(y)), type="n", asp=1, xlab="x", ylab="y", ...)
  # The polygon itself.
  polygon(p, density=-1, col="#00000040")
  # The dilated points and edges.
  plot.list <- function(l, c) lapply(l, function(p) 
                  polygon(p, density=-1, col=c, border="#00000040"))
  plot.list(d.points, "#ff000020")
  plot.list(d.edges, "#0000ff20")
  invisible(NULL) # Doesn't return anything
}
#
# Create a sector of a circle.
# `n` is the number of vertices to use for approximating its outer arc.
#
sector <- function(radius, arg1, arg2, n=1, origin=c(0,0)) {
  t(cbind(origin, radius*sapply(seq(arg1, arg2, length.out=n), 
                  function(a) c(cos(a), sin(a)))))
}
#
# Create a polygon represented as an array of rows.
#
n.vertices <- 60 # Inscribes an `n.vertices`-gon in the unit circle.
angles <- seq(2*pi, 0, length.out=n.vertices+1)
angles <- angles[-(n.vertices+1)]
polygon.the <- cbind(cos(angles), sin(angles))
if (n.vertices==1) polygon.the <- rbind(polygon.the, polygon.the)
#
# Dilate the polygon in various ways to illustrate.
#
system.time({
  radius <- 1/3
  par(mfrow=c(1,3))
  q <- sector(radius, pi/12, 2*pi/3, n=120)
  d.points <- dilate.points(polygon.the, q)
  d.edges <- dilate.edges(polygon.the, q)
  display(polygon.the, d.points, d.edges, main="-30 to 75 degrees")

  q <- sector(radius, pi/3, 4*pi/3, n=180)
  d.points <- dilate.points(polygon.the, q)
  d.edges <- dilate.edges(polygon.the, q)
  display(polygon.the, d.points, d.edges, main="-150 to 30 degrees")

  q <- sector(radius, -9/20*pi, -9/20*pi)
  d.points <- dilate.points(polygon.the, q)
  d.edges <- dilate.edges(polygon.the, q)
  display(polygon.the, d.points, d.edges, main="171 degrees")
})

Le temps de calcul pour cet exemple (d'après la dernière figure), avec N = 60 et M = 121 (à gauche), M = 181 (au milieu) et M = 2 (à droite), était d'un quart de seconde. Cependant, la plupart de cela était pour l'affichage. En règle générale, ce Rcode gérera environ N M = 1,5 million par seconde (en ne prenant que 0,002 seconde environ pour effectuer tous les exemples de calcul illustrés). Néanmoins, l'apparition du produit M N implique des dilatations de nombreuses figures ou des figures compliquées via un voisinage détaillé pouvant prendre un temps considérable, alors attention! Benchmarkez le timing des petits problèmes avant d'en aborder un gros. Dans de telles circonstances, on pourrait envisager une solution basée sur un raster (qui est beaucoup plus facile à implémenter, ne nécessitant essentiellement qu'un seul calcul de voisinage.)

whuber
la source

Wow, c'est très détaillé et fascinant. Je n'en attendais pas moins.

Paul