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

[Agustin Lobo]

Roger,

This part of your message:
"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."

is critical, could not be included in the spTransform() help page?
Actually a whole paragraph regarding changes of datum
would be very helpful, including the example.

Also, I understand that I must search the magic towgs84 parameters in
the "Grids & Datums" essays
 http://www.asprs.org/resources/grids/

I'll try and let you know. For the future, no way to compile the
towgs84 parameter so that could be automatically used
by spTransform? This must be what programs like GlobalMapper do, as
you just define input and output datum.

Thanks for your help!

Agus

2009/7/10 Roger Bivand <Roger.Bivand at nhh.no>:

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


--
Dr. Agustin Lobo
Institut de Ciencies de la Terra "Jaume Almera" (CSIC)
Lluis Sole Sabaris s/n
08028 Barcelona
Spain
Tel. 34 934095410
Fax. 34 934110012
e-mail Agustin.Lobo at ija.csic.es
http://www.ija.csic.es/gt/obster



-- 
Dr. Agustin Lobo
Institut de Ciencies de la Terra "Jaume Almera" (CSIC)
Lluis Sole Sabaris s/n
08028 Barcelona
Spain
Tel. 34 934095410
Fax. 34 934110012
e-mail Agustin.Lobo at ija.csic.es
http://www.ija.csic.es/gt/obster