2017-03-01 9 views
0

Ich fragte this question vorher, aber habe keine Antwort erhalten, also werde ich versuchen, eine bessere Arbeit diesmal zu machen!Erstelle Puffer und zähle Punkte in R

Ich möchte die räumliche Dichte von Tankstellenpunkten mit R analysieren. Ich muss einen Puffer (sagen wir 1.000 m) um die Tankstellen herum erstellen und die Anzahl der Tankstellen im Puffer zählen. Ich muss dann mit Pufferdistanzen herumspielen, um zu sehen, was ein vernünftiger Puffer ist, um etwas Interessantes zu sehen. Ich werde die gesamte Form Datei nicht veröffentlichen, weil es ziemlich chaotisch, aber das ist, was die Daten wie folgt aussehen:

all <- readShapePoints("sbc_gas.shp") 
all.df <- as(all, "data.frame") 
head(all) 

OBJECTID Fuellocati    Name Latitude  Longitude  
     1  34828  WORLD OIL #104 34.44190  -119.8304  
     2  48734 STOP AND SHOP GAS 34.41962  -119.6768  
     3  51276 EL RANCHERO MARKET 34.41911  -119.7162  
     4  52882 EDUCATED CAR WASH 34.44017  -119.7439  
     5  74038   CIRCLE K 34.63925  -120.4406  
     6  103685 7-ELEVEN #23855 34.40506  -119.5296  

ich in der Lage war, einen Puffer um die Punkte mit dem folgenden Code zu erstellen, aber jetzt, wie ich Zählen Sie die Anzahl der Punkte im Puffer?

require(sp) 
require(rdgal) 
require(geosphere) 

coordinates(all) <- c("Longitude", "Latitude") 
pc <- spTransform(all, CRS("+init=epsg:3347")) 
distInMeters <- 1000 
pc100km <- gBuffer(pc, width=100*distInMeters, byid=TRUE) 
# Add data, and write to shapefile 
pc100km <- SpatialPolygonsDataFrame(pc100km, [email protected]) 
writeOGR(pc100km, "pc100km", "pc100km", driver="ESRI Shapefile") 

plot(pc100km) 

enter image description here

Ich bin offen für andere Möglichkeiten, um dies zu realisieren.

+1

Können Sie weitere Informationen über geben Erhalten Sie eine Fehlermeldung? Können Sie in den Debug-Modus wechseln und prüfen, was die Variablen enthalten? –

+0

Fehlermeldung zu Frage hinzugefügt. – JAG2024

+1

Für den Anfang sieht es aus, als wäre 'distm' nicht definiert. Meinst du stattdessen 'distInMeters'? –

Antwort

1

kam mit einer (einfachen) Lösung für meine eigene Frage geosphere mit bis: „kann auch nicht an der Arbeit“

cbind(coordinates(all.df), X=rowSums(distm (coordinates(all.df)[,1:2], fun = distHaversine)/1000 <= 10)) # number of points within distance 10 km 
    Longitude Latitude X 
0 -119.8304 34.44190 25 
1 -119.6768 34.41962 29 
2 -119.7162 34.41911 34 
3 -119.7439 34.44017 39 
4 -120.4406 34.63925 13 
5 -119.5296 34.40506 7 
6 -120.4198 34.93860 26 
7 -119.8221 34.43598 30