Error when re-projecting spatial points using spTransform in rgdal R

Question

G'day,
I have a large number of lon/lat coordinates that are in the CRS Australian Geodetic Datum 66/84 (AGD66 for brevity). I want to change these coordinates from AGD66 to WGS84 because there is about a 200m difference between them and I have other coordinates and layers in WGS84. I have tried to do this by:

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"

Does anyone know why I get this error or have any suggestions for how I can change these coordinates from AGD66 to WGS84? Thanks for your help in advance.

Cheers,
Adam

Answer 1

I have removed part of an incorrect answer.

The function project() cannot do datum conversions so you may have a problem there and I think what you have is wrong.

The problem is that you can only project() from/to longlat on WGS84, so your first use of project is incorrect. If I guess this right, you have coordinates that are in AGD66 so you must first assign that projection and then you can transform. You cannot do datum transformations with project(), but spTransform() can.

I think you need this:

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 

This means that pts66 are not changed from their original values, but they have the metadata correct for the next step which transforms them to your target (which is Lambert Conformal Conic btw). You might need a bit more investigation to figure out what's required.

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 

Your original project() was incorrectly attempting to transform from longlat on WGS84 to longlat on AGD66, but that function cannot do that so it's just added confusion in the mix. A datum is not a projection, it is a critical part of a projection's definition and in this sense "longlat on AGD66" is a projection just as "Lambert Conformal Conic on WGS84" is.