纬度和经度，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坐标是用地圈获得的。该地图使用ggmap绘制。

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

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

Answer 1

我尽量避免在任何坐标系中使用度数而不是距离。在美国，我不断使用我们的国家飞机系统。看来瑞典使用RT系统。一旦你的坐标超出度数并与基准距离，那么你就可以使用更传统的距离来建立你的网格。如果你喜欢，你可以从那里把你的坐标恢复到度数。

我使用spTranform函数进行坐标转换，并使用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)