[R-sig-Geo] spTransform problem

[Roger Bivand]

On Thu, 16 Nov 2006, Roger Bivand wrote:

> 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/)? 

So far:

pts <- SpatialPoints(cbind(c(541987,541970), c(1846598,1846799)), 
 proj4string=CRS(
 "+proj=utm +zone=15 +ellps=clrk66 +towgs84=-12,130,190"))
coordinates(pts) - coordinates(spTransform(pts, CRS(
 "+proj=utm +zone=15 +datum=WGS84")))

gives:

     coords.x1 coords.x2
[1,]  18.65260 -200.7308
[2,]  18.65267 -200.7314

which is near your ballpark 20m E  and 200m N (well, it's -200m, but the 
absolute number agrees). I find that is very often necessary to search 
extensively to find relevant +towgs= arguments.

Roger


> 
> 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