Wrong reprojection with datum shift (Germany/ETRS)

[stoeckmann]

I updated my MapServer from 5.0.0 to 5.4.2, and now transformations between two projected coordinates systems (Germany and ETRS) are generally calculated wrong. The following PHP MapScript example shows it:

# --- direct reprojection ---
$Prj1 = ms_newProjectionObj('init=epsg:31468'); $Prj2 = ms_newProjectionObj('init=epsg:25833'); $pt = ms_newPointObj();
$pt->setXY(4588000, 5942000);
$pt->project($Prj1, $Prj2);
print $pt->x." ".$pt->y;

# --- indirect way using WGS ---
$PrjWGS = ms_newProjectionObj('init=epsg:4326'); $pt = ms_newPointObj();
$pt->setXY(4588000, 5942000);
$pt->project($Prj1, $PrjWGS);
$pt->project($PrjWGS, $Prj2);
print "<br>".$pt->x." ".$pt->y;

While MS 5.0.0 returns two identical coordinates, MS 5.4.2 returns this:
389470.239637 5940711.84385
389348.972139 5940544.41103

The upper coordinates are wrong.

[sdlime]

The version of Proj.4 is the same between the two tests? Moving to the Proj component...

Steve

[stoeckmann]

Proj had been automatically updated from version 4.4.9 to 4.6.1. I replaced the new libproj.so.0.5.5 with the old libproj.so.0.5.0 by hand, and now the calculation is correct again!

So I compiled the latest Proj version 4.7.1 and tested the libproj.so.0.6.6, but it calculates WRONG, too!

I read at the homepage of Proj that the API has been changed, especially it separates now between simple transformation (from geographic to projected coordinates - one datum only) and complex reprojection (transformation from any projected to any projected coordinate system - two datums). Maybe the MapServer sources has to be modified for the new Proj API to support two-datum-reprojetion?

[sdlime]

Really need to get Frank's input here...

Steve

[warmerdam]

From the PROJ.4 FAQ:

​http://trac.osgeo.org/proj/wiki/FAQ#HowdoIdebugproblemswithNAD27NAD83datumshifting

In particular note:

""" The cs2cs command and the underlying pj_transform() API know how to do a grid shift as part of a more complex coordinate transformation; however, it is imperative that both the source and destination coordinate system be defined with appropriate datum information. That means that implicitly or explicitly there must be a +datum= clause, a +nadgrids= clause or a +towgs84= clause. For instance "cs2cs +proj=latlong +datum=NAD27 +to +proj=latlong +ellps=WGS84" won't work because defining the output coordinate system as using the ellipse WGS84 isn't the same as defining it to use the datum WGS84 (use +datum=WGS84). If either the input or output are not identified as having a datum, the datum shifting (and ellipsoid change) step is just quietly skipped! """

In this case EPSG 31468 normally expands to:

+proj=tmerc +lat_0=0 +lon_0=12 +k=1 +x_0=4500000 +y_0=0 +ellps=bessel +datum=potsdam +units=m +no_defs

which has a datum definition, but EPSG 25833 normally expands to:

+proj=utm +zone=33 +ellps=GRS80 +units=m +no_defs

which does *not* have a datum definition. So if you go directly from 31468 to 25833 you will now get no datum shift, but if you go from 31468 to 4326 you will get a datum shift (since 4326 also has a valid datum setting) but then the shift from 4326 to 25833 does not have a shift.

Note that the behavior with regard to when a datum shift is applied changed in PROJ 4.6.0 due to many problems with datum shifts being applied in bogus ways for coordinate systems that didn't really know their datum definition. I believe what you really need to do is establish an appropriate datum shift for EPSG 25833, and update the epsg init file accordingly. Then you should get consistent results whether going via WGS84 or not, and more importantly you will get correct results.

If you wish to assume that EPSG 25833 (based on ETRS89) is effectively the same datum as WGS84 then change +ellps=GRS80 to +datum=WGS84 in the definition in the epsg init file.

PS. this is in no way a MapServer problem, and it would be dangerous to have MapServer try and work around PROJ.4's behavior.