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

[Roger Bivand]

On Fri, 10 Jul 2009, Agustin Lobo wrote:

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

I'll try to do as you ask and add details and links. The definitions of 
the CRS class also could refer to this. In addition, the website at

http://spatialreference.org/

also provides som information based on EPSG.

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

No way. Choosing the 3 or 7 parameter datum transforms has to be the 
responsibility of the user when there is no unequivocal authority. Having 
a known ground control point (a GPS reading from a known point that has a 
known CRS) can be used to see which CRS settings let you reproduce it. But 
with legacy maps, this isn't so easy. GlobalMapper must I suppose be 
making choices and using authorities that it acknowledges, but reading G&D 
does not make me optimistic.

Roger

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

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