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

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]？ 有什么想法吗？

- Mark Heckmann

---

你需要使用“mapproject”投影的原因是什么？如果你可以使用“spTransform”，那么这将变得更容易，因为你也可以使用“spTransform”来反转同样的过程。 - dww

@dww 没有什么特别的原因，我只是不知道如何使用 sp 来完成上述任务。因此，如果对于 sp 新手来说易于理解，那么 sp 解决方案也是可以接受的。 - Mark Heckmann

好的，我已经添加了一个答案，展示了如何在mapproject中完成它。如果您还需要使用spTransform完成整个过程，我也可以添加另一个答案。 - dww

Answer 1

我不知道你是否有必要使用mapproject进行投影。如果可以使用spTransform代替，则这将变得更加容易，因为您还可以使用spTransform来反转相同的过程。

假设您确实需要使用mapproject，我们仍然可以使用spTransform将点从投影坐标系转换为经纬度坐标，但需要更多的操作来处理mapproject投影的非标准格式，即将点规范化为介于-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 

现在我们可以使用这些参考资料将规范化的地图投影坐标转换为米制距离：

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

- dww

---

你真是个明星。我想从长远来看，习惯使用sp是值得的，但对于只有几个任务而言，学习曲线似乎太陡峭了... 非常感谢。 :) - Mark Heckmann

我刚刚意识到一个错误，现在已经纠正过来了。我做很多气候科学的工作，通常使用从零到360度的经度。当然，一般人使用的是-180度到+180度。因此，我改变了旋转方式，使用my.points$x[my.points$x > 180] = my.points$x[my.points$x > 180] - 360来处理标准经度。 - dww

Answer 2

如果没有现成的解决方案，你可以按照以下方式编写自己的函数：

```python def my_function(): # 在这里添加你的代码 ```

这个函数将根据你的需求进行定制。请注意保留 HTML 标签。

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

- Taeke

---

基本上，您建议为预期坐标使用查找表。虽然这是一个简单的想法，但它非常近似，导致反向转换略不精确。当然，人们可以有一个更精细的序列。也许这会奏效。对于这个简单实用的想法，点赞 :) - Mark Heckmann