+towgs84 parameter missing for epsg:28992 (Dutch coordinate system)

[janh]

The seven +towgs84 values for epsg:28992 are:

565.2369,50.0087,465.658,-0.406857330322398,0.350732676542563,-1.8703473836068,4.0812

They are derived from the epsg-database 7.1 (21 May 2009), coordinate transformation #62: Amersfoort to WGS84(3) (the most recent formula for the Netherlands).

The three rotations are converted from radians to degree and sign-reversed.

[warmerdam]

Jan,

I have upgraded the libgeotiff machinery for selecting datum shifts and now it has selected:

+towgs84=565.04,49.91,465.84,-1.9848,1.7439,-9.0587,4.0772

A review of the possible datum shifts for Amersfoort shows four possible transformations. The one you provide is operation 15934. But operation 4833 includes the remark text:

Replaces Amersfoort to WGS 84 (3) (code 15934)

So I believe it is considered the best currently available transformation. Let me know if you do not agree and feel we ought to override this in some fashion.

BTW, the full WKT for 28992 is now:

PROJCS["Amersfoort / RD New",
    GEOGCS["Amersfoort",
        DATUM["Amersfoort",
            SPHEROID["Bessel 1841",6377397.155,299.1528128,
                AUTHORITY["EPSG","7004"]],
            TOWGS84[565.04,49.91,465.84,-1.9848,1.7439,-9.0587,4.0772],
            AUTHORITY["EPSG","6289"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.01745329251994328,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4289"]],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    PROJECTION["Oblique_Stereographic"],
    PARAMETER["latitude_of_origin",52.15616055555555],
    PARAMETER["central_meridian",5.38763888888889],
    PARAMETER["scale_factor",0.9999079],
    PARAMETER["false_easting",155000],
    PARAMETER["false_northing",463000],
    AUTHORITY["EPSG","28992"],
    AXIS["X",EAST],
    AXIS["Y",NORTH]]

and the PROJ.4 definition looks like:

+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.04,49.91,465.84,-1.9848,1.7439,-9.0587,4.0772 +units=m +no_defs

The changes in libgeotiff are in r1819 and r1820.

[janh]

I spent the better part of the afternoon with the @#$% EPSG database, version 4.9, 24 nov 2009, to figure out the right Datum shift parameters for the Dutch coordinate system (EPSG:28992) to WGS84. There are 10 codes, four of which are unusable (center of rotation is Amersfoort). The six others, with the center of the ellipsoid as center, duplicate each other. The three different codes that eventually qualify for a PROJ-string (ETRS1, 3 and 5, or WGS2, 3 and 4) yield results that are less than 1 cm different. Attached an Excel spreadsheet with the computations (you may discard the macros, Tools/Macro/Security/Medium).

The resulting PROJ-string is based on ETRS5 or WGS4 (codes 4830 or 4833):

+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.999908 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +towgs84=565.4174,50.3319,465.5542,-0.398957388243134,0.343987817378283,-1.87740163998045,4.0725 +no_defs

How do you choose the right formula from the EPSG database for libgeotiff, automatically or by hand?

Jan

[warmerdam]

Jan,

The build_pcs.py script identifies four transformations from Amersfoort to WGS84 that are usable (1112,1672,4833 and 15934). Normally the script will select the first datum shift it encounters that has the largest area of use. However, I provide a mechanism for overriding the default selection in the datum_shift_pref.csv file. I was incorrectly forcing it to use 1672 but on more careful review of the notes fields your choice of 4833 is most appropriate. I have adjusted this in libgeotiff trunk (r1826).

Thanks for the careful and painful research.

[janh]

I couldn't find the scripts you mentioned, but the default is error prone:

There are ten transformations for Amersfoort:

• Usable with PROJ (rotated around center of ellipse): 1751 ETRS(=1762 WGS84), 15739 ETRS (=15734 WGS84), 4830 ETRS (=4833 WGS84)
• Not usable with PROJ (rotated around Amersfoort): 1066 (ETRS), 15740 (ETRS), 4831 (ETRS), 1112 (WGS84). These are the equivalents of the three usable ETRS transformations above.

The script seems to look only for transformations to WGS84, not to ETRS, although they are equivalent. Even so, it gets one wrong (1112) and one I couldn't find (15934).

Should the script discover the ETRS transformations, and if so, how can it distinguish between usable and non-usable versions? I don't think an automated solution is possible. Even the documentation provided by EPSG is insufficient on its own. Perhaps a WIKI for datum transformations, ordered by EPSG number?

[mikrit]

Some ideas for improvements:

Frank wrote:

Normally the script will select the first datum shift it encounters that has the largest area of use. However, I provide a mechanism for overriding the default selection in the datum_shift_pref.csv file.

It may be possible to automatically discard datum shifts that have been replaced by others. There is a Supersession table in the EPSG database, where datum shifts that have replacements are marked as Superseeded; this is a weaker form of Deprecation. (The different meanings are explained in OGP Guidance Note 7.1.)

Alternatively, one could look for the string "Replaced by" in the comment field for each datum shift, but free-format text is of course less reliable.

[warmerdam]

Mikael,

An excellent suggestion! I have incorporated it in build_pcs.py (r1829). Note that superseded datum shifts are select included in datum_shift.csv, but they are not selected as preferred if they have been superseded.

I was also able to remove half the overridden preferences, including Amersfoort with this rule in place.

I will comment on Jan's questions, and suggestions after some additional review.

[warmerdam]

Replying to janh:

I couldn't find the scripts you mentioned, but the default is error prone:

There are ten transformations for Amersfoort:

• Usable with PROJ (rotated around center of ellipse): 1751 ETRS(=1762 WGS84), 15739 ETRS (=15734 WGS84), 4830 ETRS (=4833 WGS84)
• Not usable with PROJ (rotated around Amersfoort): 1066 (ETRS), 15740 (ETRS), 4831 (ETRS), 1112 (WGS84). These are the equivalents of the three usable ETRS transformations above.

The script seems to look only for transformations to WGS84, not to ETRS, although they are equivalent. Even so, it gets one wrong (1112) and one I couldn't find (15934).

Should the script discover the ETRS transformations, and if so, how can it distinguish between usable and non-usable versions? I don't think an automated solution is possible. Even the documentation provided by EPSG is insufficient on its own. Perhaps a WIKI for datum transformations, ordered by EPSG number?

[warmerdam]

Ignore the last attempt to reply ...

Replying to janh:

I couldn't find the scripts you mentioned, but the default is error prone:

​http://svn.osgeo.org/metacrs/geotiff/trunk/libgeotiff/csv/build_pcs.py

There are ten transformations for Amersfoort:

• Usable with PROJ (rotated around center of ellipse): 1751 ETRS(=1762 WGS84), 15739 ETRS (=15734 WGS84), 4830 ETRS (=4833 WGS84)
• Not usable with PROJ (rotated around Amersfoort): 1066 (ETRS), 15740 (ETRS), 4831 (ETRS), 1112 (WGS84). These are the equivalents of the three usable ETRS transformations above.

The script seems to look only for transformations to WGS84, not to ETRS, although they are equivalent. Even so, it gets one wrong (1112) and one I couldn't find (15934).

Jan,

Can you explain what you mean about 1112 being unsuitable and rotated around Amersfoort? How do you determine this? I assume that any transformation from the GCS to WGS84 that uses methods 9603, 9606 or 9607 can be used by PROJ and other packages for a datum shift.

Should the script discover the ETRS transformations, and if so, how can it distinguish between usable and non-usable versions?

We could potentially extend the script to recognise that transformations to ETRS and some other euqivelent-to-WGS84 datums (like NAD83) could also be used, but it appears that EPSG is already providing copies of some of these as transformations to WGS84 so perhaps it is not necessary.

I don't think an automated solution is possible. Even the documentation provided by EPSG is insufficient on its own. Perhaps a WIKI for datum transformations, ordered by EPSG number?

Ideally, I'd like to ensure that such knowledge gets into "machine understandable" format - for instance to select a preferred datum shift in the datum_shift.csv table.

[janh]

Replying to warmerdam:

Can you explain what you mean about 1112 being unsuitable and rotated around Amersfoort? How do you determine this? I assume that any transformation from the GCS to WGS84 that uses methods 9603, 9606 or 9607 can be used by PROJ and other packages for a datum shift.

From literature, this can not be derived from the EPSG database. You can plug these parameters in PROJ, but the results are incorrect. The other three transformations can be used, and give results that are equal up to around a centimeter. In practice it doesn't matter very much which one you use, so use the last one.

We could potentially extend the script to recognise that transformations to ETRS and some other euqivelent-to-WGS84 datums (like NAD83) could also be used, but it appears that EPSG is already providing copies of some of these as transformations to WGS84 so perhaps it is not necessary.

True

Ideally, I'd like to ensure that such knowledge gets into "machine understandable" format - for instance to select a preferred datum shift in the datum_shift.csv table.

True, but I would like to have a place where we can see why a particular shift was prefered. Perhaps the best starting point would be Cliff Mugnier's columns (​http://www.asprs.org/resources/GRIDS/). I'll contact him for a version ordered by EPSG number.

Jan