[R-sig-Geo] spTransform: wrong results??

[Roger Bivand]

On Fri, 10 Jul 2009, Agustin Lobo wrote:

> I'm finding conflicting results at comparing
> results of converting UTM ED50 coordinates
> to lon,lat WGS84 using spTransform() (thus, proj4)
> and other programs such as GlobalMapper and Compegps.
>
> My main surprise is that spTransform() yields the same result
> whatever the ellps field for the inverse transform.

Wrong tree to bark up, I'm afraid. The tags you need and are missing are 
+datum= and +towgs84= (from memory, you'll see from the EPSG tags). The 
+ellps= tag sets the ellipsoid, but doesn't change the default +datum= 
from WGS84. Note that ED50 is a collection of datum values, and in general 
EPSG does not give 3 or 7 parameter datum transformation values unless 
they are known with certainty. On the Proj.4 side, you'll find some 
documentation on:

http://svn.osgeo.org/metacrs/proj/trunk/proj/html/gen_parms.html

otherwise use makeEPSG() in rgdal to search the bundled EPSG database, and 
do check the "Grids & Datums" essays on:

http://www.asprs.org/resources/grids/

for further help in tracking down possible +towgs84= parameter values.

Repeating, unless you give a +datum=, it will be assumed to be WGS84, so 
you shouldn't be surprised when you get similar results.

Sorry to disappoint, but surveying and cartography have their own 
traditions, some of which are localised as well.

Hope this helps,

Roger

> This is
> what I'm doing with the meuse dataset as an example, you can
> see that the results of the inverse transform to lon,lat WGS84 are identical 
> to
> those to lon,lat ED50, both starting from UTM50 or from UTMWGS84
>
> At the end
> I include results with GlobalMapper for the first point. Results
> with GlobalMapper do differ for lon,lat WGS84 and lon,lat ED50
>
> Any help would be appreciated, as I'm doing all my work
> in UTM ED50 but have to pass the coordinates as lon,lat WGS84
> for the flight.
>
> data(meuse)
> coordinates(meuse) <- c("x", "y")
> proj4string(meuse) <- CRS("+proj=stere +lat_0=52.15616055555555 
> +lon_0=5.38763888888889 +k=0.999908 +x_0=155000 +y_0=463000 +ellps=bessel 
> +units=m")
> meuse.utm <- spTransform(meuse, CRS("+proj=utm +zone=32 +ellps=WGS84"))
> coordinates(meuse.utm)[1:3,]
>            x       y
> [1,] 272571.9 5653978
> [2,] 272522.4 5653927
> [3,] 272661.2 5653899
>
> meuse.utmintl <- spTransform(meuse, CRS("+proj=utm +zone=32 +ellps=intl"))
> coordinates(meuse.utmintl)[1:3,]
>            x       y
> [1,] 272561.0 5654094
> [2,] 272511.5 5654043
> [3,] 272650.3 5654015
>
> meuse.lonlatWGS84WGS84 <- spTransform(meuse.utm, CRS("+proj=lonlat 
> +ellps=WGS84"))
> meuse.lonlatWGS84intl  <- spTransform(meuse.utm, CRS("+proj=lonlat 
> +ellps=intl"))
> coordinates(meuse.lonlatWGS84WGS84)[1:3,]
>            x        y
> [1,] 5.759053 50.99238
> [2,] 5.758380 50.99190
> [3,] 5.760373 50.99171
>
> coordinates(meuse.lonlatWGS84intl)[1:3,]
>            x        y
> [1,] 5.759053 50.99238
> [2,] 5.758380 50.99190
> [3,] 5.760373 50.99171
>
> meuse.lonlatintlWGS84 <- spTransform(meuse.utmintl, CRS("+proj=lonlat 
> +ellps=WGS84"))
> meuse.lonlatintlintl  <- spTransform(meuse.utmintl, CRS("+proj=lonlat 
> +ellps=intl"))
> coordinates(meuse.lonlatintlWGS84)[1:3,]
>            x        y
> [1,] 5.759053 50.99238
> [2,] 5.758380 50.99190
> [3,] 5.760373 50.99171
>
> coordinates(meuse.lonlatintlintl)[1:3,]
>            x        y
> [1,] 5.759053 50.99238
> [2,] 5.758380 50.99190
> [3,] 5.760373 50.99171
>
> Test with Global Mapper
>
> UTM ED50 272561.0 5654094 to lon,lat:
> WGS84   5° 45' 27.5632" E	50° 59' 29.5996" N
> ED50    5° 45' 32.5898" E	50° 59' 32.5665" N
>
> UTM WGS84 272571.9 5653978 to lon,lat:
> WGS84   5° 45' 32.5900" E	50° 59' 32.5621" N
> ED50    5° 45' 37.6165" E	50° 59' 35.5289" N
>
> Also, utmWGS84 to utmED50 is different in GlobalMapper
>
> UTM WGS84  272571.9     5653978
> UTM ED50   272662.989	5654181.171
>
>
> Thanks,
>
> Agus
>

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