纬度和经度、x 和 y，以及绘制它们时的差异：geosphere 和 ggmap

Question

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

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

output plot

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

总而言之：x和y坐标是使用geosphere获得的。该地图是使用 ggmap 绘制的。

我在 geosphere 方面做错了吗？还是经纬度地图 SO 扭曲了？当我打开 Google 地图，并使用 "measure distance" 工具大约在顶部和底部以及左侧和右侧点之间时，我得到的估计值是 16.3 和 16.9 公里，而我使用 geosphere 得到的值是 17 和分别为 32 公里（x 和 y）。

如果有人能告诉我这是怎么回事，我将不胜感激！

Answer 1

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

我使用 spTranform function for coordinate conversions and use I the Spatial Reference 指南获取坐标系的参考代码。

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)

纬度和经度、x 和 y，以及绘制它们时的差异：geosphere 和 ggmap

Latitude and Longitude, x and y, and differences when plotting them: geosphere and ggmap

r

ggplot2

latitude-longitude

ggmap

geosphere