緯度和經度，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)