reverse map-projection：投影された座標から緯度/経度座標を取得する方法

Question

私は、例えばMollweide投影を使用して投影できる緯度/経度座標のセットを持っています。今

library(mapproj) 
set.seed(0) 
n <- 100 
s <- data.frame(lon = rnorm(n, 0, 60), 
       lat = rnorm(n, 0, 40)) 
p <- mapproject(s$lon, s$lat, proj="mollweide", par=NULL, 
       orientation=c(90,200,0))     

# plot projected coors 
plot(p$x, p$y, type="n", asp=1/1, bty="n") 
map.grid(c(-180, 180, -90, 90), nx=20, ny=20, 
     font=1, col=grey(.6), labels=F) 
points(p$x, p$y, pch="x", cex = .8) 

# a point to reverse project 
points(1,0, pch=16, col="red", cex=2) 

、私は私が投影座標上のいくつかの計算を行うと、バック緯度/経度COORDSにプロジェクトに結果を逆にする必要があるシナリオを持っています。たとえば、赤い点[1,0]をどのように反転させることができますか？

どうすればいいのでしょうか？

Answer 1

はあなたがmapprojectにを使用する必要がある理由があるかどうかはわかりません最初の場所にプロジェクト。 spTransformを代わりに使用できる場合は、spTransformを使用して同じプロセスを逆転することもできるため、これは簡単になります。

我々はまだあなたの予想からのポイントは、緯度・経度座標に座標系を変換するために、spTransformを使用することができ、あなたがmapprojectを使用する必要がないと仮定すると、もう少しあいはmapprojectの非標準フォーマットに対処するために必要とされます突起ポイントは-1から1の緯度と-2から2の経度の間にあるように正規化されます。より標準的な地図投影では、緯度/経度は距離（通常はメートル）で表されます。

ので、最初私たちは、正規化mapprojectは、実際の距離座標に変換する必要があり、変換係数を見つけるためにspTransformを使用することができます。

library(rgdal) 
my.points = data.frame(x=c(0,180),y=c(90,0)) 
my.points = SpatialPoints(my.points, CRS('+proj=longlat')) 
my.points = spTransform(my.points, CRS('+proj=moll')) 
# SpatialPoints: 
#    x  y 
# [1,]  0 9020048 
# [2,] 18040096  0 
# Coordinate Reference System (CRS) arguments: +proj=moll +ellps=WGS84 

今、私たちは、正規化mapprojectから変換するために、これらの参照を使用することができますが、距離に座標

my.points = data.frame(x=p$x * 18040096/2 , y=p$y * 9020048) 
my.points = SpatialPoints(my.points, CRS('+proj=moll')) 

そして、緯度/経度の地理座標にこれらを再投影：

my.points = as.data.frame(spTransform(my.points, CRS('+proj=longlat'))) 

メートルで

最後に、mapprojectで実行された回転を元に戻すには、これらの点を経度ごとに回転させる必要があります。

my.points$x = my.points$x + 200 
my.points$x[my.points$x > 180] = my.points$x[my.points$x > 180] - 360 

そして、それが働いていることを確認することができます：あなたは拳の場所に投影する `mapproject`を使用する必要があるいくつかの理由があります

を

head(my.points) 
#   x   y 
# 1 75.77725 31.274368 
# 2 -19.57400 -31.071065 
# 3 79.78795 -24.639597 
# 4 76.34576 1.863212 
# 5 24.87848 -45.215432 
# 6 -92.39700 23.068752 

head(s) 
#   lon  lat 
# 1 75.77726 31.274367 
# 2 -19.57400 -31.071065 
# 3 79.78796 -24.639596 
# 4 76.34576 1.863212 
# 5 24.87849 -45.215431 
# 6 -92.39700 23.068751 

Answer 2

ボックスのうち何が利用できない場合、あなたはの線に沿って、あなた自身の関数を書くことができます：

ref_project <- function(x, y) { 
    long <- tibble(
     long = seq(-180, 180, 1), 
     x = mapproject(long, rep(0, length(long)), projection = 'mollweide', orientation = c(90, 200, 0))$x 
    ) 
    lat <- tibble(
     lat = seq(-90, 90, 1), 
     x = mapproject(rep(0, length(lat)), lat, projection = 'mollweide', orientation = c(90, 200, 0))$y 
    ) 

    return(c(long[which(abs(long$x - x) == min(abs(long$x - x))), 'long'], 
      lat[which(abs(lat$x - y) == min(abs(lat$x - y))), 'lat'])) 
}