Determinar si el punto está rodeado usando el procesamiento ráster

Estoy intentando mejorar un proceso de vector / pitón extremadamente engorroso para un modelo de peligro natural. Por el momento tenemos una secuencia de comandos larga que genera líneas de distancia / rumbo desde un punto determinado para determinar:

1. El tipo de polígono que se cruza (p. ej., bosque, hierba, pantano, etc.)
2. la distancia a ese polígono
3. cuántas de estas líneas se cruzan con los polígonos, para determinar qué tan 'rodeado' está.

Hay mucho más involucrado, pero esa es la esencia de esto. Estoy tratando de encontrar una manera de mejorar esto y actualmente estoy perplejo en la parte 3. La idea es determinar si un punto está completamente rodeado de polígonos, dentro de unos 200 m.

Entonces, en mi imagen adjunta, me gustaría que el punto A esté marcado como de mayor riesgo que el punto B, ya que está completamente rodeado por mis polígonos. Esto se repite durante unos 13 millones de puntos, por lo que no es una tarea pequeña y preferiría tener una superficie para obtener valores, en lugar de ejecutar nuestro script. Estoy pensando que debe haber una variación de las herramientas de hidrología o la ruta de costos para hacer esto, pero parece que no puedo entenderlo.

¿Cómo podría hacer esto?

Existe una solución de ruta de costo, pero tendrá que codificarla usted mismo. Esto es lo que podría parecer cuando se aplica a cada punto de la imagen en la pregunta (engrosado un poco para acelerar los cálculos):

Las celdas negras son partes de los polígonos circundantes. Los colores, que van desde el naranja claro (corto) hasta el azul (largo), muestran la distancia máxima (hasta un máximo de 50 celdas) que se puede alcanzar a través del recorrido de la línea de visión sin interceptar las celdas poligonales. (Cualquier celda fuera del alcance de esta imagen se trata como parte de los polígonos).

Analicemos una forma eficiente de hacerlo utilizando una representación ráster de los datos. En esta representación, todas las celdas poligonales "circundantes" tendrán, por ejemplo, valores distintos de cero y cualquier celda que se pueda "ver a través" tendrá un valor cero.

Paso 1: precomputar una estructura de datos de vecindario

Primero debe decidir qué significa que una celda bloquee a otra. Una de las reglas más justas que puedo encontrar es esta: usando coordenadas integrales para filas y columnas (y suponiendo celdas cuadradas), consideremos qué celdas podrían bloquear la celda (i, j) de la vista en el origen (0,0). Nomino la celda (i ', j') que está más cerca del segmento de línea que conecta (i, j) a (0,0) entre todas las celdas cuyas coordenadas difieren de i y j como máximo 1. Porque esto no siempre producir una solución única (por ejemplo, con (i, j) = (1,2) ambos (0,1) y (1,1) funcionarán igualmente bien), se necesitan algunos medios para resolver los lazos. Sería bueno para esta resolución de vínculos respetar las simetrías de vecindarios circulares en cuadrículas: negar las coordenadas o cambiar las coordenadas preserva estos vecindarios. Por lo tanto, podemos decidir qué células bloquean (i,

Ilustrando esta regla está el siguiente código prototipo escrito R. Este código devuelve una estructura de datos que será conveniente para determinar el "cerco" de celdas arbitrarias en una cuadrícula.

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

El valor de screen(12)se utilizó para producir esta representación de esta relación de detección: las flechas apuntan desde las celdas a las que las detectan de inmediato. Los tonos están proporcionados por la distancia al origen, que se encuentra en el medio de este vecindario:

Este cálculo es rápido y debe hacerse solo una vez para un vecindario determinado. Por ejemplo, al mirar 200 m en una cuadrícula con celdas de 5 m, el tamaño del vecindario será 200/5 = 40 unidades.

Paso 2: aplicar el cálculo a los puntos seleccionados

El resto es sencillo: para determinar si una celda ubicada en (x, y) (en coordenadas de fila y columna) está "rodeada" con respecto a esta estructura de datos de vecindad, realice la prueba de forma recursiva comenzando con un desplazamiento de (i, j) = (0,0) (el origen del barrio). Si el valor en la cuadrícula del polígono en (x, y) + (i, j) no es cero, entonces la visibilidad se bloquea allí. De lo contrario, tendremos que considerar todos los desplazamientos que podrían haberse bloqueado en el desplazamiento (i, j) (que se encuentran en el tiempo O (1) utilizando la estructura de datos devuelta por screen). Si no hay ninguno bloqueado, hemos alcanzado el perímetro y concluimos que (x, y) no está rodeado, por lo que detenemos el cálculo (y no nos molestamos en inspeccionar los puntos restantes en el vecindario).

Podemos recopilar información aún más útil al realizar un seguimiento de la distancia de la línea de visión más remota alcanzada durante el algoritmo. Si esto es menor que el radio deseado, la celda está rodeada; de lo contrario no lo es.

Aquí hay un Rprototipo de este algoritmo. Es más largo de lo que parece, porque Rno admite de forma nativa la estructura de pila (simple) necesaria para implementar la recursión, por lo que una pila también debe codificarse. El algoritmo real comienza aproximadamente dos tercios del proceso y solo necesita una docena de líneas más o menos. (Y la mitad de ellos simplemente maneja la situación alrededor del borde de la cuadrícula, buscando índices fuera de rango dentro del vecindario. Esto podría hacerse más eficiente simplemente expandiendo la cuadrícula de polígonos por kfilas y columnas alrededor de su perímetro, eliminando cualquier necesidad de verificar el rango del índice a costa de un poco más de RAM para mantener la cuadrícula poligonal).

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

En este ejemplo, las celdas poligonales son negras. Los colores proporcionan la distancia máxima de la línea de visión (hasta 50 celdas) para celdas no poligonales, que van desde el naranja claro para distancias cortas hasta el azul oscuro para las distancias más largas. (Las celdas son de una unidad de ancho y alto). Las rayas visiblemente evidentes son creadas por las pequeñas "islas" poligonales en el medio del "río": cada una bloquea una larga línea de otras celdas.

Análisis del algoritmo.

La estructura de la pila implementa una búsqueda en profundidad del gráfico de visibilidad del vecindario para buscar evidencia de que una celda no está rodeada. Cuando las celdas están lejos de cualquier polígono, esta búsqueda requerirá la inspección de solo celdas O (k) para una vecindad circular de radio-k. Los peores casos ocurren cuando hay un pequeño número de celdas poligonales dispersas dentro del vecindario, pero aun así no se puede alcanzar el límite del vecindario: esto requiere inspeccionar casi todas las celdas en cada vecindario, lo cual es un O (k ^ 2) operación.

El siguiente comportamiento es típico de lo que se encontrará. Para valores pequeños de k, a menos que los polígonos llenen la mayor parte de la cuadrícula, la mayoría de las celdas no poligonales obviamente no estarán rodeadas y el algoritmo se escalará como O (k). Para valores intermedios, la escala comienza a parecerse a O (k ^ 2). A medida que k se vuelve realmente grande, la mayoría de las celdas se rodearán y ese hecho se puede determinar mucho antes de que se inspeccione todo el vecindario: el esfuerzo computacional del algoritmo alcanza así un límite práctico. Este límite se alcanza cuando el radio de vecindad se acerca al diámetro de las regiones no poligonales conectadas más grandes en la cuadrícula.

Como ejemplo, uso la countingopción codificada en el prototipo de screenpara devolver el número de operaciones de pila utilizadas en cada llamada. Esto mide el esfuerzo computacional. El siguiente gráfico traza el número medio de operaciones de pila en función del radio de vecindad. Exhibe el comportamiento predicho.

Podemos usar esto para estimar el cálculo necesario para evaluar 13 millones de puntos en una cuadrícula. Suponga que se usa una vecindad de k = 200/5 = 40. Luego, se necesitarán unos cientos de operaciones de apilamiento en promedio (dependiendo de la complejidad de la cuadrícula del polígono y de la ubicación de los 13 millones de puntos en relación con los polígonos), lo que implica que en un lenguaje compilado eficiente, a lo sumo unos pocos miles de operaciones numéricas simples será requerido (sumar, multiplicar, leer, escribir, compensar, etc.). La mayoría de las PC podrán evaluar la cercanía de aproximadamente un millón de puntos a esa velocidad. (LosRla implementación es mucho, mucho más lenta que eso, porque es pobre en este tipo de algoritmo, por lo que solo puede considerarse un prototipo.) En consecuencia, podríamos esperar que una implementación eficiente en un lenguaje razonablemente eficiente y apropiado - C ++ y Python me viene a la mente: podría completar la evaluación de 13 millones de puntos en un minuto o menos, suponiendo que toda la cuadrícula de polígonos reside en la RAM.

Cuando una cuadrícula es demasiado grande para caber en la RAM, este procedimiento se puede aplicar a porciones de mosaico de la cuadrícula. Solo tienen que superponerse por kfilas y columnas; tome los máximos en las superposiciones al mosaizar los resultados.

Otras aplicaciones

La "búsqueda" de un cuerpo de agua está estrechamente relacionada con la "cercanía" de sus puntos. De hecho, si usamos un radio de vecindad igual o mayor que el diámetro del cuerpo de agua, crearemos una cuadrícula de la búsqueda (no direccional) en cada punto del cuerpo de agua. Al usar un radio de vecindad más pequeño, al menos obtendremos un límite inferior para la búsqueda en todos los puntos de búsqueda más altos, que en algunas aplicaciones pueden ser lo suficientemente buenos (y pueden reducir sustancialmente el esfuerzo computacional). Una variante de este algoritmo que limita la relación "apantallado por" a direcciones específicas sería una forma de calcular la recuperación de manera eficiente en esas direcciones. Tenga en cuenta que tales variantes requieren modificar el código para screen; el código para panvisibilityno cambia en absoluto.

Definitivamente puedo ver cómo uno podría querer hacer esto con una solución de trama, pero dado incluso un número reducido de puntos, esperaría una resolución muy grande / alta y, por lo tanto, difícil de procesar la cuadrícula o el conjunto de cuadrículas. Dado eso, me pregunto si explotar la topología en un gdb podría ser más eficiente. Podrías encontrar todos los vacíos internos con algo como:

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')`

entonces puedes trabajar con topoErrorsVoidPolystu patrón normal Intersect_analysis()o lo que sea. Es posible que deba perder el tiempo extrayendo los polys de interés topoErrorsVoidPolys. @whuber tiene una cantidad de publicaciones bastante excelentes sobre este tipo de cosas en otros lugares aquí en gis.stackexchange.com.

Si realmente quieres ir a la trama ... Haría algo similar a este pseudocódigo (¡no te avergüences solo porque es obvio que soy un retroceso de AML!: P)

1. rasterizar puntos ("pts_g") y polys ("polys_g" (
2. voids = regiongroup (con (isnull (polys_g), 1))
3. podría necesitar hacer algo para refinar huecos para eliminar polígonos externos no deseados / área de universo abierto
4. pts_suritated = con (vacíos, pts_g)

Simplemente inventando eso, por lo que podría necesitar refinamiento.

Si usa un valor de distancia umbral (aquí habla de unos 200 m), entonces la mejor solución es usar el análisis vectorial:

1) cree un búfer de 200 m alrededor de cada punto (en negro en la ilustración)

2) use la herramienta de intersección (análisis) entre el búfer y los polígonos (en azul en la ilustración). Se verá mejor si hace esto entre los límites de los polígonos circundantes y el búfer, pero el resultado final es el mismo.

3) use la función de polígono (gestión) para crear polígonos donde sus puntos estén completamente rodeados (en rojo en la ilustración)

4) seleccione capas por ubicación (gestión) o unión espacial (análisis) para identificar los puntos que están rodeados. El uso de la unión espacial le permite tener una información sobre el polígono de incrustación (área del polígono, estadísticas zonales ...) que podría ser útil para su posterior procesamiento.

Alternativas 2b) Dependiendo de sus necesidades, puede seleccionar por ubicación los polígonos circundantes dentro de una distancia de 200 m, luego puede identificar algunos tipos de "cerramiento" pero no tan estrictamente como en 2).

Teniendo en cuenta el "caso del laberinto", esto podría ayudar: evaluar cuánto tiempo es "escapar" de la ubicación.

• Ya puede excluir del análisis los puntos que están completamente incluidos o totalmente libres

• luego convierte sus obstáculos en un ráster y establece los valores en NoData donde tiene un polígono, y en el tamaño de celda en metro donde no lo tiene (esto hará que su costo sea ráster).

• tercero, puede calcular la distancia de costo utilizando el ráster de costo recién generado

• finalmente, utiliza una estadística zonal como tabla basada en los límites del búfer convertidos a ráster (formando un anillo). Si puede escapar en todas las direcciones, el mínimo debe ser de aproximadamente 200 (dependiendo del tamaño de celda de su análisis). Pero si está en un laberinto, el máximo será mayor que 200. Por lo tanto, el máximo de las estadísticas zonales menos 200 será un valor continuo que indica lo difícil que es "escapar".