1

我正在使用纬度和经度数据，并尝试以图形方式显示斯德哥尔摩地图上每个点的指标（基于与该点的接近程度）。我对图像上等距的点更感兴趣，而不是在实际距离上等距：从这个意义上说，我理解赤道纬度点之间的距离比它们沿极圈的距离长，这可能是一个至关重要的问题。

我的目标是将地图划分为在 x 和 y 方向上大约 1 公里增量的网格。因此，我取了最小和最大纬度和经度，计算了它们与斯德哥尔摩中心的 x 和 y 距离，然后将经纬度的跨度除以 x 和 y 坐标的跨度（使用 geosphere）。我这样做是因为我希望这些点在绘制它们时等间距（否则，由于靠近赤道，与地图底部相比，顶部的 x 点之间的距离会更小）。

然后我在地图上绘制了这些点（使用 ggmap），并观察到 ​​y 方向上的点之间的距离比 x 方向上的距离更大。我想地图可以简单地以扭曲的方式绘制，但它似乎有点太扭曲了，难以置信。我怀疑我可能做错了什么，但找不到可能是什么。

下面的代码示例：

library("ggmap")
library("RgoogleMaps")
library("geosphere")

stockholm <- get_map("stockholm", zoom=11)
ggmap(stockholm)

places <- c('Tensta', 'Hanviken')

pos <- data.frame(Places = places, lat = NA, lon = NA, x = NA, y = NA)
reflatlon = getGeoCode('Stockholm, Sweden')

for(i in 1:length(places)) {
  latlon <- getGeoCode(paste0(places[i], ', Stockholm'))
  pos$lat[i] <- as.numeric(latlon[1])
  pos$lon[i] <- as.numeric(latlon[2])

  dist_y <- distGeo(c(latlon[1], reflatlon[2]), reflatlon) * sign(latlon[1] - reflatlon[1]) # same longitude
  dist_x <- distGeo(c(reflatlon[1], latlon[2]), reflatlon) * sign(latlon[2] - reflatlon[2]) # same latitude

  pos$x[i] <- dist_x
  pos$y[i] <- dist_y
}

deglatperm <- ( max(pos$lat) - min(pos$lat) ) / ( max(pos$y) - min(pos$y) ) # degrees latitude per metre
deglonperm <- ( max(pos$lon) - min(pos$lon) ) / ( max(pos$x) - min(pos$x) ) # degrees longitude per metre

seqlat <- seq(min(pos$lat), max(pos$lat), by = deglatperm*1000) # sequence with a point every ~1km
seqlon <- seq(min(pos$lon), max(pos$lon), by = deglonperm*1000) # sequence with a point every ~1km

seqlatlon <- expand.grid(seqlat, seqlon)
names(seqlatlon) <- c('lat', 'lon')

ggmap(stockholm) + geom_point(aes(x = lon, y = lat), data=seqlatlon)

输出图

从输出图中可以看出，与 x 方向相比，y 方向上的点之间的距离至少是 x 方向上的两倍。

总结一下：x 和 y 坐标是使用 geosphere 获得的。该地图是使用 ggmap 绘制的。

我对地圈做错了吗？还是纬度和经度图如此扭曲？当我打开谷歌地图并使用“测量距离”工具大约在顶部和底部以及左右点之间时，我得到的估计值为 16.3 和 16.9 公里，而我得到的地球圈值是 17 和 32 公里（x 和 y ） 分别。

如果有人能告诉我这里发生了什么，我将非常感激！

4

1 回答 1

2

我尽量避免在使用度数而不是距离的任何坐标系中工作。在美国，我经常使用我们的 State Plane 系统。瑞典似乎使用 RT 系统。一旦您将坐标从度数中移到与基准面的距离中，您就可以使用更传统的距离来构建网格。如果你愿意，你可以从那里把你的坐标放回度数。

我使用spTranform函数进行坐标转换，并使用空间参考指南来获取坐标系的参考代码。

library("ggmap")
library("RgoogleMaps")
library("geosphere")
library("sp")

stockholm <- get_map("stockholm", zoom=11)
ggmap(stockholm)

places <- c('Tensta', 'Hanviken')

pos <- data.frame(Places = places, lat = NA, lon = NA)
reflatlon <- getGeoCode('Stockholm, Sweden')

for(i in 1:length(places)) {
  latlon <- getGeoCode(paste0(places[i], ', Stockholm'))
  pos$lat[i] <- as.numeric(latlon[1])
  pos$lon[i] <- as.numeric(latlon[2])
}

p <- SpatialPoints(data.frame(pos$lon, pos$lat), proj4string = CRS("+init=epsg:4326"))
p <- spTransform(p, CRS("+init=epsg:3022"))

seqx <- seq(min(p@coords[,1]), max(p@coords[,1]), by = 1000)
seqy <- seq(min(p@coords[,2]), max(p@coords[,2]), by = 1000)

pgrid <- expand.grid(seqx, seqy)
pgrid <- SpatialPoints(pgrid, proj4string = CRS("+init=epsg:3022"))
pgrid <- spTransform(pgrid, CRS("+init=epsg:4326"))
pgrid <- data.frame(pgrid@coords)

names(pgrid) <- c('lon', 'lat')

ggmap(stockholm) + geom_point(aes(x = lon, y = lat), data=pgrid)
于 2016-04-08T01:28:07.540 回答