7

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
2
  • It's nice to see you here, Peter. Your question is answered at gis.stackexchange.com/questions/18940/…. You may have follow-up questions, though, such as how precisely to specify the source and target coordinate systems as a "proj4string". If not, please allow us to close this question as a duplicate.
    – whuber
    Feb 15, 2012 at 21:32
  • Thanks, have added to the question - exactly that follow-up... Feb 15, 2012 at 23:11

2 Answers 2

6

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

2
  • Thanks, fantastic. That works. I only needed one of either NZTM or NZMG so I'm happy with the NZTM specification, which works great. Feb 16, 2012 at 3:17
  • Do I simply append the +towgs84 transformation to the +proj string? May 13, 2021 at 21:55
6

One of the most overlooked aspects of transforming coordinates from NZGD1949 to NZGD2000, including NZMG1949 to NZTM2000, is that it should use a transformation method. There are three methods:

  • 3 parameter similarity, with 5 metres nominal accuracy
  • 7 parameter similarity, with 4 metres nominal accuracy
  • Distortion grid, with 0.1 - 1.0 metres nominal accuracy

These are available to most GIS, including ArcGIS and PROJ-based. The distortion grid approach is clearly the best method, and can be used with a "nzgd2kgrid0005.gsb" grid shift file, available in PROJ's proj-datumgrid-1.5.zip (or newer). If building PROJ from source, this zip file needs to be unzipped to the data subdirectory before ./configure.

The following PROJ.4 string is the correct one for NZMG1949 (SRID=27200):

+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

For NZ coordinates, it's also a good idea to check them using the LINZ Online Converter, otherwise you may have points that are hundreds of meters off.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.