Determinando se o ponto é circundado usando o processamento de varredura

Estou tentando melhorar um processo de vetor / python atualmente extremamente complicado para um modelo de risco natural. No momento, temos um longo script que gera distância / linhas de rumo a partir de um determinado ponto para determinar:

1. o tipo de polígono que cruza (por exemplo: floresta, grama, pântano, etc.)
2. a distância para esse polígono
3. quantas dessas linhas cruzam polígonos, para determinar o quão 'cercado' é.

Há muito mais envolvido, mas essa é a essência. Estou tentando encontrar uma maneira de melhorar isso e estou atualmente perplexo na parte 3. A idéia é determinar se um ponto está completamente cercado por polígonos, a cerca de 200m

Portanto, na minha imagem anexada, eu gostaria que o ponto A fosse marcado como de maior risco que o ponto B, pois está completamente cercado pelos meus polígonos. Isso é repetido por cerca de 13 milhões de pontos, portanto não é uma tarefa pequena e eu preferiria ter uma superfície para derivar valores, em vez de executar nosso script. Estou pensando que deve haver uma variedade de ferramentas de hidrologia ou caminho de custo para fazer isso, mas não consigo entender o que estou pensando.

Como eu pude fazer isso?

Existe uma solução de caminho de custo, mas você precisará codificá-la. Aqui está o que pode parecer quando aplicado a todos os pontos da imagem na pergunta (um pouco grosseiro para acelerar os cálculos):

As células pretas são partes dos polígonos circundantes. As cores, variando de laranja claro (curto) a azul (longo), mostram a distância máxima (até um máximo de 50 células) que pode ser alcançada pela passagem da linha de visão sem interceptar as células poligonais. (Qualquer célula fora da extensão desta imagem é tratada como parte dos polígonos.)

Vamos discutir uma maneira eficiente de fazer isso usando uma representação rasterizada dos dados. Nesta representação, todas as células poligonais "circundantes" terão, digamos, valores diferentes de zero e qualquer célula que possa ser "vista através" terá um valor zero.

Etapa 1: pré-computando uma estrutura de dados de vizinhança

Você primeiro precisa decidir o que significa para uma célula bloquear outra. Uma das regras mais justas que posso encontrar é a seguinte: usando coordenadas integrais para linhas e colunas (e assumindo células quadradas), vamos considerar quais células podem bloquear a célula (i, j) da vista na origem (0,0). Nomeio a célula (i ', j') que está mais próxima do segmento de linha que conecta (i, j) a (0,0) entre todas as células cujas coordenadas diferem de iej por no máximo 1. Porque isso nem sempre Para obter uma solução única (por exemplo, com (i, j) = (1,2) ambos (0,1) e (1,1) funcionarão igualmente bem), são necessários alguns meios para resolver os vínculos. Seria bom para essa resolução de vínculos respeitar as simetrias de vizinhanças circulares em grades: negar coordenadas ou alternar as coordenadas preserva essas vizinhanças. Portanto, podemos decidir quais células bloqueiam (i,

Ilustrando esta regra está o seguinte código de protótipo escrito em R. Esse código retorna uma estrutura de dados que será conveniente para determinar o "ambiente" de células arbitrárias em uma grade.

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

O valor de screen(12)foi usado para produzir essa representação dessa relação de triagem: as setas apontam das células para as que as travam imediatamente. Os tons são proporcionais à distância da origem, que fica no meio desse bairro:

Esse cálculo é rápido e precisa ser feito apenas uma vez para um determinado bairro. Por exemplo, ao observar 200 m em uma grade com células de 5 m, o tamanho da vizinhança será 200/5 = 40 unidades.

Etapa 2: aplicando a computação aos pontos selecionados

O restante é direto: para determinar se uma célula localizada em (x, y) (nas coordenadas da linha e da coluna) está "cercada" com relação a essa estrutura de dados da vizinhança, execute o teste recursivamente começando com um deslocamento de (i, j) = (0,0) (a origem do bairro). Se o valor na grade do polígono em (x, y) + (i, j) for diferente de zero, a visibilidade será bloqueada lá. Caso contrário, teremos que considerar todas as compensações que poderiam ter sido bloqueadas no deslocamento (i, j) (que são encontradas no tempo O (1) usando a estrutura de dados retornada por screen). Se não houver nenhum bloqueado, alcançamos o perímetro e concluímos que (x, y) não está cercado; portanto, paramos o cálculo (e não nos importamos em inspecionar quaisquer pontos restantes na vizinhança).

Podemos coletar informações ainda mais úteis, acompanhando a maior distância da linha de visão alcançada durante o algoritmo. Se este for menor que o raio desejado, a célula será cercada; caso contrário, não é.

Aqui está um Rprotótipo desse algoritmo. É mais longo do que parece, porque Rnão suporta nativamente a estrutura de pilha (simples) necessária para implementar a recursão; portanto, uma pilha também precisa ser codificada. O algoritmo real começa cerca de dois terços do caminho e precisa de apenas uma dúzia de linhas. (E metade deles simplesmente lida com a situação na borda da grade, verificando índices fora do intervalo dentro da vizinhança. Isso poderia ser mais eficiente simplesmente expandindo a grade do polígono por klinhas e colunas em torno de seu perímetro, eliminando qualquer é necessário verificar o intervalo do índice com o custo de um pouco mais de RAM para manter a grade de polígonos.)

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

Neste exemplo, as células poligonais são pretas. As cores fornecem a distância máxima da linha de visão (até 50 células) para células não poligonais, variando de laranja claro para distâncias curtas a azul escuro nas distâncias mais longas. (As células têm uma unidade de largura e altura.) As faixas visivelmente evidentes são criadas pelas pequenas "ilhas" poligonais no meio do "rio": cada uma delas bloqueia uma longa fila de outras células.

Análise do algoritmo

A estrutura da pilha implementa uma pesquisa profunda do gráfico de visibilidade da vizinhança em busca de evidências de que uma célula não está cercada. Onde as células estão longe de qualquer polígono, essa pesquisa exigirá a inspeção de apenas células O (k) para uma vizinhança circular de raio-k. Os piores casos ocorrem quando há um pequeno número de células poligonais dispersas na vizinhança, mas mesmo assim o limite da vizinhança não pode ser alcançado: requer inspeção de quase todas as células de cada vizinhança, que é um O (k ^ 2) Operação.

O comportamento a seguir é típico do que será encontrado. Para pequenos valores de k, a menos que os polígonos preencham a maior parte da grade, a maioria das células não poligonais será obviamente sem volta e o algoritmo será escalado como O (k). Para valores intermediários, a escala começa com O (k ^ 2). À medida que k se torna realmente grande, a maioria das células é cercada e esse fato pode ser determinado bem antes de toda a vizinhança ser inspecionada: o esforço computacional do algoritmo atinge, assim, um limite prático. Esse limite é atingido quando o raio da vizinhança se aproxima do diâmetro das maiores regiões não poligonais conectadas na grade.

Como exemplo, usei a countingopção codificada no protótipo de screenpara retornar o número de operações de pilha usadas em cada chamada. Isso mede o esforço computacional. O gráfico a seguir mostra o número médio de operações da pilha em função do raio da vizinhança. Ele exibe o comportamento previsto.

Podemos usar isso para estimar o cálculo necessário para avaliar 13 milhões de pontos em uma grade. Suponha que uma vizinhança de k = 200/5 = 40 seja usada. Então, algumas centenas de operações de pilha serão necessárias, em média (dependendo da complexidade da grade de polígonos e da localização dos 13 milhões de pontos em relação aos polígonos), o que implica que, em uma linguagem compilada eficiente, no máximo alguns milhares de operações numéricas simples será necessário (adicionar, multiplicar, ler, escrever, compensar etc.). A maioria dos PCs será capaz de avaliar a área circundante de um milhão de pontos nesse ritmo. (ORa implementação é muito, muito mais lenta que isso, porque é pobre nesse tipo de algoritmo, e é por isso que só pode ser considerado um protótipo.) Dessa forma, podemos esperar que uma implementação eficiente em uma linguagem razoavelmente eficiente e apropriada - C ++ e Python vem à mente - poderia concluir a avaliação de 13 milhões de pontos em um minuto ou menos, assumindo que toda a grade de polígonos reside na RAM.

Quando uma grade é muito grande para caber na RAM, esse procedimento pode ser aplicado às partes lado a lado da grade. Eles só precisam se sobrepor por klinhas e colunas; tome os máximos nas sobreposições ao aplicar mosaicos nos resultados.

Outras aplicações

A "busca" de um corpo de água está intimamente relacionada à "envolvência" de seus pontos. De fato, se usarmos um raio de vizinhança igual ou maior que o diâmetro do corpo de água, criaremos uma grade da busca (não direcional) em todos os pontos do corpo de água. Usando um raio de vizinhança menor, obteremos pelo menos um limite inferior para a busca em todos os pontos de busca mais altos, o que em algumas aplicações pode ser bom o suficiente (e pode reduzir substancialmente o esforço computacional). Uma variante desse algoritmo que limita a relação "rastreada por" a direções específicas seria uma maneira de calcular a busca com eficiência nessas direções. Observe que essas variantes exigem a modificação do código para screen; o código para panvisibilitynão muda nada.

Definitivamente, posso ver como alguém pode querer fazer isso com uma solução raster, mas, mesmo com um número reduzido de pontos, esperaria uma resolução muito grande / alta e, portanto, difícil de processar ou conjunto de grades. Dado isso, gostaria de saber se explorar a topologia em um gdb pode ser mais eficiente. Você pode encontrar todos os vazios internos com 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')`

então você pode trabalhar com topoErrorsVoidPolysseu padrão normal Intersect_analysis()ou o que quer que seja. Você pode precisar mexer na extração de políticas de interesse topoErrorsVoidPolys. O @whuber tem várias postagens excelentes sobre esse tipo de coisa em outros lugares aqui no gis.stackexchange.com.

Se você realmente quer ficar mais raster ... Eu faria algo como esse pseudo-código (não se estremece só porque é óbvio que sou um retrocesso de AML!: P)

1. pontos de rasterização ("pts_g") e polys ("polys_g" (
2. voids = regiongroup (con (isnull (polys_g), 1))
3. pode ser necessário fazer algo para refinar os vazios para eliminar polígonos externos indesejados / área de universo aberto
4. pts_surrounded = con (vazios, pts_g)

Apenas inventando isso, pode precisar de refinamento.

Se você usar um valor de distância limite (aqui você fala sobre 200 m), a melhor solução é usar a análise vetorial:

1) crie um buffer de 200 m em torno de cada ponto (em preto na ilustração)

2) use a ferramenta de interseção (análise) entre o buffer e os polígonos (em azul na ilustração). Ficará melhor se você fizer isso entre os limites dos polígonos circundantes e o buffer, mas o resultado final será o mesmo.

3) use o recurso para polígono (gerenciamento) para criar polígonos onde seus pontos estão completamente cercados (em vermelho na ilustração)

4) selecione camadas por localização (gerenciamento) ou junção espacial (análise) para identificar os pontos que estão circundados. O uso da junção espacial permite que você tenha informações sobre o polígono de incorporação (área do polígono, estatísticas zonais ...) que podem ser úteis para processamento adicional.

Alternativas 2b) Dependendo de suas necessidades, você pode selecionar por localização os polígonos ao redor a uma distância de 200 m, e então identificar alguns tipos de "compartimento", mas não tão estritamente quanto em 2).

Considerando o "caso do labirinto", isso pode ajudar: avaliar quanto tempo é necessário "escapar" do local.

• Você já pode excluir da análise os pontos que estão totalmente incluídos ou totalmente gratuitos

• então você converte seus obstáculos em uma varredura e define os valores para NoData onde você possui um polígono e para o tamanho da célula em metros onde você não possui (isso aumentará o custo).

• terceiro, você pode calcular a distância de custo usando o recém-gerado custo raster

• finalmente, você usa estatísticas zonais como tabela com base nos limites do buffer convertidos em raster (formando um espaço anular). Se você puder escapar em todas as direções, o mínimo deve ser aproximadamente 200 (dependendo do tamanho da célula da sua análise). Mas se você estiver em um labirinto, o máximo será maior que 200. Portanto, o máximo das estatísticas zonais menos 200 será um valor contínuo, indicando o quão difícil é "escapar".