NZMG (of NZTM) omzetten in breedtegraad/lengtegraad voor gebruik met R-kaartbibliotheek

Ik ben een statisticus zonder achtergrond in GIS, dus excuses voor de newbie-vraag.

We gebruiken routinematig R voor een reeks statistische en grafische analyses en zijn geïnteresseerd in het gebruik van de kaartenbibliotheek om statistische afbeeldingen op een kaart van Nieuw-Zeeland te plaatsen. Een kaart van Nieuw-Zeeland tekenen is eenvoudig, net als het toevoegen van steden. Maar voor zover ik kan zien om meer gegevens op de kaarten te projecteren, heb ik informatie over de breedte- en lengtegraad nodig.

Ik heb een database met voor ons belangrijke belangrijke locaties met rasterreferenties in NZMG- of NZTM-projecties. Ik heb de technische instructies gevonden voor het converteren naar lat en long hier , en zou mogelijk code kunnen schrijven in R die dit proces toepast. Maar ik heb ook verschillende webcalculators gevonden die dit voor een specifiek punt doen, dus er moet relatief toegankelijke software zijn die dit automatisch doet (ik moet het doen voor ongeveer 2500 locaties).

Dus, één versie van mijn vraag is: hoe converteer ik eenvoudig, in R of Excel, een lijst met NZMG- of NZTM-projecties naar breedtegraad en lengtegraad.

Een alternatieve versie van de vraag is: als ik ze niet echt hoef te converteren, begrijp ik dan op een of andere manier niet dat de kaartprojectiepakketten in R deze locaties automatisch op een kaart kunnen projecteren?

Excuses als ik een fundamenteel punt heb gemist.

EDIT/ADDITION na de handige link van @ whuber naar een bijna identieke vraag

Ok, dus ik heb de proj4-bibliotheek geïnstalleerd, maar ja, ik heb problemen met het specificeren van mijn proj4-ring. Als mijn kennis van deze pagina klopt, is het NZTM-systeem correct is in feite een universele dwarse mercator op basis van WGS-84. Als ik begrijp deze kaart correct zou ik het moeten doen in zone 60 (of 59 als ze behandelen allebei Nieuw-Zeeland, maar geen van beide werken toch). En als ik proj4 goed begrijp moet de onderstaande code in de buurt zijn. Een van die "ifs" moet fout zijn ...

1781304 E 5485722 N is een "NZTM" -coördinaat in mijn database die een buitenwijk van Wellington zou moeten zijn, in het zuiden van het noordelijke eiland.

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
5
Het is fijn je hier te zien, Peter. Uw vraag wordt beantwoord op gis.stackexchange.com/questions/18940/… . U kunt echter vervolgvragen hebben, zoals hoe u de bron- en doelcoördinatenstelsels precies kunt specificeren als een "proj4string". Als dit niet het geval is, staat u ons toe deze vraag als een duplicaat te sluiten.
toegevoegd de auteur whuber, de bron
Bedankt, hebben toegevoegd aan de vraag - precies die follow-up ...
toegevoegd de auteur James Williams, de bron

2 antwoord

NZTM is geen UTM-zone 60 Zuid. De parameterwaarden moeten als volgt worden gespecificeerd:

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

(het is niet echt WGS84, maar NZGD2000 maar dichtbij genoeg)

NZMG gebruikt:

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

NZMG is op NZGD 1949 wat geen WGS84/NZGD 2000 is. Dat betekent een datumtransformatie. Ik geloof dat je cs2cs moet gebruiken in plaats van + proj. Er is wat discussie hier of hier en het is zeker eerder besproken. Een mogelijke transformatie gebruikt deze parameters:

+ Towgs84 = 59.47, -5.04,187.44,0.47, -0.1,1.024, -4,5993

5
toegevoegd
Bedankt, fantastisch. Dat werkt. Ik had alleen een van NZTM of NZMG nodig, dus ik ben blij met de NZTM-specificatie, die prima werkt.
toegevoegd de auteur James Williams, de bron

Een van de meest over het hoofd geziene aspecten van transformerende coördinaten van NZGD1949 naar NZGD2000 , inclusief NZMG1949 naar NZTM2000, is dat er een transformatiemethode moet worden gebruikt. Er zijn drie methoden:

  • 3 parameterovereenkomst, met 5 meter nominale nauwkeurigheid
  • 7 parameterovereenkomst, met 4 meter nominale nauwkeurigheid
  • Vervormingsraster met een nominale nauwkeurigheid van 0,1 - 1,0 meter

Deze zijn beschikbaar voor de meeste GIS, inclusief op ArcGIS en PROJ.4 gebaseerd. Het distorsieraster is duidelijk de beste methode en kan worden gebruikt met een "nzgd2kgrid0005.gsb" roosterverschuivingsbestand, beschikbaar in PROJ.4's proj-datumgrid-1.5.zip . Als PROJ.4 wordt gebouwd vanaf de bron, moet dit zipbestand worden uitgepakt in de submap nad vóór ./ configure .

De volgende PROJ.4-reeks is de juiste voor 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

Voor NZ-coördinaten, controleer altijd met behulp van de LINZ Online Converter , anders heb je misschien punten die zijn honderden meters uit.

4
toegevoegd