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.biv...@nhh.no
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
