Como criar um buffer orientado usando o arcpy?

Sumário

Essa resposta coloca a questão em um contexto maior, descreve um algoritmo eficiente aplicável à representação de recursos de shapefile (como "vetores" ou "cadeias de linhas" de pontos), mostra alguns exemplos de sua aplicação e fornece código de trabalho para usar ou portar um ambiente GIS.

fundo

Este é um exemplo de dilatação morfológica. Em geral, uma dilatação "espalha" os pontos de uma região em seus bairros; a coleção de pontos onde eles terminam é a "dilatação". As aplicações no SIG são numerosas: modelando a propagação do fogo, o movimento das civilizações, a disseminação das plantas e muito mais.

Matematicamente, e em uma generalidade muito grande (mas útil), uma dilatação espalha um conjunto de pontos em uma variedade Riemanniana (como um plano, esfera ou elipsóide). A expansão é estipulada por um subconjunto do feixe tangente nesses pontos. Isso significa que em cada um dos pontos é dado um conjunto de vetores (direções e distâncias) (eu chamo isso de "vizinhança"); cada um desses vetores descreve um caminho geodésico começando em seu ponto base. O ponto base é "espalhado" até os fins de todos esses caminhos. (Para a definição muito mais limitada de "dilatação" convencionalmente empregada no processamento de imagens, consulte o artigo da Wikipedia . A função de espalhamento é conhecida como mapa exponencial em geometria diferencial.)

O "buffer" de um recurso é um dos exemplos mais simples dessa dilatação: um disco de raio constante (o raio do buffer) é criado (pelo menos conceitualmente) em torno de cada ponto do recurso. A união desses discos é o buffer.

Essa pergunta pede o cálculo de uma dilatação um pouco mais complicada, onde a propagação é permitida apenas dentro de um determinado intervalo de ângulos (isto é, dentro de um setor circular). Isso faz sentido apenas para recursos que não envolvem uma superfície consideravelmente curva (como pequenos recursos na esfera ou elipsóide ou quaisquer recursos no plano). Quando estamos trabalhando no plano, também é significativo orientar todos os setores na mesma direção. (Se estivéssemos modelando a propagação do fogo pelo vento, gostaríamos que os setores fossem orientados com o vento e seus tamanhos também pudessem variar com a velocidade do vento: essa é uma motivação para a definição geral de dilatação que dei. ) (Em superfícies curvas como um elipsóide, é impossível, em geral, orientar todos os setores na direção "mesma".)

Nas seguintes circunstâncias, a dilatação é relativamente fácil de calcular:

O recurso está no plano (ou seja, estamos dilatando um mapa do recurso e, com sorte, o mapa é bastante preciso).
A dilatação será constante : a propagação em todos os pontos do elemento ocorrerá em vizinhanças congruentes da mesma orientação.
Este bairro comum é convexo. A convexidade simplifica e acelera muito os cálculos.

Essa questão se encaixa nessas circunstâncias especializadas: pede a dilatação de polígonos arbitrários por setores circulares cujas origens (os centros dos discos de onde vieram) estão localizadas nos pontos de base. Desde que esses setores não ultrapassem 180 graus, serão convexos. (Os setores maiores sempre podem ser divididos ao meio em dois setores convexos; a união das duas dilatações menores dará o resultado desejado.)

Implementação

Como estamos realizando cálculos euclidianos - fazendo a expansão no plano - podemos dilatar um ponto apenas traduzindo a vizinhança da dilatação para esse ponto. (Para fazer isso, o bairro precisa de uma origemisso corresponderá ao ponto base. Por exemplo, a origem dos setores nesta questão é o centro do círculo a partir do qual eles são formados. Essa origem está na fronteira do setor. Na operação padrão de buffer GIS, a vizinhança é um círculo com origem no centro; agora a origem está no interior do círculo. Escolher uma origem não é um grande problema computacionalmente, porque uma mudança de origem apenas muda toda a dilatação, mas pode ser um grande problema em termos de modelagem de fenômenos naturais. A sectorfunção no código abaixo ilustra como uma origem pode ser especificada.)

A dilatação de um segmento de linha pode ser complicada, mas para uma vizinhança convexa, podemos criar a dilatação como a união das dilatações dos dois pontos finais, juntamente com um paralelogramo cuidadosamente escolhido. (No interesse do espaço, não pararei para provar afirmações matemáticas como essa, mas encorajo os leitores a tentar suas próprias provas, porque é um exercício perspicaz.) Aqui está uma ilustração usando três setores (mostrados em rosa). Eles têm raios unitários e seus ângulos são dados nos títulos. O próprio segmento de linha tem comprimento 2, é horizontal e é mostrado em preto:

Dilatações de segmento

Os paralelogramos são encontrados localizando os pontos cor de rosa que estão o mais longe possível do segmento apenas na direção vertical . Isso fornece dois pontos inferiores e dois pontos superiores ao longo de linhas paralelas ao segmento. Nós apenas temos que juntar os quatro pontos em um paralelogramo (mostrado em azul). Observe, à direita, como isso faz sentido, mesmo quando o setor em si é apenas um segmento de linha (e não um polígono verdadeiro): lá, todos os pontos do segmento foram traduzidos em uma direção 171 graus a leste do norte por uma distância que variava de 0 a 1. O conjunto desses pontos de extremidade é o paralelogramo mostrado. Os detalhes desse cálculo aparecem na bufferfunção definida no dilate.edgescódigo abaixo.

Para dilatar uma polilinha , formamos a união das dilatações dos pontos e segmentos que a formam. As duas últimas linhas de dilate.edgesexecutar esse loop.

A dilatação de um polígono requer a inclusão do interior do polígono juntamente com a dilatação de seus limites. (Essa afirmação faz algumas suposições sobre a vizinhança de dilatação. Uma é que todas as vizinhanças contêm o ponto (0,0), o que garante que o polígono seja incluído em sua dilatação. No caso de vizinhanças variáveis, também assume que a dilatação de qualquer interior O ponto do polígono não se estende além da dilatação dos pontos de contorno. É o caso de vizinhanças constantes.)

Vejamos alguns exemplos de como isso funciona, primeiro com um nonagon (escolhido para revelar detalhes) e depois com um círculo (escolhido para corresponder à ilustração na pergunta). Os exemplos continuarão usando os mesmos três bairros, mas encolheram até um raio de 1/3.

Dilatações de um nonágono

Nesta figura, o interior do polígono é cinza, as dilatações pontuais (setores) são rosa e as dilatações das bordas (paralelogramos) são azuis.

Dilatações de um círculo

O "círculo" é realmente apenas 60 gon, mas se aproxima muito bem de um círculo.

atuação

Quando o recurso base é representado por N pontos e a vizinhança de dilatação por M pontos, esse algoritmo requer esforço de O (N M) . Ele precisa ser seguido simplificando a bagunça de vértices e arestas na união, o que pode exigir esforço O (N M log (N M)): isso é algo a ser solicitado ao SIG; não deveríamos ter que programar isso.

O esforço computacional pode ser aprimorado para O (M + N) para recursos de base convexos (porque você pode descobrir como percorrer o novo limite mesclando adequadamente as listas de vértices que descrevem os limites das duas formas originais). Isso também não precisaria de limpeza posterior.

Quando a vizinhança da dilatação muda lentamente de tamanho e / ou orientação à medida que você avança no recurso base, a dilatação da aresta pode ser aproximada do casco convexo da união das dilatações de seus pontos finais. Se as duas vizinhanças de dilatação tiverem pontos M1 e M2, isso pode ser encontrado com esforço O (M1 + M2) usando um algoritmo descrito em Shamos & Preparata, Geometria Computacional . Portanto, deixando K = M1 + M2 + ... + M (N) ser o número total de vértices nas vizinhanças de dilatação de N, podemos calcular a dilatação no tempo de O (K * log (K)).

Por que quereríamos lidar com essa generalização se tudo o que queremos é um buffer simples? Para grandes recursos na Terra, uma vizinhança de dilatação (como um disco) que possui tamanho constante na realidade pode ter tamanho variável no mapa em que esses cálculos são realizados. Assim, obtemos uma maneira de fazer cálculos precisos para o elipsóide , continuando a desfrutar de todas as vantagens da geometria euclidiana.

Código

Os exemplos foram produzidos com este Rprotótipo, que pode ser facilmente transportado para a sua linguagem favorita (Python, C ++, etc.). Na estrutura, é paralela à análise relatada nesta resposta e, portanto, não precisa de explicação separada. Os comentários esclarecem alguns dos detalhes.

(Pode ser interessante observar que os cálculos trigonométricos são usados apenas para criar os recursos de exemplo - que são polígonos regulares - e os setores. Nenhuma parte dos cálculos de dilatação requer trigonometria.)

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

O tempo de computação para este exemplo (da última figura), com N = 60 e M = 121 (esquerda), M = 181 (meio) e M = 2 (direita), foi de um quarto de segundo. No entanto, a maior parte disso foi para a exibição. Normalmente, esse Rcódigo processa cerca de N M = 1,5 milhão por segundo (levando apenas 0,002 segundos ou mais para fazer todos os cálculos de exemplo mostrados). No entanto, a aparência do produto M N implica que dilatações de muitas figuras ou figuras complicadas por uma vizinhança detalhada podem levar um tempo considerável, portanto, cuidado! Avalie o momento para problemas menores antes de enfrentar um grande problema. Em tais circunstâncias, pode-se procurar uma solução baseada em varredura (que é muito mais fácil de implementar, exigindo essencialmente apenas um único cálculo de vizinhança).

whuber
fonte

Uau, isso é muito detalhado e fascinante. Eu não esperava menos.

Paul