¿Cómo crear un búfer orientado usando arcpy?

Me gustaría crear un búfer orientado para cada polígono en mi shapefile usando arcpy. Por orientado quiero decir que tengo dos ángulos a1 y a2 que limitan la dirección del búfer. Esto se representa en el siguiente gráfico:

¿Algunas ideas?

Resumen

Esta respuesta coloca la pregunta en un contexto más amplio, describe un algoritmo eficiente aplicable a la representación de archivos de formas de características (como "vectores" o "cadenas lineales" de puntos), muestra algunos ejemplos de su aplicación y proporciona un código de trabajo para usar o portar Un entorno SIG.

Antecedentes

Este es un ejemplo de dilatación morfológica. En general, una dilatación "extiende" los puntos de una región a sus vecindarios; La colección de puntos donde terminan es la "dilatación". Las aplicaciones en SIG son numerosas: modelar la propagación del fuego, el movimiento de civilizaciones, la propagación de plantas y mucho más.

Matemáticamente, y en generalidad muy grande (pero útil), una dilatación extiende un conjunto de puntos en una variedad Riemanniana (como un plano, esfera o elipsoide). La extensión se estipula mediante un subconjunto del paquete tangente en estos puntos. Esto significa que en cada uno de los puntos se da un conjunto de vectores (direcciones y distancias) (lo llamo un "vecindario"); cada uno de estos vectores describe un camino geodésico que comienza en su punto base. El punto base se "extiende" a los extremos de todos estos caminos. (Para la definición mucho más limitada de "dilatación" que se emplea convencionalmente en el procesamiento de imágenes, consulte el artículo de Wikipedia . La función de difusión se conoce como el mapa exponencial en geometría diferencial)

El " almacenamiento en búfer " de una característica es uno de los ejemplos más simples de dicha dilatación: se crea un disco de radio constante (el radio del búfer) (al menos conceptualmente) alrededor de cada punto de la característica. La unión de estos discos es el búfer.

Esta pregunta solicita el cálculo de una dilatación un poco más complicada en la que se permite que la expansión ocurra solo dentro de un rango dado de ángulos (es decir, dentro de un sector circular). Esto tiene sentido solo para entidades que no encierran una superficie curvada apreciable (como entidades pequeñas en la esfera o elipsoide o cualquier entidad en el plano). Cuando trabajamos en el avión, también es significativo orientar todos los sectores en la misma dirección. (Sin embargo, si estuviéramos modelando la propagación del fuego por el viento, desearíamos que los sectores estén orientados con el viento y sus tamaños también podrían variar con la velocidad del viento: esa es una motivación para la definición general de dilatación que di. ) (En superficies curvas como un elipsoide, en general es imposible orientar todos los sectores en la "misma" dirección).

En las siguientes circunstancias, la dilatación es relativamente fácil de calcular:

• La característica está en el plano (es decir, estamos dilatando un mapa de la característica y esperamos que el mapa sea bastante preciso).

• La dilatación será constante : la propagación en cada punto de la característica se producirá en vecindarios congruentes de la misma orientación.

• Este barrio común es convexo. La convexidad simplifica enormemente y acelera los cálculos.

Esta pregunta se ajusta a circunstancias tan especializadas: pide la dilatación de polígonos arbitrarios por sectores circulares cuyos orígenes (los centros de los discos de los que provienen) se encuentran en los puntos base. Siempre que esos sectores no abarquen más de 180 grados, serán convexos. (Los sectores más grandes siempre se pueden dividir por la mitad en dos sectores convexos; la unión de las dos dilataciones más pequeñas dará el resultado deseado).

Implementación

Debido a que estamos realizando cálculos euclidianos, haciendo la expansión en el plano, podemos dilatar un punto simplemente traduciendo el vecindario de dilatación a ese punto. (Para poder hacer esto, el vecindario necesita un origeneso corresponderá al punto base. Por ejemplo, el origen de los sectores en esta pregunta es el centro del círculo a partir del cual se forman. Este origen pasa a estar en la frontera del sector. En la operación estándar de almacenamiento intermedio de SIG, el vecindario es un círculo con origen en su centro; ahora el origen se encuentra en el interior del círculo. Elegir un origen no es un gran problema computacional, porque un cambio de origen simplemente cambia toda la dilatación, pero puede ser un gran problema en términos de modelado de fenómenos naturales. La sectorfunción en el siguiente código ilustra cómo se puede especificar un origen).

Dilatar un segmento de línea puede ser complicado, pero para un vecindario convexo podemos crear la dilatación como la unión de dilataciones de los dos puntos finales junto con un paralelogramo cuidadosamente elegido. (En aras del espacio, no haré una pausa para probar afirmaciones matemáticas como esta, pero alentaré a los lectores a que intenten sus propias pruebas porque es un ejercicio perspicaz). Aquí hay una ilustración que usa tres sectores (que se muestran en rosa). Tienen unidades de radio y sus ángulos se dan en los títulos. El segmento de línea en sí tiene longitud 2, es horizontal y se muestra en negro:

Los paralelogramos se encuentran localizando los puntos rosados ​​que están lo más lejos posible del segmento solo en dirección vertical . Esto proporciona dos puntos inferiores y dos puntos superiores a lo largo de líneas paralelas al segmento. Solo tenemos que unir los cuatro puntos en un paralelogramo (que se muestra en azul). Observe, a la derecha, cómo tiene sentido incluso cuando el sector en sí mismo es solo un segmento de línea (y no un verdadero polígono): allí, cada punto del segmento se ha traducido en una dirección 171 grados al este del norte para una distancia que varía de 0 a 1. El conjunto de estos puntos finales es el paralelogramo que se muestra. Los detalles de este cálculo aparecen en la bufferfunción definida dilate.edgesen el siguiente código.

Para dilatar una polilínea , formamos la unión de las dilataciones de los puntos y segmentos que la forman. Las dos últimas líneas de dilate.edgesrealizan este bucle.

La dilatación de un polígono requiere incluir el interior del polígono junto con la dilatación de su límite. (Esta afirmación hace algunas suposiciones sobre el vecindario de dilatación. Una es que todos los vecindarios contienen el punto (0,0), lo que garantiza que el polígono se incluye dentro de su dilatación. En el caso de vecindarios variables también supone que la dilatación de cualquier interior el punto del polígono no se extenderá más allá de la dilatación de los puntos límite. Este es el caso de vecindades constantes).

Veamos algunos ejemplos de cómo funciona esto, primero con un nonágono (elegido para revelar detalles) y luego con un círculo (elegido para que coincida con la ilustración de la pregunta). Los ejemplos continuarán usando los mismos tres vecindarios, pero se redujeron a un radio de 1/3.

En esta figura, el interior del polígono es gris, las dilataciones de los puntos (sectores) son rosas y las dilataciones de los bordes (paralelogramos) son azules.

El "círculo" es realmente solo un 60 gon, pero se aproxima muy bien a un círculo.

Actuación

Cuando la característica base está representada por N puntos y la vecindad de dilatación por M puntos, este algoritmo requiere un esfuerzo O (N M) . Debe seguirse simplificando el desorden de vértices y bordes en la unión, lo que puede requerir un esfuerzo O (N M log (N M)): eso es algo que debe pedirle al SIG que haga; no deberíamos tener que programar eso.

El esfuerzo computacional podría mejorarse a O (M + N) para entidades base convexas (porque puede averiguar cómo viajar alrededor del nuevo límite fusionando adecuadamente las listas de vértices que describen los límites de las dos formas originales). Esto tampoco necesitaría ninguna limpieza posterior.

Cuando la vecindad de dilatación cambia lentamente de tamaño y / u orientación a medida que avanza alrededor de la entidad base, la dilatación del borde puede aproximarse estrechamente desde el casco convexo de la unión de las dilataciones de sus puntos finales. Si los dos vecindarios de dilatación tienen puntos M1 y M2, esto se puede encontrar con el esfuerzo O (M1 + M2) utilizando un algoritmo descrito en Shamos & Preparata, Computational Geometry . Por lo tanto, dejando que K = M1 + M2 + ... + M (N) sea el número total de vértices en los vecindarios de dilatación de N, podemos calcular la dilatación en el tiempo O (K * log (K)).

¿Por qué querríamos abordar tal generalización si todo lo que queremos es un búfer simple? Para grandes características en la tierra, un vecindario de dilatación (como un disco) que tiene un tamaño constante en la realidad puede tener un tamaño variable en el mapa donde se realizan estos cálculos. Por lo tanto, obtenemos una manera de hacer cálculos que son precisos para el elipsoide mientras continuamos disfrutando de todas las ventajas de la geometría euclidiana.

Código

Los ejemplos se produjeron con este Rprototipo, que se puede transportar fácilmente a su lenguaje favorito (Python, C ++, etc.). En estructura, es paralelo al análisis que se informa en esta respuesta y, por lo tanto, no necesita una explicación por separado. Los comentarios aclaran algunos de los detalles.

(Puede ser interesante observar que los cálculos trigonométricos se usan solo para crear las entidades de ejemplo, que son polígonos regulares, y los sectores. Ninguna parte de los cálculos de dilatación requiere trigonometría).

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

El tiempo de cálculo para este ejemplo (de la última figura), con N = 60 y M = 121 (izquierda), M = 181 (centro) y M = 2 (derecha), fue un cuarto de segundo. Sin embargo, la mayor parte de esto fue para la pantalla. Típicamente, este Rcódigo manejará aproximadamente N M = 1.5 millones por segundo (tomar solo 0.002 segundos más o menos para hacer todos los cálculos de ejemplo que se muestran). Sin embargo, la apariencia del producto M N implica dilataciones de muchas figuras o figuras complicadas a través de un vecindario detallado que puede llevar un tiempo considerable, ¡así que tenga cuidado! Compare el momento para problemas menores antes de abordar uno grande. En tales circunstancias, uno podría buscar una solución basada en ráster (que es mucho más fácil de implementar, que requiere esencialmente un solo cálculo de vecindad).

Esto es bastante amplio, pero podrías:

1. amortiguar el polígono original
2. encontrar el punto de origen de los rayos "orientados" que se crearán en el límite del polígono (¿algún tipo de punto de tangencia?)
3. cree / extienda una línea desde ese punto a una distancia más allá del búfer usando el ángulo discutido en los comentarios a la pregunta.
4. interseca esa línea con el búfer y el polígono original. Esto probablemente podría hacerse al mismo tiempo que 3) con argumentos adecuados para extender.
5. extraer el nuevo polígono "buffer orientado" del conjunto resultante de polígonos