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

[Agustin Lobo]

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


Roger Bivand wrote:
> 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
>>
>>
>>
>