[R-sig-Geo] spTransform problem

[Roger Bivand]

On Wed, 15 Nov 2006, Duncan Golicher wrote:

> Can anyone tell me what I am doing wrong when using spTransform?
> I need to transform points in Southern Mexico from UTM NAD27 to Long Lat 
> WGS84.

Because spTransform uses PROJ.4, this is a PROJ.4 problem. I wonder 
whether PROJ.4 provides NAD grids for Southern Mexico, or whether you will 
rather need to find the +towgs= arguments directly? If PROJ.4 doesn't 
provide the NAD transformation for your location, I'd search. Without 
knowing if it helps, this looks relevant:

http://article.gmane.org/gmane.comp.gis.proj-4.devel/1582/match=mexico

with +towgs84=-12,130,190. I don't think just putting +datum=NAD27 will 
work, as you've found. Is there a Grids&Datums for Mexico (can't see one 
on: http://www.asprs.org/resources/grids/)? 

Roger

> 
> Three small functions.
> 
> NAD27_WGS84ll<-function(points){
> coordinates(points)<-~x+y
> points at proj4string<-CRS("+proj=utm +zone=15 +ellps=clrk66 +datum=NAD27 
> +units=m")
> spTransform(points, CRS("+proj=longlat +datum=WGS84"))}
> 
> NAD27_NAD27ll<-function(points){
> coordinates(points)<-~x+y
> points at proj4string<-CRS("+proj=utm +zone=15 +datum=NAD27 +units=m")
> spTransform(points, CRS("+proj=longlat +datum=NAD27"))}
> 
> NAD27_WGS84<-function(points){
> coordinates(points)<-~x+y
> points at proj4string<-CRS("+proj=utm +zone=15 +datum=NAD27 +units=m")
> spTransform(points, CRS("+proj=utm +zone=15 +datum=WGS84 +units=m"))}
> 
> When I test them I get strange results.
> 
> x<-c(541987,541970)
> y<-c(1846598,1846799)
> 
> points<-data.frame(x,y)
> 
> NAD27_WGS84ll(points)
> 
> SpatialPoints:
>              x        y
> [1,] -92.60617 16.70272
> [2,] -92.60633 16.70454
> Coordinate Reference System (CRS) arguments:  +proj=longlat +datum=WGS84 
> +ellps=WGS84 +towgs84=0,0,0
> 
> NAD27_NAD27ll(points)
> 
> SpatialPoints:
>              x        y
> [1,] -92.60617 16.70272
> [2,] -92.60633 16.70454
> Coordinate Reference System (CRS) arguments:  +proj=longlat +datum=NAD27 
> +ellps=clrk66 +nadgrids=@conus, at alaska, at ntv2_0.gsb, at ntv1_can.dat
> 
> The results are the same, even though one function involves a datum 
> shift and the other doesn't.
> 
> The difference between NAD27 and WGS84 in Southern Mexico should be 
> around 200m N and 20m E. But I get
> 
> points-as.data.frame(NAD27_WGS84(points))
> 
>           x         y
> 1 0.5857494 -112.2658
> 2 0.5855395 -112.2768
> 
> 
> Thanks for any advice.
> 
> Duncan Golicher
> 

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no