Déterminer si le point est entouré à l'aide du traitement raster

J'essaie d'améliorer un processus vectoriel / python extrêmement lourd pour un modèle de risque naturel. À l'heure actuelle, nous avons un long script qui génère des lignes de distance / relèvement à partir d'un point donné pour déterminer:

1. le type de polygone qu'il recoupe (p. ex. forêt, herbe, marais, etc.)
2. la distance à ce polygone
3. combien de ces lignes croisent des polygones, pour déterminer à quel point il est «entouré».

Il y en a beaucoup plus, mais c'est l'essentiel. J'essaie de trouver un moyen d'améliorer cela et je suis actuellement perplexe sur la partie 3. L'idée est de déterminer si un point est complètement entouré de polygones, disons à 200 m

Donc, dans mon image jointe, je voudrais que le point A soit marqué comme étant à plus haut risque que le point B car il est complètement entouré par mes polygones. Cela se répète pour environ 13 millions de points, ce n'est donc pas une mince tâche et je préférerais avoir une surface pour dériver des valeurs plutôt que d'exécuter notre script. Je pense qu'il doit y avoir une variation des outils d'hydrologie ou de la trajectoire des coûts pour ce faire, mais je n'arrive pas à comprendre.

Comment pourrais-je m'y prendre?

Il existe une solution de chemin de coûts mais vous devrez la coder vous-même. Voici à quoi cela pourrait ressembler lorsqu'il est appliqué à chaque point de l'image dans la question (grossi un peu pour accélérer les calculs):

Les cellules noires font partie des polygones environnants. Les couleurs, allant de l'orange clair (court) au bleu (long), indiquent la distance maximale (jusqu'à un maximum de 50 cellules) qui peut être atteinte par une traversée en ligne de visée sans intercepter les cellules du polygone. (Toute cellule en dehors de l'étendue de cette image est traitée comme faisant partie des polygones.)

Voyons un moyen efficace de le faire en utilisant une représentation raster des données. Dans cette représentation, toutes les cellules polygonales "environnantes" auront, disons, des valeurs non nulles et toute cellule qui peut être "vue à travers" aura une valeur nulle.

Étape 1: précalcul d'une structure de données de voisinage

Vous devez d'abord décider ce que signifie pour une cellule d'en bloquer une autre. L'une des règles les plus justes que je puisse trouver est la suivante: en utilisant des coordonnées intégrales pour les lignes et les colonnes (et en supposant des cellules carrées), considérons quelles cellules pourraient bloquer la cellule (i, j) de la vue à l'origine (0,0). Je nomme la cellule (i ', j') qui est la plus proche du segment de ligne reliant (i, j) à (0,0) parmi toutes les cellules dont les coordonnées diffèrent de i et j d'au plus 1. Parce que cela ne signifie pas toujours produire une solution unique (par exemple, avec (i, j) = (1,2) les deux (0,1) et (1,1) fonctionneront également bien), un moyen de résoudre les liens est nécessaire. Il serait bien que cette résolution de liens respecte les symétries des voisinages circulaires dans les grilles: annuler les coordonnées ou changer les coordonnées préserve ces voisinages. Par conséquent, nous pouvons décider quelles cellules bloquent (i,

Cette règle est illustrée par le code prototype suivant écrit R. Ce code renvoie une structure de données qui sera pratique pour déterminer le "entourage" de cellules arbitraires dans une grille.

screen <- function(k=1) {
  #
  # Returns a data structure:
  #   $offset is an array of offsets
  #   $screened is a parallel array of screened offset indexes.
  #   $distance is a parallel array of distances.
  # The first index always corresponds to (0,0).
  #
  screened.by <- function(xy) {
    uv <- abs(xy)
    if (reversed <- uv[2] > uv[1]) {
      uv <- rev(uv)
    }
    i <- which.min(c(uv[1], abs(uv[1]-uv[2]), uv[2]))
    ij <- uv + c(floor((1-i)/3), floor(i/3)-1)
    if (reversed) ij <- rev(ij)
    return(ij * sign(xy))
  }
  #
  # For each lattice point within the circular neighborhood,
  # find the unique lattice point that screens it from the origin.
  #
  xy <- subset(expand.grid(x=(-k:k), y=(-k:k)), 
               subset=(x^2+y^2 <= k^2) & (x != 0 | y != 0))
  g <- t(apply(xy, 1, function(z) c(screened.by(z), z)))
  #
  # Sort by distance from the origin.
  #
  colnames(g) <- c("x", "y", "x.to", "y.to")
  ij <- unique(rbind(g[, 1:2], g[, 3:4]))
  i <- order(abs(ij[,1]), abs(ij[,2])); ij <- ij[i, , drop=FALSE]
  rownames(ij) <- 1:length(i)
  #
  # Invert the "screened by" relation to produce the "screened" relation.
  #
  # (Row, column) offsets.
  ij.df <- data.frame(ij, i=1:length(i))
  #
  # Distances from the origin (in cells).
  distance <- apply(ij, 1, function(u) sqrt(sum(u*u)))
  #
  # "Screens" relation (represented by indexes into ij).
  g <- merge(merge(g, ij.df), ij.df, 
             by.x=c("x.to", "y.to"), by.y=c("x","y"))
  g <- subset(g, select=c(i.x, i.y))
  h <- by(g$i.y, g$i.x, identity)

  return( list(offset=ij, screened=h, distance=distance) )
}

La valeur de a screen(12)été utilisée pour produire cette représentation de cette relation de criblage: les flèches pointent des cellules vers celles qui les filtrent immédiatement. Les teintes sont proportionnées par la distance à l'origine, qui est au milieu de ce quartier:

Ce calcul est rapide et ne doit être effectué qu'une seule fois pour un quartier donné. Par exemple, lorsque vous regardez 200 m sur une grille avec 5 m de cellules, la taille du quartier sera de 200/5 = 40 unités.

Étape 2: Application du calcul aux points sélectionnés

Le reste est simple: pour déterminer si une cellule située en (x, y) (en coordonnées de ligne et de colonne) est "entourée" par rapport à cette structure de données de voisinage, effectuez le test récursivement en commençant par un décalage de (i, j) = (0,0) (l'origine du voisinage). Si la valeur dans la grille polygonale en (x, y) + (i, j) est différente de zéro, la visibilité y est bloquée. Sinon, nous devrons considérer tous les décalages qui auraient pu être bloqués au décalage (i, j) (qui se trouvent en temps O (1) en utilisant la structure de données retournée par screen). S'il n'y en a pas qui sont bloqués, nous avons atteint le périmètre et conclu que (x, y) n'est pas entouré, donc nous arrêtons le calcul (et ne nous soucions pas d'inspecter les points restants dans le voisinage).

Nous pouvons collecter des informations encore plus utiles en gardant une trace de la distance en ligne de visée la plus éloignée atteinte pendant l'algorithme. S'il est inférieur au rayon souhaité, la cellule est entourée; sinon ce n'est pas le cas.

Voici un Rprototype de cet algorithme. C'est plus long qu'il n'y paraît, car Ril ne prend pas en charge nativement la structure de pile (simple) nécessaire pour implémenter la récursivité, donc une pile doit également être codée. L'algorithme proprement dit commence aux deux tiers environ et n'a besoin que d'une douzaine de lignes environ. (Et la moitié de ceux-ci se contentent de gérer la situation autour du bord de la grille, vérifiant les indices hors plage dans le quartier. Cela pourrait être rendu plus efficace simplement en étendant la grille polygonale par des klignes et des colonnes autour de son périmètre, éliminant tout nécessité de vérifier la plage d'index au prix d'un peu plus de RAM pour contenir la grille polygonale.)

#
# Test a grid point `ij` for a line-of-sight connection to the perimeter
# of a circular neighborhood.  
#   `xy` is the grid.
#   `counting` determines whether to return max distance or count of stack ops.
#   `perimeter` is the assumed values beyond the extent of `xy`.
#
# Grid values of zero admit light; all others block visibility
# Returns maximum line-of-sight distance found within `nbr`.
#
panvisibility <- function(ij, xy, nbr=screen(), counting=FALSE, perimeter=1) {
  #
  # Implement a stack for the algorithm.
  #
  count <- 0 # Stack count
  stack <- list(ptr=0, s=rep(NA, dim(nbr$offset)[1]))
  push <- function(x) {
    n <- length(x)
    count <<- count+n         # For timing
    stack$s[1:n + stack$ptr] <<- x
    stack$ptr <<- stack$ptr+n
  }
  pop <- function() {
    count <<- count+1         # For timing
    if (stack$ptr <= 0) return(NULL)
    y <- stack$s[stack$ptr]
    #stack$s[stack$ptr] <<- NA # For debugging
    stack$ptr <<- stack$ptr - 1
    return(y)
  }
  #
  # Initialization.
  #
  m <- dim(xy)[1]; n <- dim(xy)[2]
  push(1) # Stack the *indexes* of nbr$offset and nbr$screened.
  dist.max <- -1
  #
  # The algorithm.
  #
  while (!is.null(i <- pop())) {
    cell <- nbr$offset[i, ] + ij
    if (cell[1] <= 0 || cell[1] > m || cell[2] <= 0 || cell[2] > n) {
      value <- perimeter
    } else {  
      value <- xy[cell[1], cell[2]]
    }
    if (value==0) {
      if (nbr$distance[i] > dist.max) dist.max <- nbr$distance[i]
      s <- nbr$screened[[paste(i)]]
      if (is.null(s)) {
        #exited = TRUE
        break
      }
      push(s)
    }
  }
  if (counting) return ( count )
  return(dist.max)
}

Dans cet exemple, les cellules polygonales sont noires. Les couleurs donnent la distance de visibilité maximale (jusqu'à 50 cellules) pour les cellules non polygonales, allant de l'orange clair pour les courtes distances au bleu foncé pour les plus longues distances. (Les cellules ont une unité de largeur et de hauteur.) Les stries visiblement évidentes sont créées par les petits "îlots" de polygones au milieu de la "rivière": chacun bloque une longue lignée d'autres cellules.

Analyse de l'algorithme

La structure de la pile implémente une recherche en profondeur du graphe de visibilité du voisinage pour trouver la preuve qu'une cellule n'est pas entourée. Lorsque les cellules sont éloignées de tout polygone, cette recherche nécessitera l'inspection des cellules O (k) uniquement pour un voisinage circulaire de rayon k. Les pires cas se produisent quand il y a un petit nombre de cellules polygonales dispersées dans le quartier, mais même ainsi la limite du quartier ne peut pas être tout à fait atteinte: celles-ci nécessitent d'inspecter presque toutes les cellules de chaque quartier, qui est un O (k ^ 2) opération.

Le comportement suivant est typique de ce qui sera rencontré. Pour les petites valeurs de k, à moins que les polygones ne remplissent la majeure partie de la grille, la plupart des cellules non polygonales ne seront évidemment pas entourées et l'algorithme évolue comme O (k). Pour les valeurs intermédiaires, la mise à l'échelle commence à ressembler à O (k ^ 2). Comme k devient vraiment grand, la plupart des cellules seront entourées et ce fait peut être déterminé bien avant que le voisinage entier ne soit inspecté: l'effort de calcul de l'algorithme atteint ainsi une limite pratique. Cette limite est atteinte lorsque le rayon de voisinage s'approche du diamètre des plus grandes régions non polygonales connectées de la grille.

À titre d'exemple, j'utilise l' countingoption codée dans le prototype de screenpour renvoyer le nombre d'opérations de pile utilisées dans chaque appel. Cela mesure l'effort de calcul. Le graphique suivant trace le nombre moyen d'opérations de pile en fonction du rayon de voisinage. Il présente le comportement prévu.

Nous pouvons l'utiliser pour estimer le calcul nécessaire pour évaluer 13 millions de points sur une grille. Supposons qu'un voisinage de k = 200/5 = 40 soit utilisé. Ensuite, quelques centaines d'opérations de pile seront nécessaires en moyenne (en fonction de la complexité de la grille polygonale et de l'emplacement des 13 millions de points par rapport aux polygones), ce qui implique que dans un langage compilé efficace, au plus quelques milliers d'opérations numériques simples sera nécessaire (addition, multiplication, lecture, écriture, décalage, etc.). La plupart des PC seront en mesure d'évaluer l'entourage d'environ un million de points à ce rythme. (LeRl'implémentation est beaucoup, beaucoup plus lente que cela, car elle est médiocre à ce type d'algorithme, c'est pourquoi il ne peut être considéré que comme un prototype.) Par conséquent, nous pouvons espérer qu'une implémentation efficace dans un langage raisonnablement efficace et approprié - C ++ et Python me vient à l'esprit - pourrait compléter l'évaluation de 13 millions de points en une minute ou moins, en supposant que la grille polygonale entière réside dans la RAM.

Lorsqu'une grille est trop grande pour tenir dans la RAM, cette procédure peut être appliquée aux parties en mosaïque de la grille. Ils n'ont qu'à se chevaucher par des klignes et des colonnes; prendre les maxima aux chevauchements lors de la mosaïque des résultats.

Autres applications

Le "fetch" d'un plan d'eau est étroitement lié à "l'entourage" de ses points. En fait, si nous utilisons un rayon de voisinage égal ou supérieur au diamètre du plan d'eau, nous créerons une grille de l'extraction (non directionnelle) à chaque point du plan d'eau. En utilisant un rayon de voisinage plus petit, nous obtiendrons au moins une limite inférieure pour l'extraction à tous les points d'extraction les plus élevés, ce qui, dans certaines applications, peut être suffisamment bon (et peut considérablement réduire l'effort de calcul). Une variante de cet algorithme qui limite la relation "filtré par" à des directions spécifiques serait une façon de calculer efficacement la récupération dans ces directions. Notez que ces variantes nécessitent de modifier le code de screen; le code pour panvisibilityne change pas du tout.

Je peux certainement voir comment on pourrait faire cela avec une solution raster, mais étant donné même un nombre réduit de points, je m'attendrais à une très grande / haute résolution et donc difficile à traiter une grille ou un ensemble de grilles. Compte tenu de cela, je me demande si l'exploitation de la topologie dans un gdb pourrait être plus efficace. Vous pouvez trouver tous les vides internes avec quelque chose comme:

arcpy.env.workspace = 'myGDB'
arcpy.CreateTopology_management('myGDB', 'myTopology', '')    
arcpy.AddFeatureClassToTopology_management('myTopology', 'myFeatures', '1','1')    
arcpy.AddRuleToTopology_management ('myToplogy', 'Must Not Have Gaps (Area)', 'myFeatures', '', '', '')    
arcpy.ValidateTopology_management('myTopology', 'Full_Extent')
arcpy.ExportTopologyErrors_management('myTopology', 'myGDB', 'topoErrors')
arcpy.FeatureToPolygon_management('topoErrors_line','topoErrorsVoidPolys', '0.1')`

vous pouvez ensuite travailler avec topoErrorsVoidPolysvotre modèle normal Intersect_analysis()ou autre. Vous devrez peut-être jouer avec l'extraction des polys d'intérêt topoErrorsVoidPolys. @whuber a un certain nombre d'excellents articles sur ce genre de choses ailleurs ici sur gis.stackexchange.com.

Si vous voulez vraiment aller raster ... je ferais quelque chose dans le sens de ce pseudo code (ne grincez pas juste parce qu'il est évident que je suis un retour AML!: P)

1. pixelliser les points ("pts_g") et les polys ("polys_g" (
2. voids = regiongroup (con (isnull (polys_g), 1))
3. pourrait avoir besoin de faire quelque chose pour affiner les vides afin d'éliminer les polygones externes indésirables / la zone de l'univers ouvert
4. pts_surrounded = con (voids, pts_g)

Juste pour inventer ça, il pourrait donc être nécessaire de l'affiner.

Si vous utilisez une valeur de distance seuil (ici vous parlez de 200 m), alors la meilleure solution est d'utiliser l'analyse vectorielle:

1) créez un tampon de 200 m autour de chaque point (en noir sur l'illustration)

2) utilisez l'outil d'intersection (analyse) entre le tampon et les polygones (en bleu sur l'illustration). Cela sera plus joli si vous faites cela entre les limites des polygones environnants et le tampon, mais le résultat final est le même.

3) utiliser la fonction de polygone (gestion) pour créer des polygones où vos points sont complètement entourés (en rouge sur l'illustration)

4) sélectionner des couches par emplacement (gestion) ou jointure spatiale (analyse) pour identifier les points qui sont entourés. L'utilisation de la jointure spatiale vous permet d'avoir une information sur le polygone d'intégration (aire du polygone, statistiques zonales ...) qui pourrait être utile pour un traitement ultérieur.

Alternatives 2b) En fonction de vos besoins, vous pouvez sélectionner par emplacement les polygones environnants à une distance de 200 m, vous pouvez alors identifier certains types d '«enclos» mais pas aussi strictement qu'en 2).

Compte tenu du «cas du labyrinthe», cela pourrait aider: évaluer combien de temps il faut pour «s'échapper» de l'emplacement.

• Vous pouvez déjà exclure de l'analyse les points qui sont entièrement inclus ou entièrement gratuits

• puis vous convertissez vos obstacles en raster et définissez les valeurs sur NoData où vous avez un polygone, et sur la taille de cellule en mètres où vous n'en avez pas (cela rendra votre coût raster).

• troisièmement, vous pouvez calculer la distance de coût en utilisant le raster de coût nouvellement généré

• enfin, vous utilisez une statistique zonale comme table basée sur les limites du tampon converties en raster (formant un espace annulaire). Si vous pouvez vous échapper dans toutes les directions, le minimum devrait être d'environ 200 (selon la taille des cellules de votre analyse). Mais si vous êtes dans un labyrinthe, le maximum sera supérieur à 200. Ainsi, le maximum des statistiques zonales moins 200 sera une valeur continue indiquant à quel point il est difficile de "s'échapper".