[Qgis-user] Re: [Qgis-developer] Wrong reprojection if using different datums?

[Agustin Lobo]

Micha,

It would be great if just adding the towgs84
parameters were enough, but in any case this problem
has to be clearly stated for the user who
normally knows nothing about the towgs84 parameters
(i.e., me 1 week ago)
BTW, I think that all CRS need towgs84 parameters
for proj4 to actually make a shift from a non-WGS84
datum to WGS84). At appropriate places, a note should
be placed: "Warning: CRS transforms involving change of
original datum to WGS84 must have non zero parameters
in the +towgs field. User must confirm that these
parameters are correct. Possible sources:
http://www.asprs.org/resources/grids/
http://earth-info.nga.mil/GandG/coordsys/onlinedatum "

According to Roger Bivand, the note could be:
"ignoring datum declarations can seriously damage your health near cliffs"

Actually, I found this problem in R (rgdal), then in QGIS when I
looked for a reference, This is the result
of the testing I've done with different programs:

Original point set in QGIS UTM 31N ED50 (epsg 23031)
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 than spTransform() once
towgs84 is properly set.

Just to note in the passing: initially I made the error of not including
the "," in the
+towgs84 field in rgdal, which was causing wrong results and
the waist of several days...

More on this problem in rgdal
http://n2.nabble.com/spTransform%3A-wrong-results---tp3238715p3238715.html

Agus

Micha Silver wrote:
> Paolo Cavallini wrote:
> 
>> Agustin Lobo ha scritto:
>>  
>>> different in both datums? I suspect that proj4 (behind QGIS)
>>> is not using the +towgs= field.
>>>
>>> Any opinions on this problem?
>>>     
>>
>> Hi Agus.
>> I think QGIS is mostly ignoring +twgs parameters. This is the root of
>> several bugs, and I think is one of the major issues to be solved.
>> Mentioned also in our bug squashing initiative:
>> http://www.qgis.org/wiki/Bugs#Projections
>> Unfortunately, it did not attract much specific funding jet.
>> All the best.
>>   
> The Israeli CRS also uses a +towgs parameter. This did not appear in the 
> srs.db until about a year ago. After I asked on the list, the entry for 
> EPSG:2039 was fixed in the srs.db (now equivalent to the proj4 
> definition) and since then QGIS has done reprojection with datum 
> correction properly (i.e. using v.proj in GRASS gives the same results 
> as QGIS projection).
> Agus- can you check the entry for your EPSG code in 
> /usr/share/qgis/resources/srs.db (sqlitebrowser) and compare with the 
> same entry in /usr/share/proj/epsg?
>