Identifizierung aufeinanderfolgender Punkte innerhalb eines bestimmten Puffers

8

Ich habe eine Punktedatei, stündliche Umzüge eines Tieres, und ich möchte in der Lage sein, einen Puffer um jeden Punkt zu legen und die Anzahl der nachfolgenden Punkte zu berechnen, die sich innerhalb der Pufferzone befinden. Ich suche nach einer Methode, die entlang der Punktdatei funktioniert, wie ein sich bewegendes Fenster, das nur die nachfolgenden Punkte zählt, die sich innerhalb der Pufferzone befinden.

Zum Beispiel platziere ich an Punkt 10 einen Puffer von 1500 m und möchte wissen, ob Punkt 11 innerhalb des Puffers liegt und wenn ja, ob Punkt 12 innerhalb des Puffers liegt und so weiter. Ich möchte nicht wissen, ob sich Punkt 52 innerhalb der Pufferzone befindet, es sei denn, alle vorherigen Punkte befanden sich innerhalb des Puffers. Ich möchte auch nicht wissen, ob sich die Punkte 9 oder 8 usw. im Puffer befinden.

Ich habe eine zusätzliche Toolbox namens "Moving Window Point Analysis" gefunden und ausprobiert, die als Moving Window für die Punktdatei fungiert. Dies funktioniert gut, aber sehr langsam und schließt alle Punkte ein, die sich innerhalb der Pufferzone befinden, auch wenn es sich nicht um aufeinanderfolgende Punkte handelt. Ich kann keinen Weg finden, um nur aufeinanderfolgende Punkte zu betrachten.

Ich möchte eine Methode, die eine Ausgabetabelle bereitstellt, da ich viele Datenpunkte auf diese Weise betrachten kann.

Ich verwende ArcGIS 10. Jede Hilfe, die jemand leisten kann, wäre sehr dankbar.

James
quelle
Ihre Punkte stammen wahrscheinlich aus einer Liste von (x, y, Zeit) Daten. Möchten Sie diese Daten (außerhalb von ArcGIS) vorverarbeiten, um die gewünschten Informationen zu erhalten?
whuber
Wenn das es einfacher macht, dann definitiv. Ich verarbeite die Daten auch mit AdehabitatLT in R, das Abstand und Peilung usw. berechnet. Ich verstehe den von Sylvester unten vorgeschlagenen Prozess, habe aber Schwierigkeiten zu wissen, wo ich anfangen soll, da ich nicht wirklich sicher bin, welche Werkzeuge ich verwenden muss.
James
Ah! Da Sie bereits verwenden R, lassen Sie uns dann R-basierte Lösungen untersuchen.
whuber
Es gibt eine Schiebefensterfunktion in AdehabitatLT "sliwinltr", aber ich weiß nicht, wie ich sie in diesem Fall verwenden soll. Ich weiß nicht einmal, ob es auf diese Weise verwendet werden kann.
James

Antworten:

8

Bei einer Liste von Punktpositionen (vorzugsweise in projizierten Koordinaten, damit Entfernungen leicht berechnet werden können) kann dieses Problem mit fünf einfacheren Operationen gelöst werden :

  1. Berechnen Sie Punkt-Punkt-Abstände.

  2. Identifizieren Sie für jeden Punkt i, i = 1, 2, ... die Indizes dieser Punkte in Abständen, die kleiner als der Pufferradius sind (z. B. 1500).

  3. Beschränken Sie diese Indizes auf i oder höher.

  4. Behalten Sie nur die erste aufeinanderfolgende Gruppe von Indizes bei, die keine Unterbrechung aufweisen.

  5. Geben Sie die Anzahl dieser Gruppe aus.

In entspricht Rjede dieser Operationen einer Operation. Um gelten für jeden Punkt dieser Sequenz ist es bequem, die meisten der Arbeit innerhalb einer Funktion kapseln wir definieren , also:

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

(Eine effizientere Version dieser Funktion finden Sie weiter unten.)

Ich habe diese Funktion flexibel genug gemacht, um verschiedene Punktlisten ( xy) und Pufferabstände ( r) als Parameter zu akzeptieren .

Normalerweise lesen Sie eine Datei mit Punktpositionen (und sortieren Sie sie gegebenenfalls nach Zeit). Um dies in Aktion zu zeigen, werden wir nur einige Beispieldaten zufällig generieren :

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

Zahl

Ihr typischer Abstand beträgt 300 * Sqrt (2) = ungefähr 500. Wir führen die Berechnung durch, indem wir diese Funktion auf die Punkte im Array anwendenxy (und dann die Ergebnisse wieder anheften xy, da dies ein praktisches Format für den Export in ein GIS wäre ):

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

Anschließend analysieren Sie das resultArray weiter, Rindem Sie es entweder in eine Datei schreiben oder in eine andere Software importieren. Hier ist das Ergebnis für die Beispieldaten :

                        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

(Denken Sie daran , dass die Zählungen sind die Punkte , an denen sie basieren, so dass jede Zählung 1 oder größer sein müssen.)


Wenn Sie viele tausend Punkte haben, ist diese Methode zu ineffizient : Sie berechnet viel zu viele unnötige Punkt-zu-Punkt-Abstände. Da wir jedoch die Arbeit in die forwardFunktion eingekapselt haben , lässt sich die Ineffizienz leicht beheben. Hier ist eine Version, die besser funktioniert, wenn mehr als ein paar hundert Punkte beteiligt sind:

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
}

Um dies zu testen, habe ich wie zuvor zufällige Punkte erstellt, aber zwei Parameter variiert: n(die Anzahl der Punkte) und ihre Standardabweichung (oben als 300 fest codiert). Die Standardabweichung bestimmt die durchschnittliche Anzahl von Punkten in jedem Puffer ("Durchschnitt" in der folgenden Tabelle): Je mehr Punkte vorhanden sind, desto länger dauert die Ausführung dieses Algorithmus. (Bei komplexeren Algorithmen hängt die Laufzeit nicht so stark davon ab, wie viele Punkte sich in jedem Puffer befinden.) Hier einige Zeitpunkte:

 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
quelle
Dies scheint die perfekte Lösung zu sein. Der Code läuft jedoch nicht von "z <- sapply (1: n), Funktion (u) {forward (u, xy, Radius)})", da er sagt: "unerwartet ',' in" z <- sapply (1: n), "Wenn ich das Komma entferne, heißt es dann: Fehler: unerwartete 'Funktion' in" z <- sapply (1: n) Funktion "Irgendwelche Ideen, warum dies so sein könnte?
James
Es tut uns leid; Da ist ein Tippfehler: Ich werde das überflüssige ")" entfernen. (Ich habe diesen Code in letzter Minute vereinfacht. Er wurde öfter getestet, als ich zugeben
möchte
1
Das ist großartig, es läuft jetzt. Ich habe meine Daten nur als xy geladen, um zu überprüfen, ob sie funktionieren. Es dauert ein wenig, bis es läuft, wie Sie bereits erwähnt haben, aber es scheint alles richtig gemacht zu haben. Ich werde einige manuell mit meiner GIS-Karte überprüfen, aber es sieht bisher gut aus. Vielen Dank, dass Sie mir dabei geholfen haben, das herauszufinden. Ich bin sehr daran interessiert, sowohl in GIS als auch in R zu lernen, und ich bin auf der steilen Lernkurve.
James
1
Ich habe die Antwort bearbeitet, um eine Lösung mit stark verbesserter Skalierbarkeit bereitzustellen. Es ist jetzt in der Lage, Pfade mit Millionen von Punkten zu verarbeiten.
whuber
1
Ich habe den ursprünglichen Code mit Punktdateien mit 2000 Einträgen ausgeführt, was jeweils einige Stunden dauerte, da Sie sagen, dass es viel zu viele nicht verwendete Punkt-zu-Punkt-Speicherorte gab, die die Verarbeitungszeit verlängerten. Die obige Bearbeitung sieht nach einer ordentlichen Lösung aus, und ich werde dies mit denselben Daten versuchen und sehen, wie viel schneller es ist. Vielen Dank für den Aufwand beim Erstellen und Bearbeiten der Funktion.
James
1

Ich denke, Ihre beste Wette ist es, eine kleine Routine mit ArcPy zu schreiben. Ich würde so etwas wie diesen Pseudocode erstellen:

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'

Ich bin nicht sicher, was Sie mit den Informationen tun möchten, aber ich nehme an, Sie könnten ein Feld in Ihrer Tabelle erstellen und es mit der Anzahl aufeinanderfolgender Punkte aktualisieren (wenn ja, fügen Sie das Feld zuerst hinzu).

Ich würde empfehlen, Feature-Layer zu erstellen (wie eine Datenbanktabellenansicht, jedoch für Features in Arc). Machen Sie zwei aus den Originaldaten und öffnen Sie einen Aktualisierungscursor auf dem ersten Satz, der Ihre Gesamtsortierung angibt (da ESRI keine vollständigen SQL-Abfragen berücksichtigt). Verwenden Sie den zweiten, um nach Position auszuwählen, und öffnen Sie einen Suchcursor auf dem resultierenden Auswahlsatz.

[BEARBEITEN NACH James ANFRAGE] Mit Model Builder ausarbeiten. Wenn Sie Model Builder noch nicht verwendet haben, müssen Sie nur mit der rechten Maustaste in arcToolbox klicken. Wählen Sie 'Toolbox hinzufügen'. Klicken Sie mit der rechten Maustaste auf die neue Toolbox und klicken Sie auf "Neu-> Modell". Wenn Sie ein neues Modellfenster haben, ziehen Sie die benötigten Werkzeuge und Daten per Drag & Drop in das Fenster und verknüpfen Sie sie visuell (mit dem kleinen Pfeilwerkzeug). Wenn Sie so weit wie möglich gekommen sind (Sie können Ihre Cursor hier nicht hinzufügen), verwenden Sie die Option im Dateimenü des Model Builder, um nach Python zu exportieren. Das bringt Sie den größten Teil des Weges dorthin. Es ist automatisch generierter Code, wird also ein bisschen böse, aber funktional sein. Verwenden Sie dann die Links in meiner Antwort oben, um die Cursor zu verstehen und hinzuzufügen.

Wenn Sie Python noch nicht kennen, haben Sie keine Angst davor, Code zu schreiben! Python ist eine sehr einfache Skriptsprache, mit der Sie schnell Ergebnisse erzielen können. Esri hat auch eine Anleitung dazu.

Wenn Sie mit Ihrem Code nicht weiterkommen, veröffentlichen Sie ihn in diesem Forum und bitten Sie um Hilfe. Es gibt viele Leute hier, die helfen können. Eine Warnung: Stellen Sie sicher, dass Sie die richtige Hilfe von ESRI verwenden. Sie haben die Python-API zwischen den Versionen 9.x und 10 (arcgisscripting bzw. arcpy) massiv geändert. Wenn Sie also ArcGIS 9.x verwenden, finden Sie die entsprechenden Links zu meinen!

MappaGnosis
quelle
Das sieht so aus, wie ich es gerne machen würde. Derzeit verwende ich jedoch keinen Code in ArcGIS, sondern wähle lediglich aus vordefinierten Optionen aus. Wie würde ich anfangen, die oben vorgeschlagene Methode zu verwenden / zu generieren? Ich möchte, dass die Ausgabe entweder eine neue Tabelle oder ein neues Feld ist, das der Tabelle mit der Anzahl aufeinanderfolgender Punkte hinzugefügt wird.
James
Siehe meine Bearbeitung zu meinem Hauptbeitrag.
MappaGnosis
JTB, melden Sie sich bitte mit demselben Konto an, mit dem Sie diese Frage gestellt haben, damit Sie Kommentare veröffentlichen können. (Um es einfacher zu machen, habe ich das JTB-Konto mit dem James-Konto zusammengeführt.)
whuber
Entschuldigung für die Kontoänderung. Ich habe die ursprüngliche Frage als neuer Benutzer gepostet, konnte dann aber nicht mehr in dieses Konto zurückkehren, da ich kein Passwort usw. hatte. Deshalb habe ich ein weiteres Konto JTB erstellt, das ich von nun an (hoffentlich) verwenden werde. Ich werde den Modellbauer wie von Sylvester vorgeschlagen starten, aber da ich ihn noch nie benutzt habe, kann es ein wenig dauern, bis ich mich daran gewöhnt habe und herausgefunden habe, welche Werkzeuge usw. ich verwenden soll. Ich werde mit Fortschritten und Fragen zurückkommen. Vielen Dank
James
Sylvester - Ich glaube, ich verstehe den Prozess, aber ich weiß nicht, mit welchen Tools ich wirklich anfangen muss. Entfernung? Puffer? In der Nähe von? Ich weiß nicht einmal, ob es für dieses Problem ein korrektes Tool gibt, das die oben genannten Aufgaben erfüllt. Ich bin sehr lernbegierig, aber sehr am Anfang.
James
1

Sie können den Modellbauer in ArcGIS verwenden, um aufeinanderfolgende ID-Werte zu finden. Ich habe mein Modell als Python-Skript exportiert. Der Code generiert einen neuen shp mit aufeinanderfolgenden ID-Werten. !ICH WÜRDE! ist das Basis-ID-Feld. Sie müssen den Pfad, den Namen und den ID-Feldnamen von point2.shp entsprechend Ihrem Fall aktualisieren.

# 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
quelle
Ich bin nicht sicher, was der obige Code bewirkt und wie er meine Anfrage beantwortet. Ich schätze die Hilfe, aber das ist alles ziemlich neu für mich, da ich noch nie Python oder Model Builder verwendet habe. Ändere ich das ID-Feld für jeden Prozess in die ID im Datensatz?
James
@ James, es hängt davon ab, ob sich Ihre ID ändert. Um den Code zu verwenden, kopieren Sie einfach den obigen Code, fügen Sie ihn ein und speichern Sie ihn in einer leeren txt-Datei. Aktualisieren Sie den Code so, dass er mit Ihrem point.shp-Namen und -Pfad übereinstimmt. Ändern Sie dann den ID-Namen im Codeabschnitt Feld berechnen (2) entsprechend Ihrem Feld point.shp ID. Speichern Sie die txt-Datei, und benennen Sie die Datei im Windows Explorer mit der Erweiterung .py um. Klicken Sie mit der rechten Maustaste auf die Datei und öffnen Sie sie mit python.exe, um sie zu testen.
Artwork21
Idealerweise kann dieses Skript in ein Skript-Tool eingebunden werden, das auch die Pufferung und Funktionsauswahl übernimmt. help.arcgis.com/de/arcgisdesktop/10.0/help/index.html#/…
artist2121