New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Automatic datum shift when +nadgrids=@nul is not explicitly given #22

Closed
proj4-bot opened this Issue May 22, 2015 · 4 comments

Comments

Projects
None yet
2 participants
@proj4-bot

proj4-bot commented May 22, 2015

Reported by janh on 27 Nov 2008 14:27 UTC
An automatic datum shift to WGS84 seems to take place, even when not requested.

Test point 9018/375648 (Dutch RD, EPSG:28992 without WGS84 datum shift)
The source PROJ string is:
+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.999908 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs.
I tested this string both without and with +nadgrids=@null on four target projections:

  1. to latlon Bessel (+proj=longlat +ellps=bessel +no_defs)
    --> (without nadgrids in source) 3.291612 51.352060
    --> (with nadgrids null in source) 3.291612 51.351476

  2. to latlon Bessel with nadgrids null (+proj=longlat +ellps=bessel +no_defs +nadgrids=@null)
    --> (without nadgrids in source) 3.291612 51.352643
    --> (with nadgrids null in source) 3.291612 51.352060

  3. to latlon Bessel on the WGS84 datum (+proj=longlat +ellps=bessel +no_defs +datum=WGS84)
    --> (without nadgrids in source) 3.291612 51.352060
    --> (with nadgrids null in source) 3.291612 51.351476

  4. to latlon on the WGS84 datum (+proj=longlat +no_defs +datum=WGS84)
    --> (without nadgrids in source) 3.291612 51.352643
    --> (with nadgrids null in source) 3.291612 51.352060

  5. and 3) give the same results; so do 2) and 4), which seems to inicate that an automatic datum shift to WGS84 takes place, even when not requested. I'm not sure about the effects of +nadgrids=@null. Specifying +ellips=bessel without adding +nadgrids=@null seems to cause an automatic WGS84 datum shift, whether or no datum=WGS84 is specified (1a = 3a). If +nadgrids=@null is specified in the the source projection, but not in the target projection (1b = 3b), the latlon coordinates shift about 65m to the south. When +nadgrids=@null is omitted in the source projection, but added in the target projection, the resulting coordinates are shifted of 65M to the north (2a). Finally, when +nadgrids=@null is specified in both source and target projections, the datum shift is applied nevertheless (2b is equal to 1a and 3a). Finally, specifying a WGS84 datum shift for the target projection without the Bessel ellipsoide gives the same values as 2.

Migrated-From: https://trac.osgeo.org/proj/ticket/22

@proj4-bot

This comment has been minimized.

proj4-bot commented May 22, 2015

Comment by warmerdam on 27 Nov 2008 14:36 UTC
Jan,

Could you provide the exact command used for at least a couple of these? Did you use cs2cs? Also, was this 4.6.1?

@proj4-bot proj4-bot added C: Core and removed C: default labels May 22, 2015

@proj4-bot

This comment has been minimized.

proj4-bot commented May 22, 2015

Comment by janh on 27 Nov 2008 15:23 UTC
Sorry Frank, I have been mixing versions. I use Proj 4.6.1, but the tests were done with an older version of PostGIS, without realizing that it probably uses an older version of the library. I am not sure of all is well with 4.6.1 though. I reran the same test with cs2cs under proj 4.6.1 with the attached shell-file. It gives better results than my first run, but I suspect the projection still uses a WGS84 datum shift by default. Note that all transformations result in the same coordinates as when a datum shift to WGS84 is explicitly requested, except when a) a different datum is requested in the target projection and b) a +nadgrid=@null is specified in the source projection

@proj4-bot

This comment has been minimized.

proj4-bot commented May 22, 2015

Attachment added by janh on 27 Nov 2008 15:23 UTC
https://trac.osgeo.org/proj/attachment/ticket/22/testdatum.sh

@kbevers

This comment has been minimized.

Member

kbevers commented Jan 22, 2018

Tried to reproduce this with the script below.

echo "1."
echo 9018 375648 |./cs2cs -f %.6f +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.999908 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs +to +proj=longlat +ellps=bessel +no_defs

echo 9018 375648 |./cs2cs -f %.6f +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.999908 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +nadgrids=@null +no_defs +to +proj=longlat +ellps=bessel +no_defs

echo " "
echo "2."
echo 9018 375648 |./cs2cs -f %.6f +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.999908 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs +to +proj=longlat +ellps=bessel  +nadgrids=@null +no_defs

echo 9018 375648 |./cs2cs -f %.6f +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.999908 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +nadgrids=@null +no_defs +to +proj=longlat +ellps=bessel +nadgrids=@null +no_defs

echo " "
echo "3."
echo 9018 375648 |./cs2cs -f %.6f +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.999908 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs +to +proj=longlat +ellps=bessel +datum=WGS84 +no_defs

echo 9018 375648 0 |./cs2cs -f %.6f +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.999908 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +nadgrids=@null +no_defs +to +proj=longlat +ellps=bessel +datum=WGS84 +no_defs

echo " "
echo "4."
echo 9018 375648 |./cs2cs -f %.6f +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.999908 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs +to +proj=longlat +datum=WGS84 +no_defs

echo 9018 375648 |./cs2cs -f %.6f +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.999908 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +nadgrids=@null +no_defs +to +proj=longlat +ellps=bessel +no_defs

Which results in the following output:

1.
3.291612        51.352060 0.000000
3.291612        51.352060 0.000000

2.
3.291612        51.352060 0.000000
3.291612        51.352060 0.000000

3.
3.291612        51.352060 0.000000
3.291612        51.351477 699.339164

4.
3.291612        51.352060 0.000000
3.291612        51.352060 0.000000

I conclude that what is discussed in this issue is mostly fixed. The unexpected results in 3) is probably caused by the ellipsoid parameters being adjusted in pj_datum_transform() although the coordinates are not changed. Effectively this results in converting the coordinates from the WGS84 ellipsoid to the Bessel ellipsoid, which is backed up by

echo 3.291612 51.352060 0.000000 |./cct +proj=pipeline +step +proj=cart +ellps=WGS84 +step +proj=cart +ellps=bessel +inv
  3.2916120000   51.3514765406      699.3392           inf

kbevers added a commit to kbevers/proj.4 that referenced this issue Jan 22, 2018

Handle ellipsoid parameters correctly when using +nadgrids=@null. Fixes
OSGeo#22.

Make sure to not change ellipsoid parameters to WGS84 when applying the
null grid. Coordinates will still refer to the input ellipsoid so we
keep the original parameters which in turn will be used when the
coordinates are transformated to/from cartesian/geocentric space.

Adjusted regression test material in nad/proj_outIGNF.dist slightly to
accomodate numerical differences at the mm level. The transformations
in question are at best accurate to about 1m so this shouldn't change
real world usage of these transformations.

@kbevers kbevers closed this in d0dbf48 Jan 23, 2018

kbevers added a commit to kbevers/proj.4 that referenced this issue Mar 1, 2018

Revert fix to OSGeo#22
The fix in OSGeo#22 solved the problem at hand and doing what was expected
from the specified parameters. Unfortunately it also removed the slightly
hacky "feature" that makes the web mercator work in pj_transform. The
web mercator is special since the latitude is computed on the ellipsoid,
but behaves as if if was defined on a sphere. Hence it is problematic to
change the ellipsoid parameters when using the web mercator, even though
that is the geodetically correct thing to do. The web mercator is used in
more or less any web mapping application and is thus one of the most
frequently used transformations in PROJ. This justifies re-introducing
the minor bug reported in OSGeo#22.

The problem will have to be taken care of properly when pj_transform
is removed from the library in favour of the transformation pipelines
based API.

This was referenced Mar 1, 2018

kbevers added a commit that referenced this issue Mar 5, 2018

kbevers added a commit that referenced this issue Mar 5, 2018

Revert fix to #22
The fix in #22 solved the problem at hand and doing what was expected
from the specified parameters. Unfortunately it also removed the slightly
hacky "feature" that makes the web mercator work in pj_transform. The
web mercator is special since the latitude is computed on the ellipsoid,
but behaves as if if was defined on a sphere. Hence it is problematic to
change the ellipsoid parameters when using the web mercator, even though
that is the geodetically correct thing to do. The web mercator is used in
more or less any web mapping application and is thus one of the most
frequently used transformations in PROJ. This justifies re-introducing
the minor bug reported in #22.

The problem will have to be taken care of properly when pj_transform
is removed from the library in favour of the transformation pipelines
based API.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment