我有一堆NZ地图网格坐标,我想将其转换为lat/long。基于这个问题,以下是我尝试过的。
library(sp)
options(digits = 11) # to display to greater d.p.
尝试1:
proj4string <- "+proj=nzmg +lat_0=-41.0 +lon_0=173.0 +x_0=2510000.0
+y_0=6023150.0 +ellps=intl +units=m"
p <- proj4::project(c(2373200, 5718800), proj = proj4string, inverse=T)
尝试2
dat <- data.frame(id = c(1), x = c(2373200) , y = c(5718800))
sp::coordinates(dat) = ~x+y
sp::proj4string(dat) = CRS('+init=epsg:27200')
data_wgs84 <- spTransform(dat, CRS('+init=epsg:4326'))
print(data_wgs84)
如果我通过linz坐标转换工具运行坐标,我会得到一个稍微不同的结果,这就是"真实"的结果。
Results:
171.30179199 -43.72743909 # attempt 1 - ~200m off linz
171.30190004, -43.72577765 # attempt 2 - a few meters off linz
171.30189464, -43.72576664 # linz
根据Mike T的回答,我应该使用"失真网格转换方法",他链接到"nzgd2kgrid0005.gsb网格偏移文件"。
我的问题:是否可以在不下载额外文件(nzgd2kgrid0005.gsb(的情况下使用R进行转换?我想与其他人共享我的代码,而不需要他们下载任何其他文件。
非常感谢任何建议。
事实证明,这非常简单,如果您安装了rgdal
包,则会包含所需的nzgd2kgrid0005.gsb
文件,并且不需要下载任何额外的文件。
您只需要使用完整的PROJ.4字符串,如Mike T的回答所述。
dat <- data.frame(id = c(1), x = c(2373200) , y = c(5718800))
sp::coordinates(dat) = ~x+y
proj4string <- "+proj=nzmg +lat_0=-41 +lon_0=173 +x_0=2510000 +y_0=6023150
+ellps=intl +datum=nzgd49 +units=m +towgs84=59.47,-5.04,187.44,0.47,-0.1,1.024,-4.5993
+nadgrids=nzgd2kgrid0005.gsb +no_defs"
sp::proj4string(dat) = sp::CRS(proj4string)
data_wgs84 <- sp::spTransform(dat, sp::CRS('+init=epsg:4326'))
as.data.frame(data_wgs84)
id x y
1 171.3018946 -43.72576664
这与LINZ坐标转换工具的输出相同。希望这能为其他人节省一点时间。