Opened 10 years ago

Closed 9 years ago

Last modified 8 years ago

#3450 closed defect (fixed)

OSRExportToProj4() ignores towgs84 parameters if it recognises the datum name

Reported by: pkelly Owned by: warmerdam
Priority: high Milestone: 1.8.0
Component: OGR_SRS Version: 1.7.2
Severity: normal Keywords:
Cc: Markus Neteler, sgamperl

Description

A GRASS user noticed a problem whereby in a certain co-ordinate reference system manipulation using the GRASS module g.proj, towgs84 parameters he was specifying were being discarded. I've tracked this down to OSRExportToProj4() where if it recognises a datum name from the WKT it does not attempt to look for a TOWGS84 clause. There may be a good reason for this, so I submit the attached patch (against SVN) with caution.

As an example of how to reproduce the problem, consider the following WKT:

PROJCS["UTM Zone 6, Northern Hemisphere",

GEOGCS["NAD27",

DATUM["North_American_Datum_1927",

SPHEROID["Clarke 1866",6378206.4,294.978698213898,

AUTHORITY["EPSG","7008"]],

TOWGS84[-5,135,172,0,0,0,0], AUTHORITY["EPSG","6267"]],

PRIMEM["Greenwich",0,

AUTHORITY["EPSG","8901"]],

UNIT["degree",0.0174532925199433,

AUTHORITY["EPSG","9108"]],

AUTHORITY["EPSG","4267"]],

PROJECTIONTransverse_Mercator?, PARAMETER["latitude_of_origin",0], PARAMETER["central_meridian",-147], PARAMETER["scale_factor",0.9996], PARAMETER["false_easting",500000], PARAMETER["false_northing",0], UNIT["Meter",1]]

As it is, this gets converted to:

+proj=utm +zone=6 +ellps=clrk66 +datum=NAD27 +units=m +no_defs

With the attached patch it gets converted to

+proj=utm +zone=6 +ellps=clrk66 +datum=NAD27 +towgs84=-5,135,172,0,0,0,0 +units=m +no_defs

This is more useful for our purposes anyway. I can't see any obvious side effects, but I fear I'm not aware of all the implications of this change so would not want to wholeheartedly recommend it - in particular I'm not sure how it might interact with the block of code

else if( nEPSGGeogCS != -1 ) { [...] }

and what that refers to.

But I respectfully submit it for consideration anyway! (Relevant GRASS-user message: <http://lists.osgeo.org/pipermail/grass-user/2010-March/055146.html>)

Thanks,

Paul

Attachments (1)

towgs84.patch (2.0 KB) - added by pkelly 10 years ago.

Download all attachments as: .zip

Change History (10)

Changed 10 years ago by pkelly

Attachment: towgs84.patch added

comment:1 Changed 10 years ago by Markus Neteler

Cc: Markus Neteler added

comment:2 Changed 9 years ago by sgamperl

Cc: sgamperl added
Priority: normalhigh
Version: 1.6.21.7.2

In my special case I created an instance of the OGRSpatialReference class and imported a spatial reference system through the method SetFromUserInput?.

SetFromUserInput?("EPSG:31259");

The reference system was read correctly. The WKT is

PROJCS["MGI / Austria GK M34",

GEOGCS["MGI",

DATUM["Militar_Geographische_Institute",

SPHEROID["Bessel 1841",6377397.155,299.1528128,

AUTHORITY["EPSG","7004"]],

TOWGS84[577.326,90.129,463.919,5.137,1.474,5.297,2.4232], AUTHORITY["EPSG","6312"]],

PRIMEM["Greenwich",0,

AUTHORITY["EPSG","8901"]],

UNIT["degree",0.0174532925199433,

AUTHORITY["EPSG","9108"]],

AUTHORITY["EPSG","4312"]],

UNIT["metre",1,

AUTHORITY["EPSG","9001"]],

PROJECTIONTransverse_Mercator?, PARAMETER["latitude_of_origin",0], PARAMETER["central_meridian",16.33333333333333], PARAMETER["scale_factor",1], PARAMETER["false_easting",750000], PARAMETER["false_northing",-5000000], AUTHORITY["EPSG","31259"], AXIS["X",NORTH], AXIS["Y",EAST]]

Until here everything works fine. When I try to transform a coordinate from this reference system to another, the Proj4 library is used. Therefore a Proj4 string is created of the spatial reference system (OGRSpatialReference::exportToProj4() ).

The correct proj4 string would be: +proj=tmerc +lat_0=0 +lon_0=16.33333333333333 +k=1 +x_0=750000 +y_0=-5000000 +ellps=bessel +towgs84=577.326,90.129,463.919,5.137,1.474,5.297,2.4232 +units=m +no_defs

but the generated string is: +proj=tmerc +lat_0=0 +lon_0=16.33333333333333 +k=1 +x_0=750000 +y_0=-5000000 +ellps=bessel +datum=hermannskogel +units=m +no_defs

The datum "hermannskogel" is old and hardcoded in proj4, so it uses the wrong TOWGS84 parameter and I get the wrong transformation results.

The code was added in revision 17546 by the user rouault. My suggestion is not to use hardcoded datums from Proj4 because a user cannot modify the TOWGS84 parameter through a configuration file.

Is there a workaround to this issue?

Thanks,

Stefan

comment:3 Changed 9 years ago by pkelly

Hi Stefan, Well spotted that this was a recent regression in r17546 - I assumed it was there longer and we just hadn't come across it before. I believe it is incorrect behaviour of OGR to include the +datum in the PROJ.4 string if a TOWGS84 clause is present in the WKT - including the +datum is indeed useful if the datum is one recognised by PROJ.4 and *no* TOWGS84 parameters are available, otherwise it is a bad idea. The TOWGS84 parameters in the WKT should be assumed to be correct and not be overridden.

Paul

comment:4 Changed 9 years ago by sgamperl

Hey Paul!

According to the EPSG model their might be several possible transformation parameters (TOWGS84) depending on where you work in the spatial reference domain. To make a correct transformation from one datum to another you normally have to choose which parameter set is the correct one. Actually the parameter set is hardcorded in the gdal csv-files. I don't think that this is the optimal solution.

Also the hardcoded Proj4 datums has a specific parameter set (in the source), which is not the right parameters for every use cases.

Source pj_datums.c:


C_NAMESPACE_VAR struct PJ_DATUMS pj_datums[] = { /* id definition ellipse comments */ /* -- ---------- ------- -------- */ "WGS84", "towgs84=0,0,0", "WGS84", "", "GGRS87", "towgs84=-199.87,74.79,246.62", "GRS80",

"Greek_Geodetic_Reference_System_1987",

"NAD83", "towgs84=0,0,0", "GRS80",

"North_American_Datum_1983",

"NAD27", "nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat",

"clrk66",

"North_American_Datum_1927",

"potsdam", "towgs84=606.0,23.0,413.0", "bessel", "Potsdam Rauenberg 1950 DHDN", "carthage", "towgs84=-263.0,6.0,431.0", "clark80", "Carthage 1934 Tunisia", "hermannskogel", "towgs84=653.0,-212.0,449.0", "bessel", "Hermannskogel", "ire65", "towgs84=482.530,-130.596,564.557,-1.042,-0.214,-0.631,8.15", "mod_airy", "Ireland 1965", "nzgd49", "towgs84=59.47,-5.04,187.44,0.47,-0.1,1.024,-4.5993", "intl", "New Zealand Geodetic Datum 1949", "OSGB36", "towgs84=446.448,-125.157,542.060,0.1502,0.2470,0.8421,-20.4894", "airy", "Airy 1830", NULL, NULL, NULL, NULL };

Finally I have to say, that I agree to your opinion that the datum should only be used if there is no TOWGS84 parameter set specified.

I'v the GDAL version 1.7.2. Is there a workaround for this issue at the moment?

Thanks in advance,

Stefan

comment:5 Changed 9 years ago by Even Rouault

In trunk in r19903, add +towgs84= instead of +datum= if both information are available. This behaviour can be turned off by setting OVERRIDE_PROJ_DATUM_WITH_TOWGS84=NO

Test adapted in r19904

I'm a bit hesitant to backport this in 1.7 branch for now.

comment:6 Changed 9 years ago by sgamperl

Thank you for the prompt correction. Now the OSRExportToProj4 function works fine.

Stefan

comment:7 Changed 9 years ago by Even Rouault

Milestone: 1.8.0
Resolution: fixed
Status: newclosed

comment:8 Changed 9 years ago by warmerdam

Some related tweaking to this logic done for #3737.

comment:9 Changed 8 years ago by warmerdam

FYI, this change has caused some issues for the PROJ.4 project clients, as described in:

http://trac.osgeo.org/proj/ticket/122

Proj.4's epsg init file will in the future be generated with the config option set to YES.

Note: See TracTickets for help on using tickets.