反向地图投影：如何从投影坐标获取纬度/经度坐标

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来颠倒相同的过程。

假设你需要使用mapproject，我们仍然可以使用spTransform来分，从您的投影坐标系统转换为纬度和经度坐标，却多了几分摆弄来处理与mapproject的非标准格式预测Ie这些点被标准化为-1至1纬度和-2至2经度。在更标准的地图投影中，纬度/长度用距离表示（通常为米）。

所以，首先我们可以用spTransform找出我们需要的标准化mapproject转换的转换系数坐标转换成实际距离：

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 

并让检查它的工作：

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'])) 
}