反向 map 投影 : how to get lat/lon coordinates from projected coordinates

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)

现在，我有一个场景，我需要对投影坐标进行一些计算，并将结果反向投影回纬度/经度坐标。例如，如何反向投影红点[1,0]？

有什么想法可以做到吗？

Answer 1

我不知道您是否需要首先使用 mapproject 进行投影。如果您可以使用 spTransform 来代替，那么这会变得更容易，因为您还可以使用 spTransform 来反转相同的过程。

假设您确实需要使用mapproject，我们仍然可以使用spTransform将点从投影坐标系转换为经纬度坐标，但需要更多的摆弄需要处理 mapproject 投影的非标准格式，即这些点被标准化为位于 -1 到 1 纬度和 -2 到 2 经度之间。在更标准的 map 投影中，纬度/经度以距离(通常是米)表示。

因此，首先我们可以使用 spTransform 找出将归一化的 map 投影坐标转换​​为实际距离所需的转换因子:

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 

现在我们可以使用这些引用将标准化 map 项目坐标转换为以米为单位的距离:

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