2017-03-01 30 views
0

我之前詢問過this question,但沒有得到回覆,所以我會盡力做好這一次的工作!在R中創建緩衝區和計數點

我想用R分析加油站點的空間密度。我需要在加油站周圍創建一個緩衝區(比方說1000米)並計算緩衝區內的加油站數量。然後,我需要充分利用緩衝距離來查看什麼是合理的緩衝區,以查看有趣的內容。我不會發布整體形狀文件,因爲它是相當凌亂,但這是該數據的樣子:

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  

我能夠創建點周圍的緩衝用下面的代碼,但現在我該怎麼辦計算緩衝區內的點數?

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

我打開其他的方式去了解這一點。

+1

你能給更多信息你是否收到任何錯誤輸出?你能進入調試模式並檢查變量持有的是什麼嗎? –

+0

添加錯誤消息以提問。 – JAG2024

+1

對於初學者來說,它看起來像'distm'沒有被定義。你的意思是把'distInMeters'放在那裏嗎? –

回答

1

想出了我自己的問題一個(簡單的)解決方案使用geosphere:「不能讓任何工作」

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