2017-03-01 3 views
0

J'ai déjà demandé this question auparavant, mais je n'ai pas reçu de réponse, alors je vais essayer de faire un meilleur travail cette fois-ci! Je veux analyser la densité spatiale des points de station-service en utilisant R. Je dois créer un tampon (disons 1 000 m) autour des stations-service et compter le nombre de stations-service dans le tampon. Je devrai alors jouer avec les distances de tampon pour voir ce qui est un tampon raisonnable pour voir quelque chose d'intéressant. Je ne vais pas poster le dossier complet de forme car il est assez désordonné, mais c'est ce que les données ressemblent à:Créer un tampon et compter les points dans R

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  

j'ai pu créer un tampon autour des points avec le code suivant, mais maintenant comment puis-je compte le nombre de points dans le tampon?

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

Je suis ouvert à d'autres façons d'aller à ce sujet.

+1

Pouvez-vous donner plus d'informations sur Obtenez-vous une sortie d'erreur? Pouvez-vous entrer en mode débogage et inspecter ce que les variables tiennent? –

+0

Ajout de messages d'erreur à interroger. – JAG2024

+1

Pour commencer, il semble que 'distm' ne soit pas défini. Vouliez-vous y placer 'distInMeters' à la place? –

Répondre

1

est venu avec une solution (simple) pour ma propre question en utilisant geosphere: « ne peut pas obtenir soit au travail »

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