spTransform: wrong results??

Roger,

First at all, I think that the problem is solved.
Note that I was setting the CRS in QGIS, thus it is unequivocally
specified as EPSG 23031, then read into R through readOGR()

Anyway, the problem is that I was setting the towgs84 field as:
proj4string(wpcast300SPDF) <- CRS(paste(" +proj=utm +zone=31 +ellps=intl 
+units=m +towgs84=-84 -107 -120 +no_defs"))

while it must be:
proj4string(wpcast300SPDF) <- CRS(paste(" +proj=utm +zone=31 +ellps=intl 
+units=m +towgs84=-84,-107,-120 +no_defs"))

Note the "," !!!

With this I get:

Origianl point set in QGIS UTM 31N ED50
POINT(449919.586269 4631874.109627)	

TRANSFORMS TO LONGITUDE,LATITUDE:

Official reference (on-line tool of the Catalan Institute of Cartography)
ICC WGS84 2.395717586,41.835328397

GlobalMapper  ED50  2.39686371? E	41.83643863? N
GlobalMapper  WGS84 2.39572870? E	41.83535625? N

R (spTransform() in rgdal)
ED50 no towgs84 field specified
WP1	2.39686370617169	41.8364386271064

WGS84 no towgs84 field specified
WP1	2.39686370617169	41.8364386271064

WGS84 with towgs84=-87, -98,-121
(http://earth-info.nga.mil/GandG/coordsys/onlinedatum/CountryEuropeTable.html)
WP1  2.395729  41.83536

WGS84 with towgs84=-84,-107,-120
(http://www.asprs.org/resources/grids/)
WP1  2.395619  41.83535

The equivalent syntax in cs2cs:
cs2cs -v -f "%.9f" +init=epsg:23031 +towgs84=-87,-98,-121 +to 
+init=epsg:4326

provides identical output

Also note that the parameters found in http://earth-info.nga.mil
are closer to ICC output and identical to GlobalMapper output

I think that without the "," spTransform() was probably using just
the first parameter.

Thanks for your help!

Agus

On Sun, 12 Jul 2009, Agustin Lobo wrote:

Here you have some results with
different software for the transformation
of one single  UTM 31N ED50 coordinate to
lon,lat ED50 and lon,lat WGS84.
My understanding of this problem seems to
be decreasing as I do more testing.

QGIS UTM 31N ED50

This isn't helpful, I'm afraid. You are starting from an unknown, not a 
known. It isn't imposible that QGIS doesn't give a +towgs84 for 
so-called ED50, because ED50 is *not* one datum but many.

Start with an unequivocally known point in a known datum, such as a GPS 
reading in geographical coordinates and WGS84. Work from there, and it 
should get clearer.

Datum definitions are not easy at all, because organisations typically 
always worked within a single datum, and never needed to transform 
(until oil was found across two different jurisdictions using different 
datum standards, hence EPSG).

Roger

POINT(449919.586269 4631874.109627)

TRANSFORMS TO LONGITUDE,LATITUDE:

QGIS ED50
POINT(2.396864 41.836439)
QGIS WGS84
POINT(2.396864 41.836439)

GlobalMapper  ED50  2.39686371? E    41.83643863? N
GlobalMapper  WGS84 2.39572870? E    41.83535625? N

TNT WGS84
POINT(2.396864 41.836439)
TNT ED50
POINT(2.396864 41.836439)

rgdal ED50 no towgs84 field specified
WP1    2.39686370617169    41.8364386271064

rgdal WGS84 no towgs84 field specified
WP1    2.39686370617169    41.8364386271064

rgdal WGS84 with towgs84=-87, -98,-121
WP1    2.39690750837097    41.8361433350555



Agus