R：GIS 数据中点与线之间的距离（以米为单位）

Question

我有一堆点和线（可能还有多边形）的经度和纬度，我需要测量 R 中的距离。

大多数人建议我使用 rgeos 库中的 gDistance 函数。该函数需要对经度和纬度进行一系列转换，然后才能以米为单位输出结果。我已经浏览了许多教程，但仍然出现问题。希望您能帮我找出错误。

首先我们创建一些点和线

# Points
points <- data.frame(long = c(12.5633074637037,12.54638671875,12.6039819633632,12.54638671875,12.5668119504436,12.54638671875,12.5482921600342,12.5428380966187,12.5983709560864,12.5914064335047),
                     lat = c(55.6730208606487,55.6685371398926,55.6592116097919,55.6685371398926,55.6855954585358,55.6685371398926,55.7007255554199,55.6902847290039,55.663807868529,55.684380959963))

# Lines
lines <- data.frame(id = c("a", "b"))
lines$matrices <- list(matrix(c(12.5737695648244,12.5736645937496,12.5729168988113,12.5722100725459,12.5720446280546,55.6793201903946,55.6792991790095,55.6791495067552,55.6790112547884,55.6789788981105), ncol = 2),
                       matrix(c(12.5763890840661,12.57598090855,12.5759575726618,12.5757666379295,12.5757392134412,55.6799504343614,55.6797510847791,55.6797426062619,55.6796732345451,55.6796625397541), ncol = 2))

其次，使用 sp 库将数据转换为空间坐标。对于转换，我使用了丹麦的投影，因为这是我的数据所在的位置。我还将单位设置为米。 proj.4 是从https://epsg.io/23032 获取的。

# Libraries
library(sp)
library(rgeos)

# Converting to spatial coordinates
points$sp <- lapply(1:10, function(i) SpatialPoints(points[i,c("long", "lat")], proj4string=CRS("+proj=utm +zone=32 +ellps=intl +towgs84=-87,-98,-121,0,0,0,0 +units=m +no_defs")))
lines$sp <- lapply(lines$matrices, function(x) SpatialLines(list(Lines(Line(x), ID="a")), proj4string=CRS("+proj=utm +zone=32 +ellps=intl +towgs84=-87,-98,-121,0,0,0,0 +units=m +no_defs")))

为了进行完整性检查，我使用 leaflet 绘制了数据。不要介意丑陋的代码，我知道它可以使用空间数据框更智能地完成。中间的蓝色小线是我需要测量到的距离。

# Library
library(leaflet)

# Plotting to visualize
leaflet() %>% addTiles() %>%
  addLabelOnlyMarkers(data = points$sp[[1]], label = "1", labelOptions = labelOptions(noHide = T)) %>%
  addLabelOnlyMarkers(data = points$sp[[2]], label = "2", labelOptions = labelOptions(noHide = T)) %>%
  addLabelOnlyMarkers(data = points$sp[[3]], label = "3", labelOptions = labelOptions(noHide = T)) %>%
  addLabelOnlyMarkers(data = points$sp[[4]], label = "4", labelOptions = labelOptions(noHide = T)) %>%
  addLabelOnlyMarkers(data = points$sp[[5]], label = "5", labelOptions = labelOptions(noHide = T)) %>%
  addLabelOnlyMarkers(data = points$sp[[6]], label = "6", labelOptions = labelOptions(noHide = T)) %>%
  addLabelOnlyMarkers(data = points$sp[[7]], label = "7", labelOptions = labelOptions(noHide = T)) %>%
  addLabelOnlyMarkers(data = points$sp[[8]], label = "8", labelOptions = labelOptions(noHide = T)) %>%
  addLabelOnlyMarkers(data = points$sp[[9]], label = "9", labelOptions = labelOptions(noHide = T)) %>%
  addLabelOnlyMarkers(data = points$sp[[10]], label = "10", labelOptions = labelOptions(noHide = T)) %>%
  addPolylines(data = lines$sp[[1]]) %>%
  addPolylines(data = lines$sp[[2]])

由于数据看起来不错，我继续计算点和线之间的距离

# Calculating distances to the lines
## Distance to line 1
lapply(1:10, function(i) gDistance(points$sp[[i]], lines$sp[[1]]))

[[1]]
[1] 0.01057527

[[2]]
[1] 0.02770124

[[3]]
[1] 0.03629248

[[4]]
[1] 0.02770124

[[5]]
[1] 0.008435626

[[6]]
[1] 0.02770124

[[7]]
[1] 0.03220399

[[8]]
[1] 0.03131842

[[9]]
[1] 0.02908369

[[10]]
[1] 0.01834859

显然输出不是以米为单位，而是可能以度为单位。有人能指出哪里出了问题吗？

Answer 1

我发现如果我使用不同的投影，然后再使用 spTransform，我会得到想要的结果。我不确定，我理解这个过程，但它确实符合我的要求。

# Converting to spatial coordinates
points$sp <- lapply(1:10, function(i) SpatialPoints(points[i,c("long", "lat")], proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")))
points$sp <- lapply(1:10, function(i) spTransform(points$sp[[i]], CRS("+proj=utm +zone=32 +ellps=intl +towgs84=-87,-98,-121,0,0,0,0 +units=m +no_defs")))

lines$sp <- lapply(lines$matrices, function(x) SpatialLines(list(Lines(Line(x), ID="a")), proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")))
lines$sp <- lapply(1:2, function(i) spTransform(lines$sp[[i]], CRS("+proj=utm +zone=32 +ellps=intl +towgs84=-87,-98,-121,0,0,0,0 +units=m +no_defs")))

以米为单位的距离

lapply(1:10, function(i) gDistance(points$sp[[i]], lines$sp[[1]]))
[[1]]
[1] 861.6909

[[2]]
[1] 1989.797

[[3]]
[1] 2937.754

[[4]]
[1] 1989.797

[[5]]
[1] 807.0334

[[6]]
[1] 1989.797

[[7]]
[1] 2845.563

[[8]]
[1] 2227.443

[[9]]
[1] 2319.78

[[10]]
[1] 1244.597