To give this more meat, see below for the development:
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
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.
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
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