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

[Roger Bivand]

On Fri, 10 Jul 2009, Roger Bivand wrote:

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

To give this more meat, see below for the development:

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

library(rgdal)
data(meuse)
meuse_no_towgs84 <- meuse
meuse_towgs84 <- meuse
coordinates(meuse_no_towgs84) <- c("x", "y")
coordinates(meuse_towgs84) <- c("x", "y")
EPSG <- make_EPSG()
EPSG[grep("28992", EPSG$code),]
proj4string(meuse_no_towgs84) <- CRS("+init=epsg:28992")
proj4string(meuse_towgs84) <- CRS(paste("+init=epsg:28992",
  "+towgs84=565.237,50.0087,465.658,-0.406857,0.350733,-1.87035,4.0812"))
# chunk 36, http://www.asdar-book.org/book/die_mod.R
# see ASDAR p. 96 for authority
# and G&D February 2003, for similar, alternative 7-parameters in right, 
# not left-hand rotation
proj4string(meuse_no_towgs84)
proj4string(meuse_towgs84)
utm_crs_ellps <- CRS("+proj=utm +zone=32 +ellps=WGS84")
utm_crs_datum <- CRS("+proj=utm +zone=32 +datum=WGS84")
coordinates(spTransform(meuse_no_towgs84, utm_crs_ellps))[1:3,]
coordinates(spTransform(meuse_towgs84, utm_crs_ellps))[1:3,]
coordinates(spTransform(meuse_no_towgs84, utm_crs_datum))[1:3,]
coordinates(spTransform(meuse_towgs84, utm_crs_datum))[1:3,]

Only the final version works. As expected, when +towgs84= is not set, no 
datum shift occurs, and the same happens when the target datum is not 
specified. In general +datum= is to be prefered to +ellps=, because the 
datum always fixes the ellipsoid, but the ellipsoid never fixes the datum.

Note that your choice of projection is not correct, there are two:

projs <- projInfo()
projs[grep("stere", projs$name),]

and yours is "stere", not the correct "sterea" for this data. This is why

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

Please always use +proj=longlat, not the abbreviation, as a lot of code 
assumes only longlat and will misinterpret the other.

meuse_utm_towgs84 <- spTransform(meuse_towgs84, utm_crs_datum)
ll_datum <- CRS("+proj=longlat +datum=WGS84")
ll_ellps <- CRS("+proj=longlat +ellps=WGS84")
coordinates(spTransform(meuse_utm_towgs84, ll_datum))[1:3,]
coordinates(spTransform(meuse_utm_towgs84, ll_ellps))[1:3,]

Both are OK, but both were starting from WGS84 as datum. But if we try 
other datum values:

projInfo("datum")[5,]
ll_potsdam_datum <- CRS("+proj=longlat +datum=potsdam")
coordinates(spTransform(meuse_utm_towgs84, ll_potsdam_datum))[1:3,]

we get other coordinate values. This is discussed in ASDAR on pp. 82-87.

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

Since the data you are giving to GlobalMapper (I don't have it) are 
arbitrary, I cannot reproduce this. I said earlier that ED50 is not a 
standard, so it may well be that you are using one of the Spanish ED50 
datum values (or rather some +towgs84 parameter values), rather than the 
ED50 values for the Netherlands (which are given with an interval in G&D). 
ED50 is typically very arbitrary, and conversions to WGS84 may depend on 
local and unpublished decisions taken in mapping agencies, not 
infrequently military mapping agencies. This is what is described in G&D 
by Cliff Mugnier in his excellent forensic archeology. Reading the 
November 2006 G&D on Denmark is like a detective story, and most countries 
have well-kept secrets!

The actual WGS84 geographical coordinates for the first three points in 
meuse in dd:mm:ss are:

res <- spTransform(meuse_towgs84, ll_datum)
x <- as(dd2dms(coordinates(res)[1:3, 1]), "character")
y <- as(dd2dms(coordinates(res)[1:3, 2], TRUE), "character")
cbind(x, y)

Try:

res$long <- coordinates(res)[, 1]
res$lat <- coordinates(res)[, 2]

writeOGR(res, "meuse.kml", "meuse", driver="KML")

and check in GE (assuming that their coordinates are OK).

Hope this helps,

Roger

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