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)

第一个任务：已经完成！

第一步是使用over() 的rgeos 加入我的SpatialPoints 和SpatialPolygonsDataFrame。 SpatialPolygonsDataFrame是通过getData('GADM', country='ITA', level=3)的rgeos获得的。
对于第一个完成的任务，目标是将它们所属的有关City 和Region 的信息与每个 GPS 坐标相关联。
我能够获得的结果的一个例子是：

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 的坐标不在shapefile_roads 中Lines 的确切位置，因此，我应该需要某种半径参数。但是，我不确定。

您能否建议我在我的SpatialPoints 和从 OpenStreetMap 的road_shapefile 获得的SpatialLinesDataFrame 的属性之间执行这种空间连接的正确方法？

如果有不清楚的地方，请不要犹豫。

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 

Answer 2

您需要将多边形与点连接，而不是线。为此，您可以使用rgeos::gBuffer() 在您的线条周围创建一个缓冲区。小心，因为缓冲区将在您的 Lines 的坐标系中。在您的情况下可能是度数（wgs84）（验证它）。根据您的情况选择正确的距离（width）。

LinesBuffer <- rgeos::gBuffer(shapefile_roads, width = 0.01)

然后您将能够使用over 将点与“LinesBuffer”连接起来（如果它们在同一坐标系中）。

result <- over(my_spdf, LinesBuffer)