R Map Library – Converting NZMG or NZTM to Latitude/Longitude

coordinate systemr

I am a statistician with no background in GIS so apologies for the newbie question.

We routinely use R for a range of statistical and graphic analysis and are interested in using its map library to superimpose statistical graphics onto a map of New Zealand. Drawing a map of New Zealand is easy, as is adding cities. But as far as I can tell to project more data onto the maps I need latitude and longitude information.

I have a database of key locations of interest to us with grid references in NZMG or NZTM projections. I've found the technical instructions on converting to lat and long here, and conceivably could write code in R that applies this process. But I've also found various web calculators that will do this for a specific point, so there must be relatively accessible software that does this automatically (I need to do it for about 2500 locations).

So, one version of my question is: how do I easily convert, in R or Excel, a list of NZMG or NZTM projections into latitude and longitude.

An alternative version of the question is: if I don't really need to convert them, is there some way I don't understand that the map projection packages in R can automatically project these locations onto a map?

Apologies again if I have missed some fundamental point.

EDIT / ADDITION after @whuber's helpful link to an almost identical question

Ok, so I've installed the proj4 library but yes I am having difficulty in specifying my proj4string. If my understanding of this page is correct the NZTM system is basically a universal transverse mercator based on WGS-84. If I understand this map correctly I should by in zone 60 (or 59, as they both cover New Zealand? but neither work anyway). And if I understand proj4 properly the code below should be close. One of those "ifs" must be wrong…

1781304 E 5485722 N is a "NZTM" coordinate in my database that should be a suburb of Wellington, in the south of the north island.

library(proj4)
library(maps)
library(mapproj)
proj4string <- "+proj=utm +zone=60 +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs "
p <- project(c(1781304, 5485722), proj=proj4string, inverse=T) # should be wellington
map('nz')
points(p[1], p[2], cex=5, pch=19, col=2)
points(-p[1], p[2], cex=5, pch=19, col=2) # In case I've misunderstood the sign. Appears on map but in the sea

Best Answer

NZTM is not UTM zone 60 South. Its parameter values should be specified like this:

+proj=tmerc +lat_0=0.0 +lon_0=173.0 +k=0.9996 +x_0=1600000.0 +y_0=10000000.0 +datum=WGS84 +units=m

(it's not really WGS84, but NZGD2000 but close enough)

NZMG uses:

+proj=nzmg +lat_0=-41.0 +lon_0=173.0 +x_0=2510000.0 +y_0=6023150.0 +ellps=intl +units=m

NZMG is on NZGD 1949 which is not WGS84/NZGD 2000. That means a datum transformation. I believe you need to use cs2cs rather than +proj. There's some discussion here or here and it's definitely been discussed before. A possible transformation uses these parameters:

+towgs84=59.47,-5.04,187.44,0.47,-0.1,1.024,-4.5993

Related Question