Identificação de pontos consecutivos dentro de um buffer especificado

8

Eu tenho um arquivo de pontos, realocações horárias de um animal e quero poder colocar um buffer em torno de cada ponto e calcular o número de pontos subseqüentes que estão dentro da zona de buffer. Estou procurando um método que funcione ao longo do arquivo de ponto, como uma janela em movimento, que contará apenas os pontos subsequentes que estão dentro da zona de buffer.

Por exemplo, no ponto 10, coloco um buffer de 1500m e quero saber se o ponto 11 está dentro do buffer e, se sim, então se o ponto 12 está dentro do buffer e assim por diante. Não quero saber se o ponto 52 está dentro da zona do buffer, a menos que todos os pontos anteriores estejam dentro do buffer. Também não quero saber se os pontos 9 ou 8 etc. estão dentro do buffer.

Eu encontrei e tentei uma caixa de ferramentas adicional chamada "análise de pontos em movimento da janela", que funciona como uma janela em movimento no arquivo do ponto. Isso funciona bem, mas muito lentamente, e inclui todos os pontos que estão dentro da zona de buffer, mesmo que não sejam pontos consecutivos. Não consigo encontrar uma maneira de fazer com que pareça pontos consecutivos.

Gostaria de um método que forneça uma tabela de saída, pois tenho muitos pontos de dados para analisar dessa maneira.

Estou usando o ArcGIS 10. Qualquer ajuda que alguém possa fornecer seria muito apreciada.

James
fonte
Seus pontos provavelmente se originaram como uma lista de (x, y, hora) dados. Você estaria aberto para pré-processar esses dados (fora do ArcGIS) para obter as informações desejadas?
whuber
Se isso torna mais fácil, então definitivamente. Também estou processando os dados usando o AdehabitatLT no R, que calcula distância e rumo, etc. Entendo o processo sugerido por Sylvester abaixo, mas estou lutando para saber por onde começar, pois não tenho muita certeza de quais ferramentas preciso usar.
James
Ah! Como você já está usando R, vamos explorar soluções baseadas em R então.
whuber
Existe uma função de janela deslizante no AdehabitatLT "sliwinltr", mas não sei como usá-lo neste caso. Nem sei se pode ser usado dessa maneira.
James

Respostas:

8

Dada uma lista de locais dos pontos (de preferência nas coordenadas projetadas, para facilitar a computação das distâncias), esse problema pode ser resolvido com cinco operações mais simples :

  1. Calcular distâncias ponto a ponto.

  2. Para cada ponto i, i = 1, 2, ..., identifique os índices desses pontos a distâncias menores que o raio do buffer (como 1500).

  3. Restrinja esses índices a serem i ou superiores.

  4. Reter apenas o primeiro grupo consecutivo de índices sem interrupção.

  5. Saída a contagem desse grupo.

Em R, cada um deles corresponde a uma operação. Para aplicar essa sequência a cada ponto, é conveniente encapsular a maior parte do trabalho dentro de uma função que definimos , assim:

#
# forward(j, xy, r) counts how many contiguous rows in array xy, starting at index j,
#                   are within (Euclidean) distance r of the jth row of xy.
#
forward <- function(j, xy, r) {
  # Steps 1 and 2: compute an array of indexes of points within distance r of point j.
  i <- which(apply(xy, 1, function(x){sum((x-xy[j,])^2) <= r^2}))
  # Step 3: select only the indexes at or after j.
  i <- i[i >= j]
  # Steps 4 and 5: retain only the first consecutive group and count it.
  length(which(i <= (1:length(i) + j)))
}

(Veja abaixo uma versão mais eficiente dessa função.)

Tornei essa função flexível o suficiente para aceitar várias listas de pontos ( xy) e distâncias de buffer ( r) como parâmetros.

Normalmente, você lia um arquivo de locais dos pontos (e, se necessário, os classifica por tempo). Aqui, para mostrar isso em ação, apenas geraremos alguns dados de amostra aleatoriamente :

# Create sample data
n<-16                                     # Number of points
set.seed(17)                              # For reproducibility
xy <- matrix(rnorm(2*n) + 1:n, n, 2) * 300
#
# Display the track.
plot(xy, xlab="x", ylab="y")
lines(xy, col="Gray")

Figura

Seu espaçamento típico é de 300 * Sqrt (2) = cerca de 500. Fazemos o cálculo aplicando essa função aos pontos da matrizxy (e depois voltando a seus resultados xy, pois esse seria um formato conveniente para exportar para um GIS ):

radius <- 1500
z <- sapply(1:n, function(u){forward(u,xy,radius)})
result <- cbind(xy, z)                              # List of points, counts

Você analisaria ainda mais a resultmatriz, Rgravando-a em um arquivo e importando-a para outro software. Aqui está o resultado para os dados da amostra :

                        z
  [1,]   -4.502615  551.5413 4
  [2,]  576.108979  647.8110 3
  [3,]  830.103893 1087.7863 4
  [4,]  954.819620 1390.0754 3
...
 [15,] 4977.361529 4146.7291 2
 [16,] 4783.446283 4511.9500 1

(Lembre-se de que as contagens incluem os pontos em que se baseiam, de modo que cada contagem deve ser 1 ou maior.)


Se você tem muitos milhares de pontos, esse método é muito ineficiente : calcula distâncias ponto a ponto demais desnecessárias. Mas, como encapsulamos o trabalho dentro da forwardfunção, a ineficiência é fácil de corrigir. Aqui está uma versão que funcionará melhor quando mais de algumas centenas de pontos estiverem envolvidos:

forward <- function(j, xy, r) {
   n <- dim(xy)[1]     # Limit the search to the number of points in xy
   r2 <- r^2           # Pre-compute the squared distance threshold
   z <- xy[j,]         # Pre-fetch the base point coordinates
   i <- j+1            # Initialize an index into xy (just past point j)

   # Advance i while point i remains within distance r of point j.
   while(i <= n && sum((xy[i,]-z)^2) <= r2) i <- i+1

   # Return the count (including point j).
   i-j
}

Para testar isso, criei pontos aleatórios como anteriormente, mas variei dois parâmetros: n(o número de pontos) e seu desvio padrão (codificado como 300 acima). O desvio padrão determina o número médio de pontos dentro de cada buffer ("média" na tabela abaixo): quanto mais houver, mais tempo esse algoritmo demora para ser executado. (Com algoritmos mais sofisticados, o tempo de execução não depende tanto de quantos pontos existem em cada buffer.) Aqui estão alguns horários:

 Time (sec)    n    SD  Average  Distances checked per minute
 1.30       10^3     3  291      13.4 million
 1.72       10^4    30   35.7    12.5
 2.50       10^5   300    3.79    9.1
16.4        10^6  3000    1.04    3.8
whuber
fonte
Parece a solução perfeita. No entanto, o código não é executado a partir de "z <- sapply (1: n), função (u) {forward (u, xy, raio)})" como diz: "inesperado ',' em" z <- sapply (1: n)," Se eu remover a vírgula então diz: Erro: 'função' inesperado no "z <- sapply: function (1 n)" Todas as ideias por que isso pode ser assim?
James
Desculpa; há um erro de digitação lá: vou remover o ")" estranho. (Fiz uma simplificação de última hora desse código. Ele foi testado mais vezes do que gostaria de admitir!)
whuber
1
Isso é ótimo, está funcionando agora. Acabei de carregar meus dados como xy para simplificar a verificação de que funciona. Demora um pouco para ser executado como você mencionou, mas parece ter feito tudo corretamente. Vou verificar manualmente alguns com o meu mapa GIS, mas parece bom até agora. Obrigado por me ajudar a resolver isso, muito interessado em aprender em GIS e R e estou na íngreme curva de aprendizado.
James
1
Editei a resposta para fornecer uma solução com escalabilidade muito aprimorada. Agora ele é capaz de lidar com caminhos contendo milhões de pontos.
whuber
1
Executei o código original com arquivos de pontos de 2000 entradas, que levavam algumas horas cada, como você diz que havia muitos locais ponto a ponto não utilizados, prolongando o tempo de processamento. A edição acima parece uma solução interessante e tentarei isso com os mesmos dados e ver quanto mais rápido é. Obrigado pelo esforço em produzir e editar a função.
James
1

Acho que sua melhor aposta é criar um pouco de rotina usando o ArcPy. Eu criaria algo como este pseudo-código:

Select all points sorted by location-id    
Iterate for each point:
    Select points by location using a distance (no need to create buffer geometry)
    Sort points by location-id
    Set a variable to the value of your reference id
    consecutive-counter = 0
        Iterate your selection:
            Is the location-id of the first (or next) point equal to variable + 1?
            if 'yes' increment consecutive-counter by 1
            Repeat until 'no'

Não tenho certeza do que você deseja fazer com as informações, mas suponho que você possa criar um campo em sua tabela e atualizá-lo com a contagem de pontos consecutivos (se houver, adicione o campo primeiro).

Eu recomendaria criar camadas de recursos (como uma exibição de tabela de banco de dados, mas para recursos no Arc). Faça dois a partir dos dados originais e abra um cursor de atualização no primeiro conjunto, especificando sua classificação geral (porque o ESRI não atende a consultas SQL completas). Use o segundo para selecionar por local e abrir um Cursor de Pesquisa no conjunto de seleção resultante.

[EDITAR POR PEDIDO DO Jame] Faça o detalhamento usando o Model Builder. Se você não usou o Model Builder antes, basta clicar com o botão direito do mouse no arcToolbox. Selecione 'Adicionar caixa de ferramentas'. Clique com o botão direito do mouse na nova caixa de ferramentas e clique em 'Novo-> modelo'. Depois de ter uma nova janela de modelo, arraste e solte as ferramentas e os dados necessários na janela e vincule-os visualmente (usando a pequena ferramenta de seta). Quando você chegar o mais longe possível (não poderá adicionar seus cursores aqui), use a opção no menu Arquivo do Model Builder para exportar para o Python. Isso o levará a maior parte do caminho até lá. Como o código é gerado automaticamente, será um pouco desagradável, mas funcional. Em seguida, use os links na minha resposta acima para entender e adicionar os cursores.

Se você é novo no Python, não tenha medo de escrever código! Python é uma linguagem de script muito fácil de obter resultados rapidamente. Esri também tem algumas orientações .

Se você ficar preso ao seu código, publique-o neste fórum e peça ajuda. Há muitas pessoas aqui que podem ajudar. Um aviso - certifique-se de usar a ajuda certa da ESRI. Eles mudaram massivamente a API do Python entre as versões 9.xe 10 (respectivamente arcgisscripting e arcpy). Portanto, se você estiver usando o ArcGIS 9.x, encontre os links equivalentes aos meus!

MappaGnosis
fonte
Parece o que eu gostaria de fazer. No entanto, atualmente não estou usando código no ArcGIS, apenas selecionando opções pré-definidas. Como eu começaria a usar / gerar o método sugerido acima? Gostaria que a saída fosse uma nova tabela ou um novo campo adicionado à tabela com a contagem de pontos consecutivos.
James
Veja minha edição no meu post principal.
MappaGnosis
JTB, faça logon usando a mesma conta que você usou para postar esta pergunta para poder postar comentários. (Para facilitar,
fundei
Desculpe pela alteração da conta. Eu postei a pergunta original como um novo usuário, mas não consegui voltar para a conta porque não tinha uma senha etc. Por isso, criei outra conta JTB que usarei a partir de agora (espero). Iniciarei o construtor de modelos conforme sugerido por Sylvester, mas nunca o usei antes, pode demorar um pouco para me acostumar e para descobrir quais ferramentas etc. usar. Voltarei com progresso e perguntas. Obrigado
James
Sylvester - Acho que entendi o processo, mas não sei com quais ferramentas realmente preciso começar. Distância? Amortecedor? Perto? Nem sei se existe uma ferramenta correta para esse problema que fará o que foi mencionado acima. Estou interessado em aprender, mas estou muito no começo.
James
1

Você pode usar o construtor de modelos no ArcGIS para encontrar valores de ID consecutivos. Eu exportei meu modelo como um script python. O código irá gerar um novo shp com valores de ID consecutivos. !EU IRIA! é o campo de identificação base. Você precisará atualizar o caminho point2.shp, o nome e o nome do campo de ID para corresponder ao seu caso.

# Import arcpy module
import arcpy


# Local variables:
point2_shp = "C:\\temp\\point2.shp"
point2_shp__2_ = "C:\\temp\\point2.shp"
point2_shp__4_ = "C:\\temp\\point2.shp"
Freq_dbf__2_ = "C:\\temp\\Freq.dbf"
point2_shp__5_ = "C:\\temp\\point2.shp"
point2__2_ = "C:\\temp\\point2.shp"
point2__4_ = "C:\\temp\\point2.shp"
Freq_dbf = "C:\\temp\\Freq.dbf"
PointConsecutive_shp = "C:\\temp\\PointConsecutive.shp"

# Process: Add Field
arcpy.AddField_management(point2_shp, "AUTOID", "LONG", "", "", "", "", "NON_NULLABLE", "NON_REQUIRED", "")

# Process: Calculate Field
arcpy.CalculateField_management(point2__2_, "AUTOID", "autoIncrement()", "PYTHON_9.3", "rec=0\\ndef autoIncrement():\\n global rec\\n pStart = 1 #adjust start value, if req'd \\n pInterval = 1 #adjust interval value, if req'd\\n if (rec == 0): \\n  rec = pStart \\n else: \\n  rec = rec + pInterval \\n return rec\\n")

# Process: Add Field (2)
arcpy.AddField_management(point2_shp, "DIFF", "LONG", "", "", "", "", "NON_NULLABLE", "NON_REQUIRED", "")

# Process: Calculate Field (2)
arcpy.CalculateField_management(point2__4_, "DIFF", "!ID! - !AUTOID!", "PYTHON_9.3", "")

# Process: Frequency
arcpy.Frequency_analysis(point2_shp__2_, Freq_dbf, "DIFF", "")

# Process: Join Field
arcpy.JoinField_management(point2_shp__4_, "DIFF", Freq_dbf__2_, "DIFF", "")

# Process: Select
arcpy.Select_analysis(point2_shp__5_, PointConsecutive_shp, "\"FREQUENCY\" >1")
artwork21
fonte
Não sei ao certo o que o código acima faz e como ele responde à minha consulta. Agradeço a ajuda, mas isso tudo é bastante novo para mim, pois nunca usei o Python ou o construtor de modelos. Eu altero o campo ID para cada processo listado para o ID no conjunto de dados?
James
@ James, depende se o seu ID mudar. Para usar o código, basta copiar e colar o código acima e salve-o em um arquivo txt em branco. Atualize o código para corresponder ao seu nome e caminho point.shp. Em seguida, altere o nome do ID na seção de código Calcular campo (2) para corresponder ao seu campo point.shp ID. Salve o arquivo txt e, no Windows Explorer, renomeie o arquivo com uma extensão .py. Clique com o botão direito do mouse no arquivo e abra com python.exe para testar.
artwork21
Idealmente, esse script pode ser conectado a uma ferramenta de script, que também lida com o buffer e a seleção de recursos. help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#/...
artwork21