1

緯度と経度のデータを使用して、ストックホルムのマップの各ポイントでメトリックをグラフィカルに表示しようとしています (そのポイントへの近さに基づいて)。私は、実際の距離ではなく、画像上で点が等間隔に配置されていることに関心があります。この意味で、赤道での緯度点間の距離は、極円に沿った距離よりも長いことを理解しています。質問にとって重要です。

私の目標は、マップを x 方向と y 方向の両方で約 1 km 刻みのグリッドに分割することでした。そのため、緯度と経度の最小値と最大値を取得し、ストックホルムの中心からの x と y の距離を計算し、緯度と経度の範囲を x と y 座標の範囲で割りました (地球圏を使用)。これを行ったのは、ポイントをプロットするときにポイントを等間隔にしたかったからです (そうしないと、赤道に近いため、マップの下部に比べて上部の x ポイント間の距離が短くなります)。

次に、これらの点を (ggmap を使用して) マップ上にプロットし、点間の距離が x 方向よりも y 方向の方が大きいことを観察しました。地図は単純に歪んで描かれているのかもしれませんが、少し歪んでいるとは思えません。何か間違ったことをしている可能性があると思いますが、それが何であるかを見つけることができません。

以下のコード例:

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 方向のポイント間の距離は少なくとも 2 倍です。

要約すると、x 座標と y 座標は地圏を使用して取得されます。マップはggmapを使用してプロットされます。

私は地球圏で何か間違ったことをしていますか? それとも、緯度と経度の地図が歪んでいますか? Google マップを開き、「距離を測定」ツールを使用して、上下左右の点の間で約 16.3 km と 16.9 km の推定値を取得しますが、地球圏で取得した値は 17 km と 32 km (x と y) です。 ） それぞれ。

誰かがここで何が起こっているのか教えてくれたら、とても感謝しています!

4

1 に答える 1

2

距離の代わりに度を使用する座標系での作業は避けようとしています。米国では、私はステート プレーン システムを常に使用しています。スウェーデンは 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 に答える