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

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.

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

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

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!

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