Ошибка при перепроецировании пространственных точек с использованием spTransform в rgdal R - PullRequest
Новая
       34

Ошибка при перепроецировании пространственных точек с использованием spTransform в rgdal R

3 голосов
/ 02 апреля 2012

G'day,
У меня есть большое количество координат lon / lat, которые находятся в CRS Australian Geodetic Datum 66/84 (AGD66 для краткости). Я хочу изменить эти координаты с AGD66 на WGS84, потому что между ними есть разница примерно в 200 м, и у меня есть другие координаты и слои в WGS84. Я пытался сделать это с помощью: 

lon        lat
147.1428   -43.49083

library(rgdal)
pts<-read.table(file.choose(),header=TRUE,sep=',')  
# I first project the pts in their original CRS
pts66<-project(cbind(pts$lon,pts$lat), "+init=epsg:4202")
# Then to transform it to WGS84
pts84 = spTransform(pts66,CRS("+init=epsg:3033"))

Error in function (classes, fdef, mtable)  : 
unable to find an inherited method for function "spTransform", for signature "matrix", "CRS"

Кто-нибудь знает, почему я получаю эту ошибку, или есть какие-либо предложения о том, как я могу изменить эти координаты с AGD66 на WGS84? Заранее спасибо за помощь.

Приветствия
Адам

1 Ответ

6 голосов
/ 02 апреля 2012

Я удалил часть неправильного ответа .

Функция project () не может выполнить преобразование базовых данных, поэтому у вас может возникнуть проблема, и я думаю, что у вас что-то не так.

Проблема в том, что вы можете только проецировать () из / в longlat на WGS84, поэтому ваше первое использование проекта неверно. Если я думаю, что это правильно, у вас есть координаты, которые есть в AGD66, поэтому вы должны сначала назначить эту проекцию, а затем вы можете преобразовать. Вы не можете выполнять базовые преобразования с project (), но spTransform () может.

Я думаю, тебе нужно это: 

pts = read.table(text = "lon        lat
147.1428   -43.49083", header = TRUE)

## assign original coordinate system
pts66 = SpatialPoints(cbind(pts$lon,pts$lat), CRS("+init=epsg:4202"))

## Then to transform it to WGS84
pts84 = spTransform(pts66, CRS("+init=epsg:3033"))


pts66
SpatialPoints:
     coords.x1 coords.x2
[1,]  147.1428 -43.49083
Coordinate Reference System (CRS) arguments: +init=epsg:4202 +proj=longlat +ellps=aust_SA
+no_defs 

pts84
SpatialPoints:
     coords.x1 coords.x2
[1,]  11126605   2971806
Coordinate Reference System (CRS) arguments: +init=epsg:3033 +proj=lcc +lat_1=-68.5     +lat_2=-74.5
+lat_0=-50 +lon_0=70 +x_0=6000000 +y_0=6000000 +ellps=WGS84 +datum=WGS84 +units=m +no_defs
+towgs84=0,0,0

Это означает, что pts66 не изменились по сравнению с их исходными значениями, но у них есть правильные метаданные для следующего шага, который преобразует их в вашу цель (что, кстати, является конформно-конической Ламберта). Вам может потребоваться немного больше расследования, чтобы выяснить, что требуется. 

CRS("+init=epsg:4202")
CRS arguments:
+init=epsg:4202 +proj=longlat +ellps=aust_SA +no_defs 

CRS("+init=epsg:3033")

CRS arguments:
+init=epsg:3033 +proj=lcc +lat_1=-68.5 +lat_2=-74.5 +lat_0=-50
+lon_0=70 +x_0=6000000 +y_0=6000000 +ellps=WGS84 +datum=WGS84 +units=m
+no_defs +towgs84=0,0,0

Ваш первоначальный проект () неправильно пытался преобразовать longlat на WGS84 в longlat на AGD66, но эта функция не может этого сделать, поэтому это просто добавило путаницу в миксе. Элемент данных не является проекцией, он является критической частью определения проекции, и в этом смысле «longlat on AGD66» является проекцией, как и «конформная коническая проекция Ламберта на WGS84».

...