Search code examples
rcoordinate-transformationproj

Convert NZMG coordinates to lat/long

I have a bunch of NZ Map Grid coordinates, which I want convert to lat/long. Based on this question, here is what I tried.

library(sp)
options(digits = 11) # to display to greater d.p.

Attempt 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)

Attempt 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)

If I run my coordinates through the linz coordinate conversion tool I get a slightly different result, which is the "true" result. 

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

Based on Mike T's answer I should be using a "distortion grid transformation method" and he links to a "nzgd2kgrid0005.gsb grid shift file".

My Question: Is it possible to do this conversion using R without downloading additional files (nzgd2kgrid0005.gsb)? I want to share my code with others without them having to download any additional files.

Any advice much appreciated.

Solution

  • Turns out it is pretty simple, if you have the rgdal package installed, the required nzgd2kgrid0005.gsb file is included and you don't need to download anything extra.

    You just need to use the full PROJ.4 string as outlined in Mike T's answer.

    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

    Which is the same as the output from the LINZ coordinate conversion tool. Hopefully this saves someone else a bit of time.