反向地图投影:如何从投影坐标中获取纬度/经度坐标
Posted
技术标签:
【中文标题】反向地图投影:如何从投影坐标中获取纬度/经度坐标【英文标题】:reverse map-projection: how to get lat/lon coordinates from projected coordinates 【发布时间】:2017-06-29 16:26:12 【问题描述】:我有一组纬度/经度坐标,我可以使用例如 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]
?
有什么想法可以做到吗?
【问题讨论】:
是否有某些原因需要使用mapproject
进行第一处投影。如果您可以改用spTransform
,那么这会变得更容易,因为您还可以使用 spTransform 来反转相同的过程。
@dww 没有原因,我只是不知道如何使用sp
进行上述操作。因此,如果sp
newbees 可以理解,也欢迎sp
解决方案:)
好的,我添加了一个已经显示如何为 mapproject 执行此操作的答案。我还可以在没有 mapproject 的情况下在 spTransform
中添加另一个答案 - 如果您也需要这个
【参考方案1】:
如果没有现成可用的东西,您可以编写自己的函数,如下所示:
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']))
【讨论】:
基本上,您建议使用查找表来查找投影坐标。虽然是一个简单的想法,但它非常近似,导致反向转换稍微不精确。当然,可以有一个查找器粒度序列。也许这会奏效。这个简单的动手想法的道具:)【参考方案2】:我不知道您是否有理由需要使用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
现在我们可以使用这些引用将归一化的 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
【讨论】:
你是明星。我认为从长远来看,习惯 sp 是值得的,但是对于一些任务来说,学习曲线似乎太陡峭了……非常感谢。 :) 我刚刚意识到我现在已经纠正了一个错误。我做了很多气候科学,我们通常使用从零到 360 的经度。当然,通常人们使用 -180 到 +180。所以我改变了旋转以使用my.points$x[my.points$x > 180] = my.points$x[my.points$x > 180] - 360
来处理标准经度以上是关于反向地图投影:如何从投影坐标中获取纬度/经度坐标的主要内容,如果未能解决你的问题,请参考以下文章