反向地图投影:如何从投影坐标中获取 lat/lon 坐标
reverse map-projection: how to get lat/lon coordinates from projected coordinates
我有一组 lat/lon 坐标,我可以使用它们进行投影,例如 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)
现在,我有一个场景需要对投影坐标进行一些计算并将结果反向投影回 lat/lon 坐标。比如我如何反投影红点[1,0]
?
知道如何做到这一点吗?
如果没有现成可用的东西,您可以按照以下行编写自己的函数:
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']))
}
我不知道你是否有理由需要使用 mapproject
在第一个地方进行投影。如果您可以改用 spTransform
,那么这会变得更容易,因为您还可以使用 spTransform
来反转相同的过程。
假设您确实需要使用 mapproject
,我们仍然可以使用 spTransform
将点从您的投影坐标系转换为 lat-long 坐标,但是需要更多的操作处理 mapproject
投影的 non-standard 格式,即这些点被标准化为位于 -1 到 1 纬度和 -2 到 2 经度之间。在更标准的地图投影中,lat/long 以距离(通常为米)表示。
因此,首先我们可以使用 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'))
并将这些重新投影到 lat/long 地理坐标中:
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
我有一组 lat/lon 坐标,我可以使用它们进行投影,例如 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)
现在,我有一个场景需要对投影坐标进行一些计算并将结果反向投影回 lat/lon 坐标。比如我如何反投影红点[1,0]
?
知道如何做到这一点吗?
如果没有现成可用的东西,您可以按照以下行编写自己的函数:
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']))
}
我不知道你是否有理由需要使用 mapproject
在第一个地方进行投影。如果您可以改用 spTransform
,那么这会变得更容易,因为您还可以使用 spTransform
来反转相同的过程。
假设您确实需要使用 mapproject
,我们仍然可以使用 spTransform
将点从您的投影坐标系转换为 lat-long 坐标,但是需要更多的操作处理 mapproject
投影的 non-standard 格式,即这些点被标准化为位于 -1 到 1 纬度和 -2 到 2 经度之间。在更标准的地图投影中,lat/long 以距离(通常为米)表示。
因此,首先我们可以使用 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'))
并将这些重新投影到 lat/long 地理坐标中:
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