Identificación de puntos consecutivos dentro de un buffer especificado

Tengo un archivo de puntos, reubicaciones por hora de un animal, y quiero poder colocar un búfer alrededor de cada punto y calcular el número de puntos posteriores que están dentro de la zona de búfer. Estoy buscando un método que funcione a lo largo del archivo de puntos, como una ventana en movimiento, que solo contará los puntos posteriores que están dentro de la zona de amortiguamiento.

Por ejemplo, en el punto 10 coloco un búfer de 1500 my quiero saber si el punto 11 está dentro del búfer y, en caso afirmativo, si el punto 12 está dentro del búfer y así sucesivamente. No quiero saber si el punto 52 está dentro de la zona de seguridad a menos que todos los puntos anteriores hayan estado dentro de la memoria intermedia. Tampoco quiero saber si los puntos 9 u 8, etc. están dentro del búfer.

Encontré y probé una caja de herramientas adicional llamada "análisis de puntos de ventana móvil" que funciona como una ventana móvil en el archivo de puntos. Esto funciona bien, pero muy lentamente, e incluye todos los puntos que están dentro de la zona de amortiguamiento, incluso si no son puntos consecutivos. No puedo encontrar una manera de hacer que mire puntos consecutivos.

Me gustaría un método que proporcione una tabla de salida, ya que tengo muchos puntos de datos para ver de esta manera.

Estoy usando ArcGIS 10. Cualquier ayuda que alguien pueda proporcionar sería muy apreciada.

Dada una lista de ubicaciones de puntos (preferiblemente en coordenadas proyectadas, para que las distancias sean fáciles de calcular), este problema se puede resolver con cinco operaciones más simples :

1. Calcular distancias punto-punto.

2. Para cada punto i, i = 1, 2, ..., identifique los índices de esos puntos a distancias inferiores al radio del búfer (como 1500).

3. Restrinja esos índices para que sean i o mayores.

4. Retenga solo el primer grupo consecutivo de índices sin interrupción.

5. Salida de la cuenta de ese grupo.

En R, cada uno de estos corresponde a una operación. Para aplicar esta secuencia a cada punto, es conveniente encapsular la mayor parte del trabajo dentro de una función que definimos , por lo tanto:

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

(Consulte a continuación una versión más eficiente de esta función).

He hecho esta función lo suficientemente flexible como para aceptar varias listas de puntos ( xy) y distancias de búfer ( r) como parámetros.

Normalmente, leería un archivo de ubicaciones de puntos (y, si fuera necesario, los ordenaría por hora). Aquí, para mostrar esto en acción, solo generaremos algunos datos de muestra al azar :

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

Su espaciado típico es 300 * Sqrt (2) = aproximadamente 500. Hacemos el cálculo aplicando esta función a los puntos en la matrizxy (y luego añadiendo sus resultados nuevamente xy, porque este sería un formato conveniente para exportar a un SIG ):

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

Luego analizaría más a fondo la resultmatriz, ya sea Rescribiéndola en un archivo e importándola a otro software. Aquí está el resultado para los datos de muestra :

                        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

(Recuerde que los recuentos incluyen los puntos en los que se basan, por lo que cada recuento debe ser 1 o mayor).

Si tiene muchos miles de puntos, este método es demasiado ineficiente : calcula demasiadas distancias punto a punto que son innecesarias. Pero debido a que hemos encapsulado el trabajo dentro de la forwardfunción, la ineficiencia es fácil de solucionar. Aquí hay una versión que funcionará mejor cuando estén involucrados más de unos pocos cientos de puntos:

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 probar esto, creé puntos aleatorios como anteriormente, pero varié dos parámetros: n(el número de puntos) y su desviación estándar (codificada como 300 arriba). La desviación estándar determina el número promedio de puntos dentro de cada búfer ("promedio" en la tabla a continuación): cuantos más haya, más tardará en ejecutarse este algoritmo. (Con algoritmos más sofisticados, el tiempo de ejecución no dependerá tanto de cuántos puntos haya en cada búfer). Aquí hay algunos tiempos:

 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

Creo que su mejor opción es crear una pequeña rutina con ArcPy. Crearía algo como este pseudocó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'

No estoy seguro de qué desea hacer con la información, pero supongo que podría crear un campo en su tabla y actualizarlo con el recuento de puntos consecutivos (si es así, agregue el campo primero).

Recomendaría crear capas de entidades (como una vista de tabla de base de datos pero para entidades en Arc). Haga dos a partir de los datos originales y abra un cursor de actualización en el primer conjunto que especifique su clasificación general (porque ESRI no acepta consultas SQL completas). Use el segundo para seleccionar por ubicación y abra un Cursor de búsqueda en el conjunto de selección resultante.

[EDITAR DE ACUERDO CON LA SOLICITUD DE JAMÉ] Ábrelo usando Model Builder. Si no ha utilizado Model Builder antes, todo lo que tiene que hacer es hacer clic derecho en arcToolbox. Seleccione 'Agregar cuadro de herramientas'. Haga clic derecho en la nueva caja de herramientas y haga clic en 'Nuevo-> modelo'. Una vez que tenga una nueva ventana de modelo, arrastre y suelte las herramientas y los datos que necesita en la ventana y vincúlelos visualmente (usando la pequeña herramienta de flecha). Cuando haya llegado lo más lejos posible (no podrá agregar sus cursores aquí), use la opción en el menú Archivo del Generador de modelos para exportar a Python. Eso te llevará a la mayor parte del camino. Es un código generado automáticamente, por lo que será un poco desagradable pero funcional. Luego use los enlaces en mi respuesta anterior para comprender y agregar los cursores.

Si eres nuevo en Python, ¡no tengas miedo de escribir código! Python es un lenguaje de script muy fácil para obtener resultados rápidamente. Esri también tiene alguna guía al respecto.

Si se queda atascado con su código, publíquelo en este foro y solicite ayuda. Hay MUCHAS personas aquí que pueden ayudar. Una advertencia: asegúrese de utilizar la ayuda correcta de ESRI. Cambiaron masivamente la API de Python entre las versiones 9.xy 10 (respectivamente arcgisscripting y arcpy). Entonces, si está utilizando ArcGIS 9.x, ¡encuentre los enlaces equivalentes a los míos!

Puede usar el generador de modelos en ArcGIS para buscar valores de ID consecutivos. Exporté mi modelo como un script de Python. El código generará un nuevo shp que tiene valores de ID consecutivos. !¡CARNÉ DE IDENTIDAD! es el campo de ID base. Deberá actualizar la ruta point2.shp, el nombre y el nombre del campo ID para que coincida con su 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")