2017-02-14 13 views
0

私は2つの別々のしかし関連する質問があります。gDistance rgeosを使用して2つのSpatialPointsDataframesの間の最近接距離を見つける?

まず、subset_original_data.csvファイル内のすべてのデータポイントについて、最も近い建設現場(construction_layer.csv)までの距離を決定したいと考えています。私は最も近い隣人を計算するためにgDistance()関数を使用しようとしていますが、私は他のアイデアも開いています。

私のsubset_original_data.csvデータフレームに、construction_layer.csvからの最近隣の距離のこの新しいベクトルを追加したいとします。つまり、私のsubset_original_data.csvデータフレームのすべての行について、最も近い建設現場との最小距離が必要です。

第2の目標は、各subset_original_data.csv行から高速道路形状ファイル(fwy.shp)までの最寄り距離を決定することです。この新しいベクトルをsubset_original.csvデータフレームに追加したいと思います。

construction_layer.csvsubset_original_data.csvSpatialPointsDataFrameに変換しました。また、readOGR()機能を持つ形状ファイルを読み込んで、fwy.shpファイルをSpatialLinesDataFrameに変換しました。私は次に行くべき場所が分からない。あなたのご意見は大変ありがとうございます!

〜$ spacedSparking

ここに私のデータです:ここに私のコードですが、 construction_layer.csvfwy.shpsubset_original_data.csv

#requiring necessary packages: 
library(rgeos) 
library(sp) 
library(rgdal) 

#reading in the files: 
mydata <- read.csv("subset_original_data.csv", header = T) 
con <- read.csv("construction_layer.csv", header = T) 
fwy <- readOGR(dsn = "fwy.shp") 

#for those who prefer not to download any files: 
data.lat <- c(45.53244, 45.53244, 45.53244, 45.53244, 45.53245, 45.53246) 
data.lon <- c(-122.7034, -122.7034, -122.7034, -122.7033, -122.7033, -122.7032) 
data.black.carbon <- c(187, 980, 466, 826, 637, 758) 
mydata <- data.frame(data.lat, data.lon, data.black.carbon) 

con.lat <- c(45.53287, 45.53293, 45.53299, 45.53259, 45.53263, 45.53263) 
con.lon <- c(-122.6972, -122.6963, -122.6952, -122.6929, -122.6918, -122.6918) 
con <- data.frame(con.lat, con.lon) 

#I am not sure how to include the `fwy.shp` in a similar way, 
#so don't worry about trying to solve that problem if you would prefer not to download the file. 

#convert each file to SpatialPoints or SpatialLines Dataframes: 
mydata.coords <- data.frame(lon = mydata[,2], lat = mydata[,1], data = mydata) 
mydata.sp <- sp::SpatialPointsDataFrame(mydata.coords, data = data.frame(BlackCarbon = mydata[,3])) #appending a vector containing air pollution data 

con.coords <- data.frame(lon = con[,2], lat = con[,1]) 
con.sp <- sp:SpatialPointsDataFrame(con.coords, data = con) 

str(fwy) #already a SpatialLinesDataFrame 

#Calculate the minimum distance (in meters) between each observation between mydata.sp and con.sp and between mydata.sp and fwy objects. 

#Create a new dataframe appending these two nearest distance vectors back to the original mydata file. 

#Desired output: 
head(mydata.appended) 
    LATITUDE LONGITUDE BC6. NEAREST_CON (m) NEAREST_FWY (m) 
1 45.53244 -122.7034 187 ???    ??? 
2 45.53244 -122.7034 980 ???    ??? 
3 45.53244 -122.7034 466 ???    ??? 
4 45.53244 -122.7033 826 ???    ??? 
5 45.53245 -122.7033 637 ???    ??? 
6 45.53246 -122.7032 758 ???    ??? 

はEDIT:

解決策: 疑問がある場合は、Rウィザードの友人にお尋ねください!彼は地図を作った。

library(rgeos) 
library(rgdal) 
library(leaflet) 
library(magrittr) 

#Define Projections 
wgs84<-CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0") 
utm10n<-CRS("+proj=utm +zone=10 +ellps=GRS80 +datum=NAD83 +units=m +no_defs +towgs84=0,0,0") 

#creating example black carbon data by hand: 
lat <- c(45.5324, 45.5325, 45.53159, 45.5321, 45.53103, 45.53123) 
lon <- c(-122.6972, -122.6963, -122.6951, -122.6919, -122.6878, -122.6908) 
BlackCarbon <- c(187, 980, 466, 826, 637, 758) 
bc.coords <- data.frame(lat, lon, BlackCarbon) 
bc<-SpatialPointsDataFrame(data.frame(x=lon,y =lat),data=data.frame(BlackCarbon),proj4string = wgs84) 

# Project into something - Decimal degrees are no fun to work with when measuring distance! 
bcProj<-spTransform(bc,utm10n) 

#creating example construction data layer: 
con.lat <- c(45.53287, 45.53293, 45.53299, 45.53259, 45.53263, 45.53263) 
con.lon <- c(-122.6972, -122.6963, -122.6952, -122.6929, -122.6918, -122.6910) 
con.coords <- data.frame(con.lat, con.lon) 
con<-SpatialPointsDataFrame(data.frame(x=con.lon,y =con.lat),data=data.frame(ID=1:6),proj4string = wgs84) 
conProj<-spTransform(con,utm10n) 

#All at once (black carbon points on top, construction on the y-axis) 
dist<-gDistance(bcProj,conProj,byid=T) 

min_constructionDistance<-apply(dist, 2, min) 

# make a new column in the WGS84 data, set it to the distance 
# The distance vector will stay in order, so just stick it on! 
[email protected]$Nearest_Con<-min_constructionDistance 
[email protected]$Near_ID<-as.vector(apply(dist, 2, function(x) which(x==min(x)))) 

#Map the original WGS84 data 
pop1<-paste0("<b>Distance</b>: ",round(bc$Nearest_Con,2),"<br><b>Near ID</b>: ",bc$Near_ID) 
pop2<-paste0("<b>ID</b>: ",con$ID) 
m<-leaflet()%>% 
    addTiles()%>% 
    addCircleMarkers(data=bc,radius=8,fillColor = 'red',fillOpacity=0.8,weight=1,color='black',popup=pop1)%>% 
    addCircleMarkers(data=con,radius=8,fillColor = 'blue',fillOpacity=0.8,weight=1,color='black',popup=pop2) 
m 
+0

は、一般的には、むしろそれを取り付けるよりも、あなたの質問に例のデータセットを作成するために、より適切だ - 私は他人のために話すことはできませんが、私は理解しやすいです未知のソース – SymbolixAU

+0

からファイルをダウンロードするに熱心ではありませんよ。私は、ポスト内にサンプルデータを含める方法を試してみる。 – spacedSparking

+0

また、 'sp'パッケージのmeuseのような、組み込みのデータセットを利用することもできます:' data( "meuse") ' – SymbolixAU

答えて

1

haversine distance関数を使用し、関数型プログラミングを使用して目的の結果を得ることができます。

library(geosphere) 
find_min_dist <- function(site, sites) { 
    min(distHaversine(site, sites)) 
} 
#X is the data id, split into a list so you can iterate through each site point 
data <- split(mydata[ , 3:2], mydata$X) 
sapply(data, find_min_dist, sites = con.coords) 
+0

私を 'geosphere'パッケージに紹介してくれてありがとう。私は、あなたが使用している 'mydata'オブジェクトが少し残っています。それは 'mydata < - read.csv( 'subset_original_data.csv')か' mydata < - data.frame(data.lat、data.lon、data.black.carbon) 'ですか? – spacedSparking

+0

また、 'mydata $ X'引数の代わりに何をすべきかわかりません。私はそれが 'BC6.'変数であると仮定していますが、それは私に奇妙な見た目の出力を与えます:' head(data1) -773 -374 -288 -272 -212 -127 174.4300 229.2947 146。6889 204.5449 146.9234 132.5356 ' – spacedSparking

+0

あなたのコードであなたのcsvを読むと、Xは行名です。各行をリストに分割し、建設現場からの距離をすべて計算できるように、各行にIDとIDだけが必要です。私もあなたのコードで形状ファイルを読むことができませんでしたが、形状ファイルからすべての緯度、経度のペアを取得したら同じ原理です – troh

関連する問題