Convert NZMG coordinates to lat/long

Question

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.

The linked answer is pretty clear: for best accuracy you need a file. Note that you have already done the conversion in R without the file, you just have to decide whether the reduced accuracy is good enough for you. — Hong Ooi, Oct 28 '19 at 02:58

score 1 · Accepted Answer · answered Nov 15 '19 at 01:45

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.

Convert NZMG coordinates to lat/long

1 Answers1