r-纬度和经度，x和y，以及绘制它们时的差异:地圈和ggmap

Question

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

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

我是不是做错了地球圈？或者纬度和经度SO的地图失真了吗？当我打开谷歌地图，使用"测量距离"工具大约在顶部和底部、左侧和右侧点之间时，我得到的估计值分别为16.3和16.9公里，而我使用geosphere得到的值分别为17和32公里（x和y）。

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

Answer 1

我尽量避免在任何使用度数而不是距离的坐标系中工作。在美国，我经常使用我们的国家飞机系统。瑞典似乎使用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)