Identificatie van opeenvolgende punten binnen een opgegeven buffer

Ik heb een puntbestand, elk uur verplaatsingen van een dier en ik wil in staat zijn om een ​​buffer rond elk punt te plaatsen en het aantal opeenvolgende punten binnen de bufferzone te berekenen. Ik ben op zoek naar een methode die langs het puntbestand werkt, zoals een bewegend venster, dat alleen de volgende punten telt die zich binnen de bufferzone bevinden.

Bijvoorbeeld op punt 10 plaats ik een buffer van 1500 meter en ik wil weten of punt 11 zich in de buffer bevindt en zo ja, of punt 12 binnen de buffer ligt enzovoort. Ik wil niet weten of punt 52 zich in de bufferzone bevindt, tenzij alle vorige punten zich in de buffer bevinden. Ik wil ook niet weten of punten 9 of 8 enz. Zich in de buffer bevinden.

Ik heb een extra toolbox gevonden en geprobeerd, genaamd "moving window point analysis", die werkt als een bewegend venster op het puntbestand. Dit werkt goed, maar heel langzaam, en omvat alle punten die zich binnen de bufferzone bevinden, ook als dit geen opeenvolgende punten zijn. Ik kan geen manier vinden om alleen naar opeenvolgende punten te kijken.

Ik zou een methode willen die een uitvoertabel zal verstrekken aangezien ik veel gegevenspunten heb om op deze manier te bekijken.

Ik gebruik ArcGIS 10. Alle hulp die iemand kan bieden, wordt zeer op prijs gesteld.

8
Uw punten zijn waarschijnlijk ontstaan ​​als een lijst met (x, y, tijd) gegevens. Zou je openstaan ​​om deze gegevens (buiten ArcGIS) voor te verwerken om de gewenste informatie te verkrijgen?
toegevoegd de auteur whuber, de bron
Ah! Omdat u R al gebruikt, laten we dan R-gebaseerde oplossingen verkennen.
toegevoegd de auteur whuber, de bron
Er is een schuifraamfunctie binnen AdehabitatLT "sliwinltr" maar ik weet niet hoe ik het in dit geval moet gebruiken. Ik weet niet eens of het op deze manier kan worden gebruikt.
toegevoegd de auteur Core, de bron
Als dat het gemakkelijker maakt dan zeker. Ik ben ook de gegevens aan het verwerken met behulp van AdehabitatLT in R, die afstanden en peilingen enz. Berekent. Ik begrijp het proces dat hieronder door Sylvester wordt gesuggereerd, maar ik heb moeite om te weten waar ik moet beginnen, omdat ik niet echt zeker weet welke hulpmiddelen ik moet gebruiken.
toegevoegd de auteur Core, de bron

3 antwoord

Gegeven een lijst met puntlocaties (bij voorkeur in geprojecteerde coördinaten, zodat afstanden eenvoudig kunnen worden berekend), kan dit probleem worden opgelost met vijf eenvoudigere bewerkingen :

  1. Bereken punt-puntafstanden.

  2. Bepaal voor elk punt i, i = 1, 2, ... de indexen van die punten op afstanden kleiner dan de bufferradius (bijvoorbeeld 1500).

  3. Beperk die indexen tot i of groter.

  4. Behoud alleen de eerste opeenvolgende groep indexen zonder pauze.

  5. Voer het aantal van die groep uit.

In R komt elk van deze overeen met één bewerking. Als u deze reeks op elk punt wilt toepassen, is het handig om het meeste werk in te voegen in een functie die we definiëren , dus:

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

(Zie hieronder voor een efficiëntere versie van deze functie.)

Ik heb deze functie flexibel genoeg gemaakt om verschillende puntenlijsten ( xy ) en bufferafstanden ( r ) als parameters te accepteren.

Normaal gesproken zou u een bestand met puntlocaties lezen (en, indien nodig, sorteren op tijd). Hier, om dit in actie te tonen, zullen we willekeurig enkele voorbeeldgegevens genereren :

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

Figure

Hun typische spatiëring is 300 * Sqrt (2) = ongeveer 500. We doen de berekening door deze functie toe te passen op de punten in de array xy (en de resultaten vervolgens terug te halen) op naar xy , omdat dit een handig formaat zou zijn voor export naar een GIS):

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

Vervolgens analyseert u de matrix result in R of schrijft u deze naar een bestand en importeert u deze in andere software. Dit is het resultaat van de voorbeeldgegevens :

                        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

(Houd er rekening mee dat de tellingen de punten bevatten waarop ze zijn gebaseerd, zodat elke telling 1 of hoger moet zijn.)


Als je veel duizenden punten hebt, is deze methode te inefficiënt : het berekent veel te veel punt-tot-punt-afstanden die niet nodig zijn. Maar omdat we het werk in de functie doorsturen hebben ingekapseld, is de inefficiëntie eenvoudig te herstellen. Dit is een versie die beter werkt als er meer dan een paar honderd punten bij betrokken zijn:

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
}

Om dit te testen, creëerde ik willekeurige punten zoals eerder, maar varieerde ik twee parameters: n (het aantal punten) en hun standaardafwijking (hard-gecodeerd als 300 hierboven). De standaarddeviatie bepaalt het gemiddelde aantal punten binnen elke buffer ("gemiddelde" in de onderstaande tabel): hoe meer er zijn, hoe langer dit algoritme duurt om te worden uitgevoerd. (Met meer geavanceerde algoritmen zal de looptijd niet zo veel afhankelijk zijn van het aantal punten in elke buffer.) Dit zijn enkele tijdstippen:

 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
8
toegevoegd
Ik heb het antwoord bewerkt om een ​​oplossing te bieden met een sterk verbeterde schaalbaarheid. Het is nu in staat paden aan te pakken die miljoenen punten bevatten.
toegevoegd de auteur whuber, de bron
Ja. (Zie de opmerking op de regel waar n is gedefinieerd.) U kunt het berekenen vanuit de matrix xy als n <- dim (xy) [1] . Herhaal uw vorige opmerking: om deze code te laten werken, moet u deze uitvoeren van begin tot eind zonder iets over te slaan. Uiteindelijk vervangt u het maken van voorbeeldgegevens door uw eigen gegevens te lezen; ergo, u moet alle van dat blok code vervangen, wat betekent dat u vervangende items voor zowel n als xy . Dezelfde principes gelden voor alle andere wijzigingen die u wilt maken. (StackOverflow is een goede plaats om algemene R vragen te stellen.)
toegevoegd de auteur whuber, de bron
Sorry; er is een typfout: ik verwijder de externe ")". (Ik maakte een last-minute vereenvoudiging van deze code. Het is meermaals getest dan ik wil toegeven!)
toegevoegd de auteur whuber, de bron
Ik heb de originele code uitgevoerd met puntbestanden van 2000 vermeldingen die elk een paar uur duurden, omdat je zegt dat er veel te veel ongebruikte punt-naar-punt-locaties zijn die de verwerkingstijd verlengen. De bovenstaande bewerking ziet eruit als een mooie oplossing en ik zal dit op dezelfde gegevens proberen en zien hoe veel sneller het is. Bedankt voor de moeite met het produceren en bewerken van de functie.
toegevoegd de auteur Core, de bron
Dat is geweldig, het draait nu. Ik heb zojuist mijn gegevens in xy geladen voor eenvoud om te controleren of het werkt. Het kost wat tijd om te rennen zoals je al zei maar lijkt het allemaal goed gedaan te hebben. Ik zal er handmatig een paar controleren met mijn GIS-kaart, maar het ziet er tot nu toe goed uit. Bedankt dat je me hebt geholpen om dit uit te werken, heel graag wilt leren in zowel GIS als R en ik sta op de steile leercurve.
toegevoegd de auteur Core, de bron
n komt uit uw voorbeeldgegevens. Wat heb ik nodig om het te vervangen, het totale aantal gegevenspunten binnen xy?
toegevoegd de auteur Core, de bron
Ik probeerde uit te vinden welke haak extra was. Nu vindt het niet "n" en ik kan niet bepalen waar n ook vandaan komt
toegevoegd de auteur Core, de bron
Dit lijkt de perfecte oplossing. De code loopt echter niet van "z <- sapply (1: n), function (u) {forward (u, xy, radius)})" zoals het zegt: "onverwacht", "in" z: sappig (1: n), "Als ik de komma verwijder, zegt het: Fout: onverwachte 'functie' in 'z <- sapply (1: n) -functie' Hebt u enig idee waarom dit zo is?
toegevoegd de auteur Core, de bron

Ik denk dat je het beste kunt doen om een ​​beetje routine te schrijven met ArcPy. Ik zou zoiets als deze pseudo-code maken:

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'

Ik weet niet zeker wat u met de informatie wilt doen, maar ik veronderstel dat u een veld in uw tabel kunt maken en dit kunt bijwerken met het aantal opeenvolgende punten (als dat het geval is, voegt u eerst het veld toe).

Ik zou aanraden om feature-layers te maken (zoals een database-tabelweergave maar voor functies in Arc). Maak twee van de oorspronkelijke gegevens en open een update-cursor op de eerste set die uw algemene sortering specificeert (omdat ESRI geen volledige SQL-query's honoreert). Gebruik de tweede om selecteer op locatie vanaf en open een zoekcursor op de resulterende selectie-set.

[EDIT AS PER Jame's REQUEST] Rough it out using Model Builder. If you have not used Model Builder before, all you have to do is right-click in arcToolbox. Select 'Add Tool box'. The right-click on the new toolbox and click 'New->model'. Once you have a new model window, drag and drop the tools and data you need into the window and visually link them together (using the little arrow tool). When you have got as far as you can (you'll not be able to add your cursors here), use the option in the Model Builder's File Menu to export to Python. That will get you most of the way there. It is auto-generated code so will be a bit nasty but functional. Then use the links in my answer above to understand and add the cursors.

Als je nieuw bent bij Python, wees dan niet bang om code te schrijven! Python is een zeer eenvoudige scripttaal om snel resultaten te krijgen. Esri heeft ook een aantal begeleiding .

Als je vast komt te zitten met je code, plaats deze dan op dit forum en vraag om hulp. Er is VEEL mensen hier die kunnen helpen. Eén waarschuwing: zorg ervoor dat u de juiste hulp van ESRI gebruikt. Ze veranderden de Python API massaal tussen versies 9.x en 10 (respectievelijk arcgisscripting en arcpy). Dus als u ArcGIS 9.x gebruikt, zoek dan de equivalente links naar de mijne!

1
toegevoegd
Zie mijn bewerking van mijn hoofdpost.
toegevoegd de auteur Nick, de bron
JTB, meld u aan met hetzelfde account dat u hebt gebruikt om deze vraag te posten, zodat u opmerkingen kunt plaatsen. (Om het gemakkelijker te maken, heb ik het JTB-account samengevoegd met het James-account.)
toegevoegd de auteur whuber, de bron
Sylvester - Ik denk dat ik het proces begrijp, maar ik weet niet welke tool (s) ik echt nodig heb om mee te beginnen. Afstand? Buffer? In de buurt? Ik weet niet eens of er een juiste tool voor dit probleem is die zal doen wat hierboven is genoemd. Ik wil graag leren, maar sta helemaal aan het begin.
toegevoegd de auteur Core, de bron
Sorry voor de wijziging van het account. Ik postte de oorspronkelijke vraag als een nieuwe gebruiker, maar toen kon ik niet terug naar dat account omdat ik geen wachtwoord had enz. Daarom heb ik een ander account JTB gemaakt dat ik vanaf nu zal gebruiken (hopelijk). Ik zal de modelbouwer starten zoals gesuggereerd door Sylvester, maar heb hem nooit gebruikt voordat het me wat tijd kost om er aan te wennen en om uit te zoeken welke hulpmiddelen ik moet gebruiken. Ik kom terug met vorderingen en vragen. Bedankt
toegevoegd de auteur Core, de bron
Dit lijkt op wat ik zou willen doen. Ik gebruik momenteel echter geen code binnen ArcGIS, alleen maar uit vooraf gedefinieerde opties. Hoe zou ik de voorgestelde methode hierboven beginnen/gebruiken? Ik wil graag dat de uitvoer een nieuwe tabel is of een nieuw veld toegevoegd aan de tabel met het aantal opeenvolgende punten.
toegevoegd de auteur Core, de bron

U kunt modelbouwer in ArcGIS gebruiken om opeenvolgende ID-waarden te vinden. Ik heb mijn model geëxporteerd als een pythonscript. De code genereert een nieuwe shp met opeenvolgende ID-waarden. !ID KAART! is het veld voor de basis-ID. U moet het point2.shp-pad, de naam en de ID-veldnaam bijwerken om aan te passen aan uw geval.

# 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")
1
toegevoegd
@James, het hangt ervan af of je ID verandert. Om de code te gebruiken, kopieert en plakt u de bovenstaande code en slaat u deze op in een leeg txt-bestand. Werk de code bij zodat deze overeenkomt met uw point.shp-naam en -pad. Verander vervolgens de ID-naam in het gedeelte Codeerveld (2) code om overeen te komen met uw point.shp ID-veld. Sla het txt-bestand op en verwijs in Windows Explorer het bestand met de extensie .py. Klik met de rechtermuisknop op het bestand en open met python.exe om te testen.
toegevoegd de auteur artwork21, de bron
Idealiter kan dit script worden aangesloten op een scripttool, die ook de buffering en functiekeuze verwerkt. help.arcgis.com/en/arcgisdesktop/ 10.0/help/index.html #/& hellip;
toegevoegd de auteur artwork21, de bron
Ik weet niet zeker wat de bovenstaande code doet en hoe deze mijn vraag beantwoordt. Ik waardeer de hulp maar dit is allemaal vrij nieuw voor mij, ik heb nog nooit Python of modelbouwer gebruikt. Verander ik het veld ID voor elk proces dat wordt vermeld in de ID in de gegevensset?
toegevoegd de auteur Core, de bron