R - SpatialPoints（GPS 坐标）和 SpatialLinesDataFrame 之间的空间连接
R - Spatial Join Between SpatialPoints (GPS coordinates) and SpatialLinesDataFrame
我正在从事一个结合了数据科学和 GIS 的大学项目。我们需要找到一个能够从海量 GPS 坐标数据集中获取额外信息的开源解决方案。显然，我不能使用任何具有每日请求限制的 API。
数据
在这里您可以找到教授提供给我们的数据集样本：
longitude <- c(10.86361, 10.96062, 10.93032, 10.93103, 10.93212)
latitude <- c(44.53355, 44.63234, 44.63470, 44.63634, 44.64559)
longlat <- data.frame(longitude, latitude)
ID <- seq.int(1, 10)
第一个任务：已经完成！
第一步是使用 rgeos 的 over() 将我的 SpatialPoints 与 SpatialPolygonsDataFrame 连接起来。 SpatialPolygonsDataFrame是通过rgeos的getData('GADM', country='ITA', level=3)获得的。
对于第一个完成的任务，objective 是将每个 GPS 坐标关联到它们所属的 City 和 Region 的信息。
我能够获得的结果示例是：
require(sp)
require(rgeos)
my_spdf <- SpatialPointsDataFrame(coords = longlat, data = ID, proj4string = CRS(" +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 "))
italy_administrative_boundaries_level3 <- getData('GADM', country='ITA', level=3)
result <- over(my_spdf, italy_administrative_boundaries_level3)[, c("NAME_0", "NAME_1", "NAME_2", "NAME_3")]
result$ID <- ID
print(result)
第二个任务：我的问题
现在事情变得棘手了，因为我需要关联更多更深入的信息，例如 road_name 和 road_type。
此信息包含在 OpenStreetMap 上创建的形状文件中，可在以下位置找到：download.geofabrik.de/europe/italy.html。
我在 R 中加载了 shapefile，获得了 SpatialLinesDataFrame:
require(rgdal)
shapefile_roads <- readOGR(dsn = "./road", layer = "roads")
然后，我天真地尝试应用与加入 SpatialPoints 和 SpatialPolygonsDataFrame:
相同的技术
result <- over(my_spdf, shapefile_roads)
很明显，结果就是NA。我想到的一个可能原因是 my_df 的坐标不在 shapefile_roads 中 Lines 的确切位置，因此，我应该需要某种半径参数。但是，我不太确定。
你能建议我在我的 SpatialPoints 和从 OpenStreetMap 的 road_shapefile 获得的 SpatialLinesDataFrame 的属性之间执行这种空间连接的正确方法吗？
如果有什么不是很清楚的地方，请不要犹豫。
您需要用点而不是线连接多边形。为此，您可以使用 rgeos::gBuffer() 在您的线条周围创建一个缓冲区。要小心，因为缓冲区将在您的线的坐标系中。在你的情况下可能是学位（wgs84）（验证它）。根据您的情况选择正确的距离 (width)。
LinesBuffer <- rgeos::gBuffer(shapefile_roads, width = 0.01)
然后您将能够使用 over 将点与 "LinesBuffer" 连接起来（如果它们在同一坐标系中）。
result <- over(my_spdf, LinesBuffer)
您的示例数据
library(raster)
longitude <- c(10.86361, 10.96062, 10.93032, 10.93103, 10.93212)
latitude <- c(44.53355, 44.63234, 44.63470, 44.63634, 44.64559)
longlat <- data.frame(longitude, latitude)
ID <- data.frame(ID=1:5)
ita_gadm3 <- getData('GADM', country='ITA', level=3)[, c("NAME_0", "NAME_1", "NAME_2", "NAME_3")]
#use `sp::over` or `raster::extract`
result <- extract(ita_gadm3, longlat)
部分道路：
road <- spLines(cbind(longitude+.1, latitude), cbind(longitude-.1, rev(latitude)), cbind(longitude-.1, latitude+1), crs=crs(ita_gadm3))
现在找到最近的路段。您可以使用 geosphere::dist2Line，因为您使用的是 angular (lon/lat) 坐标。
library(geosphere)
geosphere::dist2Line(longlat, road)
# distance lon lat ID
#[1,] 2498.825 10.83212 44.53355 2
#[2,] 5527.646 11.03032 44.63470 1
#[3,] 5524.227 10.86062 44.63634 2
#[4,] 5577.372 10.86062 44.63634 2
#[5,] 5756.113 10.86062 44.63634 2
注意变量 ID，它指回道路。问题是 dist2line 目前很慢，而且你有一个大数据集。
替代方法是将您的空间数据转换为适合意大利的平面坐标系并使用 gDistance。
library(rgeos)
library(rgeos)
sp <- SpatialPoints(longlat, proj4string=crs(ita_gadm3))
spita <- spTransform(sp, "+proj=tmerc +lat_0=0 +lon_0=15 +k=0.9996 +x_0=2520000 +y_0=0 +ellps=intl +units=m")
rdita <- spTransform(road, "+proj=tmerc +lat_0=0 +lon_0=15 +k=0.9996 +x_0=2520000 +y_0=0 +ellps=intl +units=m")
gd <- gDistance(rdita, spita, byid=TRUE)
a <- apply(gd, 1, which.min)
a
#1 2 3 4 5
#2 1 2 2 2
即点2离路1最近，其他点离路2最近。
您可能需要分批处理点或图块以避免获得太大的距离矩阵。
Sébastien 建议的缓冲解决方案原则上可行，但由于没有合适的缓冲区大小而变得非常复杂。一方面，点可能位于任何缓冲区之外，另一方面，它们可能与多个缓冲区重叠。如果您使用缓冲区，sp::over return 是一个任意匹配（如果有多个匹配），而 raster::extract 将 return 全部匹配。两者都不漂亮，我会避免这种方法。此处图示：
b <- buffer(road, width=.15, dissolve=F)
plot(b)
lines(road, col='red', lwd=2)
points(longlat, pch=20, col='blue')
extract(b, longlat)
# point.ID poly.ID
#1 1 1
#2 1 2
#3 2 2
#4 2 1
#5 3 2
#6 3 1
#7 4 2
#8 4 1
#9 5 1
#10 5 2
over(sp, b)
#1 2 3 4 5
#2 2 2 2 2
我正在从事一个结合了数据科学和 GIS 的大学项目。我们需要找到一个能够从海量 GPS 坐标数据集中获取额外信息的开源解决方案。显然，我不能使用任何具有每日请求限制的 API。
数据
在这里您可以找到教授提供给我们的数据集样本：
longitude <- c(10.86361, 10.96062, 10.93032, 10.93103, 10.93212)
latitude <- c(44.53355, 44.63234, 44.63470, 44.63634, 44.64559)
longlat <- data.frame(longitude, latitude)
ID <- seq.int(1, 10)
第一个任务：已经完成！
第一步是使用 rgeos 的 over() 将我的 SpatialPoints 与 SpatialPolygonsDataFrame 连接起来。 SpatialPolygonsDataFrame是通过rgeos的getData('GADM', country='ITA', level=3)获得的。
对于第一个完成的任务，objective 是将每个 GPS 坐标关联到它们所属的 City 和 Region 的信息。
我能够获得的结果示例是：
require(sp)
require(rgeos)
my_spdf <- SpatialPointsDataFrame(coords = longlat, data = ID, proj4string = CRS(" +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 "))
italy_administrative_boundaries_level3 <- getData('GADM', country='ITA', level=3)
result <- over(my_spdf, italy_administrative_boundaries_level3)[, c("NAME_0", "NAME_1", "NAME_2", "NAME_3")]
result$ID <- ID
print(result)
第二个任务：我的问题
现在事情变得棘手了，因为我需要关联更多更深入的信息，例如 road_name 和 road_type。
此信息包含在 OpenStreetMap 上创建的形状文件中，可在以下位置找到：download.geofabrik.de/europe/italy.html。
我在 R 中加载了 shapefile，获得了 SpatialLinesDataFrame:
require(rgdal)
shapefile_roads <- readOGR(dsn = "./road", layer = "roads")
然后，我天真地尝试应用与加入 SpatialPoints 和 SpatialPolygonsDataFrame:
result <- over(my_spdf, shapefile_roads)
很明显，结果就是NA。我想到的一个可能原因是 my_df 的坐标不在 shapefile_roads 中 Lines 的确切位置，因此，我应该需要某种半径参数。但是，我不太确定。
你能建议我在我的 SpatialPoints 和从 OpenStreetMap 的 road_shapefile 获得的 SpatialLinesDataFrame 的属性之间执行这种空间连接的正确方法吗？
如果有什么不是很清楚的地方，请不要犹豫。
您需要用点而不是线连接多边形。为此，您可以使用 rgeos::gBuffer() 在您的线条周围创建一个缓冲区。要小心，因为缓冲区将在您的线的坐标系中。在你的情况下可能是学位（wgs84）（验证它）。根据您的情况选择正确的距离 (width)。
LinesBuffer <- rgeos::gBuffer(shapefile_roads, width = 0.01)
然后您将能够使用 over 将点与 "LinesBuffer" 连接起来（如果它们在同一坐标系中）。
result <- over(my_spdf, LinesBuffer)
您的示例数据
library(raster)
longitude <- c(10.86361, 10.96062, 10.93032, 10.93103, 10.93212)
latitude <- c(44.53355, 44.63234, 44.63470, 44.63634, 44.64559)
longlat <- data.frame(longitude, latitude)
ID <- data.frame(ID=1:5)
ita_gadm3 <- getData('GADM', country='ITA', level=3)[, c("NAME_0", "NAME_1", "NAME_2", "NAME_3")]
#use `sp::over` or `raster::extract`
result <- extract(ita_gadm3, longlat)
部分道路：
road <- spLines(cbind(longitude+.1, latitude), cbind(longitude-.1, rev(latitude)), cbind(longitude-.1, latitude+1), crs=crs(ita_gadm3))
现在找到最近的路段。您可以使用 geosphere::dist2Line，因为您使用的是 angular (lon/lat) 坐标。
library(geosphere)
geosphere::dist2Line(longlat, road)
# distance lon lat ID
#[1,] 2498.825 10.83212 44.53355 2
#[2,] 5527.646 11.03032 44.63470 1
#[3,] 5524.227 10.86062 44.63634 2
#[4,] 5577.372 10.86062 44.63634 2
#[5,] 5756.113 10.86062 44.63634 2
注意变量 ID，它指回道路。问题是 dist2line 目前很慢，而且你有一个大数据集。
替代方法是将您的空间数据转换为适合意大利的平面坐标系并使用 gDistance。
library(rgeos)
library(rgeos)
sp <- SpatialPoints(longlat, proj4string=crs(ita_gadm3))
spita <- spTransform(sp, "+proj=tmerc +lat_0=0 +lon_0=15 +k=0.9996 +x_0=2520000 +y_0=0 +ellps=intl +units=m")
rdita <- spTransform(road, "+proj=tmerc +lat_0=0 +lon_0=15 +k=0.9996 +x_0=2520000 +y_0=0 +ellps=intl +units=m")
gd <- gDistance(rdita, spita, byid=TRUE)
a <- apply(gd, 1, which.min)
a
#1 2 3 4 5
#2 1 2 2 2
即点2离路1最近，其他点离路2最近。 您可能需要分批处理点或图块以避免获得太大的距离矩阵。
Sébastien 建议的缓冲解决方案原则上可行，但由于没有合适的缓冲区大小而变得非常复杂。一方面，点可能位于任何缓冲区之外，另一方面，它们可能与多个缓冲区重叠。如果您使用缓冲区，sp::over return 是一个任意匹配（如果有多个匹配），而 raster::extract 将 return 全部匹配。两者都不漂亮，我会避免这种方法。此处图示：
b <- buffer(road, width=.15, dissolve=F)
plot(b)
lines(road, col='red', lwd=2)
points(longlat, pch=20, col='blue')
extract(b, longlat)
# point.ID poly.ID
#1 1 1
#2 1 2
#3 2 2
#4 2 1
#5 3 2
#6 3 1
#7 4 2
#8 4 1
#9 5 1
#10 5 2
over(sp, b)
#1 2 3 4 5
#2 2 2 2 2