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
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
