R - SpatialPoints（GPS 坐标）和 SpatialLinesDataFrame 之间的空间连接

Question

我正在从事一个结合数据科学和 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)

第一个任务：已经完成！

第一步是连接我SpatialPoints用SpatialPolygonsDataFrame用over()的rgeos。的SpatialPolygonsDataFrame是通过获得getData('GADM', country='ITA', level=3)的rgeos。
对于第一个完成的任务，目标是将每个 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 上创建的 shapefile 中，可从以下网址获得：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不在的确切位置Lines在shapefile_roads，所以，我应该需要某种半径参数。但是，我不确定。

你能给我建议在我的SpatialPoints和SpatialLinesDataFrame从road_shapefileOpenStreetMap获得的属性之间执行这种空间连接的正确方法吗？

如果有什么不是很清楚，请不要犹豫，问。

Answer 1

您的示例数据

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因为您使用的是角度（经纬度）坐标。

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如果有多个匹配项，则返回任意匹配项，而raster::extract将全部返回。两者都不漂亮，我会避免这种方法。图示如下：

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