使用 gDistance rgeos 查找两个 SpatialPointsDataframes 之间的最近距离?
Finding the nearest distance between two SpatialPointsDataframes using gDistance rgeos?
我有两个独立但相关的问题。
首先,我想确定 subset_original_data.csv
文件中每个数据点到最近建筑工地 (construction_layer.csv
) 的距离。我正在尝试使用 gDistance()
函数来计算最近的邻居,但我也对其他想法持开放态度。
我想在我的 subset_original_data.csv
数据帧中附加这个距 construction_layer.csv
最近邻距离的新向量。也就是说,对于我的 subset_original_data.csv
数据框的每一行,我想要到最近的建筑工地的最小距离。
第二个目标是确定从每个 subset_original_data.csv
行到高速公路形状文件 (fwy.shp
) 的最近距离。我还想将这个新向量追加回 subset_original.csv
数据框。
我已经成功地将 construction_layer.csv
和 subset_original_data.csv
转换为 SpatialPointsDataFrame
。我还通过使用 readOGR()
函数读取形状文件将 fwy.shp
文件转换为 SpatialLinesDataFrame
。我不确定下一步该去哪里。非常感谢您的意见!
~ $ spacedSparking
这是我的数据:
construction_layer.csv, fwy.shp, subset_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 ??? ???
编辑:
解决方案:
如有疑问,请问一位 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!
bc@data$Nearest_Con<-min_constructionDistance
bc@data$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
可以使用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)
我有两个独立但相关的问题。
首先,我想确定 subset_original_data.csv
文件中每个数据点到最近建筑工地 (construction_layer.csv
) 的距离。我正在尝试使用 gDistance()
函数来计算最近的邻居,但我也对其他想法持开放态度。
我想在我的 subset_original_data.csv
数据帧中附加这个距 construction_layer.csv
最近邻距离的新向量。也就是说,对于我的 subset_original_data.csv
数据框的每一行,我想要到最近的建筑工地的最小距离。
第二个目标是确定从每个 subset_original_data.csv
行到高速公路形状文件 (fwy.shp
) 的最近距离。我还想将这个新向量追加回 subset_original.csv
数据框。
我已经成功地将 construction_layer.csv
和 subset_original_data.csv
转换为 SpatialPointsDataFrame
。我还通过使用 readOGR()
函数读取形状文件将 fwy.shp
文件转换为 SpatialLinesDataFrame
。我不确定下一步该去哪里。非常感谢您的意见!
~ $ spacedSparking
这是我的数据: construction_layer.csv, fwy.shp, subset_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 ??? ???
编辑:
解决方案: 如有疑问,请问一位 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!
bc@data$Nearest_Con<-min_constructionDistance
bc@data$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
可以使用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)