Преобразование NZMG (или NZTM) в широту / долготу для использования в библиотеке карт R
Я статистик без опыта работы в ГИС, поэтому прошу прощения за вопрос новичка.
Мы регулярно используем R для различных статистических и графических анализов и заинтересованы в использовании его библиотеки карт для наложения статистических графиков на карту Новой Зеландии. Нарисовать карту Новой Зеландии легко, как и добавить города. Но, насколько я могу судить, чтобы спроецировать больше данных на карту, мне нужна информация о широте и долготе.
У меня есть база данных ключевых мест, представляющих для нас интерес, с координатной сеткой в проекциях NZMG или NZTM. Я нашел технические инструкции по конвертации в широту и долготу здесь, и, возможно, мог бы написать код в R, который применяет этот процесс. Но я также нашел различные веб-калькуляторы, которые делают это для конкретной точки, поэтому должно существовать относительно доступное программное обеспечение, которое делает это автоматически (мне нужно сделать это для примерно 2500 мест).
Итак, одна из версий моего вопроса: как мне легко преобразовать в R или Excel список проекций NZMG или NZTM в широту и долготу.
Альтернативная версия вопроса: если мне не нужно конвертировать их, есть ли какой-то способ, который я не понимаю, чтобы пакеты картографических проекций в R могли автоматически проецировать эти местоположения на карту?
Еще раз прошу прощения, если я упустил какой-то фундаментальный момент.
EDIT / ДОБАВЛЕНИЕ после полезной ссылки @whuber' на почти идентичный вопрос
Итак, я установил библиотеку proj4, но у меня возникли трудности с определением строки proj4. Если я правильно понимаю эту страницу, система NZTM - это, по сути, универсальный поперечный меркатор, основанный на WGS-84. Если я правильно понимаю эту карту, я должен находиться в зоне 60 (или 59, поскольку они обе охватывают Новую Зеландию? но в любом случае ни то, ни другое не работает). И если я правильно понимаю proj4, код ниже должен быть близок. Одно из этих "если" должно быть неправильным...
1781304 E 5485722 N - это координаты "NZTM" в моей базе данных, которые должны быть пригородом Веллингтона, на юге Северного острова.
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
NZTM не является зоной UTM 60 Юг. Значения его параметров должны быть указаны следующим образом:
+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
(на самом деле это не WGS84, а NZGD2000, но достаточно близко)
NZMG использует:
+proj=nzmg +lat_0=-41.0 +lon_0=173.0 +x_0=2510000.0 +y_0=6023150.0 +ellps=intl +units=m
NZMG находится на NZGD 1949, который не является WGS84/NZGD 2000. Это означает преобразование точек отсчета. Я считаю, что вам нужно использовать cs2cs, а не +proj. Есть некоторое обсуждение здесь или здесь, и это определенно обсуждалось раньше. Возможное преобразование использует эти параметры:
+towgs84=59.47,-5.04,187.44,0.47,-0.1,1.024,-4.5993
Одним из наиболее упускаемых из виду аспектов преобразования координат из NZGD1949 в NZGD2000, включая NZMG1949 в NZTM2000, является использование метода преобразования. Существует три метода:
Они доступны для большинства ГИС, включая ArcGIS и PROJ.4. Подход distortion grid является лучшим методом, и может быть использован с "nzgd2kgrid0005.gsb" файлом сдвига сетки, доступным в PROJ.4's proj-datumgrid-1.5.zip. Если вы собираете PROJ.4 из исходников, этот zip-файл необходимо разархивировать в подкаталог
nadперед./configure.Следующая строка PROJ.4 является правильной для NZMG1949 (SRID=27200):
Для координат NZ всегда проверяйте их с помощью Онлайн-конвертера LINZ, иначе вы можете получить точки, отличающиеся на сотни метров.