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

[pkelly]

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

[sgamperl]

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

[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

[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

[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.

[sgamperl]

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

Stefan

[warmerdam]

Some related tweaking to this logic done for #3737.

[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.