Skip to content
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.

Sign up for GitHub

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

Jump to bottom

Incorrect data for EPSG:31370 (Belgian Lambert 72) / non-optimal transformation chosen for BD72 to WGS84 #1580

Closed
Wim-De-Clercq opened this issue Sep 3, 2019 · 8 comments
Closed

Incorrect data for EPSG:31370 (Belgian Lambert 72) / non-optimal transformation chosen for BD72 to WGS84 #1580

Wim-De-Clercq opened this issue Sep 3, 2019 · 8 comments
Labels
pinned Prevent stale-bot from closing an issue

Comments

@Wim-De-Clercq
Copy link

Wim-De-Clercq commented Sep 3, 2019 •
edited
Loading

I'm coming from here: pyproj4/pyproj#422 (comment)

The proj.4 string at https://epsg.io/31370 is

+proj=lcc +lat_1=51.16666723333333 +lat_2=49.8333339 +lat_0=90 +lon_0=4.367486666666666 +x_0=150000.013 +y_0=5400088.438 +ellps=intl +towgs84=-106.869,52.2978,-103.724,0.3366,-0.457,1.8422,-1.2747 +units=m +no_defs

projinfo gives me

$ projinfo EPSG:31370
+proj=lcc +lat_0=90 +lon_0=4.36748666666667 +lat_1=51.1666672333333 +lat_2=49.8333339 +x_0=150000.013 +y_0=5400088.438 +ellps=intl +towgs84=-99.059,53.322,-112.486,0.419,-0.83,1.885,-1 +units=m +no_defs +type=crs

The towsg84 part is different.

Not sure if that's all that is wrong, please read the comment at pyproj4/pyproj#422 (comment) I'm not sure if that is also proj related or pyproj.
@kbevers
Copy link
Member

kbevers commented Sep 3, 2019

The data is in fact correct. The towgs84 parameter holds parameters for a datumshift to WGS84. PROJ currently knows two mappings from EPSG:31370 to WGS84, as seen with projinfo:

projinfo -s EPSG:31370 -t EPSG:4326
Candidate operations found: 2
-------------------------------------
Operation n┬░1:

unknown id, Inverse of Belgian Lambert 72 + BD72 to WGS 84 (1), 1 m, Belgium - onshore

PROJ string:
+proj=pipeline +step +inv +proj=lcc +lat_0=90 +lon_0=4.36748666666667 +lat_1=51.1666672333333 +lat_2=49.8333339 +x_0=150000.013 +y_0=5400088.438 +ellps=intl +step +proj=push +v_3 +step +proj=cart +ellps=intl +step +proj=helmert +x=-99.059 +y=53.322 +z=-112.486 +rx=-0.419 +ry=0.83 +rz=-1.885 +s=-1 +convention=coordinate_frame +step +inv +proj=cart +ellps=WGS84 +step +proj=pop +v_3 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1

WKT2_2018 string:
CONCATENATEDOPERATION["Inverse of Belgian Lambert 72 + BD72 to WGS 84 (1)",
    SOURCECRS[
        PROJCRS["Belge 1972 / Belgian Lambert 72",
            BASEGEOGCRS["Belge 1972",
                DATUM["Reseau National Belge 1972",
                    ELLIPSOID["International 1924",6378388,297,
                        LENGTHUNIT["metre",1]]],
                PRIMEM["Greenwich",0,
                    ANGLEUNIT["degree",0.0174532925199433]],
                ID["EPSG",4313]],
            CONVERSION["Belgian Lambert 72",
                METHOD["Lambert Conic Conformal (2SP)",
                    ID["EPSG",9802]],
                PARAMETER["Latitude of false origin",90,
                    ANGLEUNIT["degree",0.0174532925199433],
                    ID["EPSG",8821]],
                PARAMETER["Longitude of false origin",4.36748666666667,
                    ANGLEUNIT["degree",0.0174532925199433],
                    ID["EPSG",8822]],
                PARAMETER["Latitude of 1st standard parallel",51.1666672333333,
                    ANGLEUNIT["degree",0.0174532925199433],
                    ID["EPSG",8823]],
                PARAMETER["Latitude of 2nd standard parallel",49.8333339,
                    ANGLEUNIT["degree",0.0174532925199433],
                    ID["EPSG",8824]],
                PARAMETER["Easting at false origin",150000.013,
                    LENGTHUNIT["metre",1],
                    ID["EPSG",8826]],
                PARAMETER["Northing at false origin",5400088.438,
                    LENGTHUNIT["metre",1],
                    ID["EPSG",8827]]],
            CS[Cartesian,2],
                AXIS["easting (X)",east,
                    ORDER[1],
                    LENGTHUNIT["metre",1]],
                AXIS["northing (Y)",north,
                    ORDER[2],
                    LENGTHUNIT["metre",1]],
            ID["EPSG",31370]]],
    TARGETCRS[
        GEOGCRS["WGS 84",
            DATUM["World Geodetic System 1984",
                ELLIPSOID["WGS 84",6378137,298.257223563,
                    LENGTHUNIT["metre",1]]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["degree",0.0174532925199433]],
            CS[ellipsoidal,2],
                AXIS["geodetic latitude (Lat)",north,
                    ORDER[1],
                    ANGLEUNIT["degree",0.0174532925199433]],
                AXIS["geodetic longitude (Lon)",east,
                    ORDER[2],
                    ANGLEUNIT["degree",0.0174532925199433]],
            ID["EPSG",4326]]],
    STEP[
        CONVERSION["Inverse of Belgian Lambert 72",
            METHOD["Inverse of Lambert Conic Conformal (2SP)",
                ID["INVERSE(EPSG)",9802]],
            PARAMETER["Latitude of false origin",90,
                ANGLEUNIT["degree",0.0174532925199433],
                ID["EPSG",8821]],
            PARAMETER["Longitude of false origin",4.36748666666667,
                ANGLEUNIT["degree",0.0174532925199433],
                ID["EPSG",8822]],
            PARAMETER["Latitude of 1st standard parallel",51.1666672333333,
                ANGLEUNIT["degree",0.0174532925199433],
                ID["EPSG",8823]],
            PARAMETER["Latitude of 2nd standard parallel",49.8333339,
                ANGLEUNIT["degree",0.0174532925199433],
                ID["EPSG",8824]],
            PARAMETER["Easting at false origin",150000.013,
                LENGTHUNIT["metre",1],
                ID["EPSG",8826]],
            PARAMETER["Northing at false origin",5400088.438,
                LENGTHUNIT["metre",1],
                ID["EPSG",8827]],
            ID["INVERSE(EPSG)",19961]]],
    STEP[
        COORDINATEOPERATION["BD72 to WGS 84 (1)",
            VERSION["IGN-Bel 7"],
            SOURCECRS[
                GEOGCRS["Belge 1972",
                    DATUM["Reseau National Belge 1972",
                        ELLIPSOID["International 1924",6378388,297,
                            LENGTHUNIT["metre",1]]],
                    PRIMEM["Greenwich",0,
                        ANGLEUNIT["degree",0.0174532925199433]],
                    CS[ellipsoidal,2],
                        AXIS["geodetic latitude (Lat)",north,
                            ORDER[1],
                            ANGLEUNIT["degree",0.0174532925199433]],
                        AXIS["geodetic longitude (Lon)",east,
                            ORDER[2],
                            ANGLEUNIT["degree",0.0174532925199433]],
                    ID["EPSG",4313]]],
            TARGETCRS[
                GEOGCRS["WGS 84",
                    DATUM["World Geodetic System 1984",
                        ELLIPSOID["WGS 84",6378137,298.257223563,
                            LENGTHUNIT["metre",1]]],
                    PRIMEM["Greenwich",0,
                        ANGLEUNIT["degree",0.0174532925199433]],
                    CS[ellipsoidal,2],
                        AXIS["geodetic latitude (Lat)",north,
                            ORDER[1],
                            ANGLEUNIT["degree",0.0174532925199433]],
                        AXIS["geodetic longitude (Lon)",east,
                            ORDER[2],
                            ANGLEUNIT["degree",0.0174532925199433]],
                    ID["EPSG",4326]]],
            METHOD["Coordinate Frame rotation (geog2D domain)",
                ID["EPSG",9607]],
            PARAMETER["X-axis translation",-99.059,
                LENGTHUNIT["metre",1],
                ID["EPSG",8605]],
            PARAMETER["Y-axis translation",53.322,
                LENGTHUNIT["metre",1],
                ID["EPSG",8606]],
            PARAMETER["Z-axis translation",-112.486,
                LENGTHUNIT["metre",1],
                ID["EPSG",8607]],
            PARAMETER["X-axis rotation",-0.419,
                ANGLEUNIT["arc-second",4.84813681109536E-06],
                ID["EPSG",8608]],
            PARAMETER["Y-axis rotation",0.83,
                ANGLEUNIT["arc-second",4.84813681109536E-06],
                ID["EPSG",8609]],
            PARAMETER["Z-axis rotation",-1.885,
                ANGLEUNIT["arc-second",4.84813681109536E-06],
                ID["EPSG",8610]],
            PARAMETER["Scale difference",-1,
                SCALEUNIT["parts per million",1E-06],
                ID["EPSG",8611]],
            OPERATIONACCURACY[1.0],
            ID["EPSG",1609]]],
    USAGE[
        SCOPE["unknown"],
        AREA["Belgium - onshore"],
        BBOX[49.5,2.5,51.51,6.4]]]

-------------------------------------
Operation n┬░2:

unknown id, Inverse of Belgian Lambert 72 + BD72 to WGS 84 (3), 1 m, Belgium - onshore

PROJ string:
+proj=pipeline +step +inv +proj=lcc +lat_0=90 +lon_0=4.36748666666667 +lat_1=51.1666672333333 +lat_2=49.8333339 +x_0=150000.013 +y_0=5400088.438 +ellps=intl +step +proj=push +v_3 +step +proj=cart +ellps=intl +step +proj=helmert +x=-106.8686 +y=52.2978 +z=-103.7239 +rx=-0.3366 +ry=0.457 +rz=-1.8422 +s=-1.2747 +convention=coordinate_frame +step +inv +proj=cart +ellps=WGS84 +step +proj=pop +v_3 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1

WKT2_2018 string:
CONCATENATEDOPERATION["Inverse of Belgian Lambert 72 + BD72 to WGS 84 (3)",
    SOURCECRS[
        PROJCRS["Belge 1972 / Belgian Lambert 72",
            BASEGEOGCRS["Belge 1972",
                DATUM["Reseau National Belge 1972",
                    ELLIPSOID["International 1924",6378388,297,
                        LENGTHUNIT["metre",1]]],
                PRIMEM["Greenwich",0,
                    ANGLEUNIT["degree",0.0174532925199433]],
                ID["EPSG",4313]],
            CONVERSION["Belgian Lambert 72",
                METHOD["Lambert Conic Conformal (2SP)",
                    ID["EPSG",9802]],
                PARAMETER["Latitude of false origin",90,
                    ANGLEUNIT["degree",0.0174532925199433],
                    ID["EPSG",8821]],
                PARAMETER["Longitude of false origin",4.36748666666667,
                    ANGLEUNIT["degree",0.0174532925199433],
                    ID["EPSG",8822]],
                PARAMETER["Latitude of 1st standard parallel",51.1666672333333,
                    ANGLEUNIT["degree",0.0174532925199433],
                    ID["EPSG",8823]],
                PARAMETER["Latitude of 2nd standard parallel",49.8333339,
                    ANGLEUNIT["degree",0.0174532925199433],
                    ID["EPSG",8824]],
                PARAMETER["Easting at false origin",150000.013,
                    LENGTHUNIT["metre",1],
                    ID["EPSG",8826]],
                PARAMETER["Northing at false origin",5400088.438,
                    LENGTHUNIT["metre",1],
                    ID["EPSG",8827]]],
            CS[Cartesian,2],
                AXIS["easting (X)",east,
                    ORDER[1],
                    LENGTHUNIT["metre",1]],
                AXIS["northing (Y)",north,
                    ORDER[2],
                    LENGTHUNIT["metre",1]],
            ID["EPSG",31370]]],
    TARGETCRS[
        GEOGCRS["WGS 84",
            DATUM["World Geodetic System 1984",
                ELLIPSOID["WGS 84",6378137,298.257223563,
                    LENGTHUNIT["metre",1]]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["degree",0.0174532925199433]],
            CS[ellipsoidal,2],
                AXIS["geodetic latitude (Lat)",north,
                    ORDER[1],
                    ANGLEUNIT["degree",0.0174532925199433]],
                AXIS["geodetic longitude (Lon)",east,
                    ORDER[2],
                    ANGLEUNIT["degree",0.0174532925199433]],
            ID["EPSG",4326]]],
    STEP[
        CONVERSION["Inverse of Belgian Lambert 72",
            METHOD["Inverse of Lambert Conic Conformal (2SP)",
                ID["INVERSE(EPSG)",9802]],
            PARAMETER["Latitude of false origin",90,
                ANGLEUNIT["degree",0.0174532925199433],
                ID["EPSG",8821]],
            PARAMETER["Longitude of false origin",4.36748666666667,
                ANGLEUNIT["degree",0.0174532925199433],
                ID["EPSG",8822]],
            PARAMETER["Latitude of 1st standard parallel",51.1666672333333,
                ANGLEUNIT["degree",0.0174532925199433],
                ID["EPSG",8823]],
            PARAMETER["Latitude of 2nd standard parallel",49.8333339,
                ANGLEUNIT["degree",0.0174532925199433],
                ID["EPSG",8824]],
            PARAMETER["Easting at false origin",150000.013,
                LENGTHUNIT["metre",1],
                ID["EPSG",8826]],
            PARAMETER["Northing at false origin",5400088.438,
                LENGTHUNIT["metre",1],
                ID["EPSG",8827]],
            ID["INVERSE(EPSG)",19961]]],
    STEP[
        COORDINATEOPERATION["BD72 to WGS 84 (3)",
            VERSION["IGN-Bel 1m"],
            SOURCECRS[
                GEOGCRS["Belge 1972",
                    DATUM["Reseau National Belge 1972",
                        ELLIPSOID["International 1924",6378388,297,
                            LENGTHUNIT["metre",1]]],
                    PRIMEM["Greenwich",0,
                        ANGLEUNIT["degree",0.0174532925199433]],
                    CS[ellipsoidal,2],
                        AXIS["geodetic latitude (Lat)",north,
                            ORDER[1],
                            ANGLEUNIT["degree",0.0174532925199433]],
                        AXIS["geodetic longitude (Lon)",east,
                            ORDER[2],
                            ANGLEUNIT["degree",0.0174532925199433]],
                    ID["EPSG",4313]]],
            TARGETCRS[
                GEOGCRS["WGS 84",
                    DATUM["World Geodetic System 1984",
                        ELLIPSOID["WGS 84",6378137,298.257223563,
                            LENGTHUNIT["metre",1]]],
                    PRIMEM["Greenwich",0,
                        ANGLEUNIT["degree",0.0174532925199433]],
                    CS[ellipsoidal,2],
                        AXIS["geodetic latitude (Lat)",north,
                            ORDER[1],
                            ANGLEUNIT["degree",0.0174532925199433]],
                        AXIS["geodetic longitude (Lon)",east,
                            ORDER[2],
                            ANGLEUNIT["degree",0.0174532925199433]],
                    ID["EPSG",4326]]],
            METHOD["Coordinate Frame rotation (geog2D domain)",
                ID["EPSG",9607]],
            PARAMETER["X-axis translation",-106.8686,
                LENGTHUNIT["metre",1],
                ID["EPSG",8605]],
            PARAMETER["Y-axis translation",52.2978,
                LENGTHUNIT["metre",1],
                ID["EPSG",8606]],
            PARAMETER["Z-axis translation",-103.7239,
                LENGTHUNIT["metre",1],
                ID["EPSG",8607]],
            PARAMETER["X-axis rotation",-0.3366,
                ANGLEUNIT["arc-second",4.84813681109536E-06],
                ID["EPSG",8608]],
            PARAMETER["Y-axis rotation",0.457,
                ANGLEUNIT["arc-second",4.84813681109536E-06],
                ID["EPSG",8609]],
            PARAMETER["Z-axis rotation",-1.8422,
                ANGLEUNIT["arc-second",4.84813681109536E-06],
                ID["EPSG",8610]],
            PARAMETER["Scale difference",-1.2747,
                SCALEUNIT["parts per million",1E-06],
                ID["EPSG",8611]],
            OPERATIONACCURACY[1.0],
            ID["EPSG",15929]]],
    USAGE[
        SCOPE["unknown"],
        AREA["Belgium - onshore"],
        BBOX[49.5,2.5,51.51,6.4]]]

As you can see, the first operation that is listed is holds the parameters you claim to be incorrect. The second operation holds the "correct" parameters. In reality both are correct. From the data that is in the PROJ database it is difficult to determine which of the two is best for your use case. They have the same accuracy and area of use.

This is a limitation in the old way of describing a CRS in terms of WGS84: There can be several mappings. The current recommendation is to noget use proj-strings to define CRS's. Instead you should use the corresponding EPSG-code.

@Wim-De-Clercq Wim-De-Clercq mentioned this issue Sep 3, 2019
Closed
@Wim-De-Clercq
Copy link
Author

But, whichever transformation it would use, both should yield the same results, right?

@busstoptaktik
Copy link
Member

busstoptaktik commented Sep 3, 2019 via email •
edited
Loading

@kbevers
Copy link
Member

kbevers commented Sep 4, 2019

But, whichever transformation it would use, both should yield the same results, right?

Let's try it out:

First transformation:

echo 50 4 0 0 | cct +proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=push +v_3 +step +proj=cart +ellps=WGS84 +step +inv +proj=helmert +x=-99.059 +y=53.322 +z=-112.486 +rx=-0.419 +ry=0.83 +rz=-1.885 +s=-1 +convention=coordinate_frame +step +inv +proj=cart +ellps=intl +step +proj=pop +v_3 +step +proj=lcc +lat_0=90 +lon_0=4.36748666666667 +lat_1=51.1666672333333 +lat_2=49.8333339 +x_0=150000.013 +y_0=5400088.438 +ellps=intl
  123563.5666     76585.0785        0.0000           0.0000

Second transformation:

echo 50 4 0 0 | cct +proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=push +v_3 +step +proj=cart +ellps=WGS84 +step +inv +proj=helmert +x=-106.8686 +y=52.2978 +z=-103.7239 +rx=-0.3366 +ry=0.457 +rz=-1.8422 +s=-1.2747 +convention=coordinate_frame +step +inv +proj=cart +ellps=intl +step +proj=pop +v_3 +step +proj=lcc +lat_0=90 +lon_0=4.36748666666667 +lat_1=51.1666672333333 +lat_2=49.8333339 +x_0=150000.013 +y_0=5400088.438 +ellps=intl
  123563.5713     76585.0788        0.0000        0.0000

In summary, yes, it yields the same result (remember that the accuracy of both transformations is 1 m). Projstrings from days past can't be expected to be the default in PROJ 6+, since the internal plumbing has been changed quite a bit, but generally you can trust PROJ to do the right thing in circumstances like this.

I'm closing this since it is not a bug.

@kbevers kbevers closed this as completed Sep 4, 2019
@rouault
Copy link
Member

rouault commented Sep 4, 2019

I had an exchange with Roger Lott from IOGP about that issue, and he gave very interested hints.

My email to him

Roger,

I'm wondering if there's some logic of preference in the EPSG 
epsg_coordoperation table when transforming between 2 CRS, when there are 
several transformations which are reported to have the same accuracy, area of 
use and scope (after removing the deprecated entries)

For example to transform between BD72 and WGS84, if I keep the ones with the 
best accuracy, I get:

EPSG:1609, BD72 to WGS 84 (1), 1m, IGN-Bel 7, For applications to an accuracy 
of 1 metre, revision_date = 2010/03/30
EPSG:15929,BD72 to WGS 84 (3), 1m, IGN-Bel 1m, For applications to an accuracy 
of 1 metre, revision_date = 2018/01/19

Both are not superseded by anything. For the average user who would have no 
clue of which one to use, would there be one "better" than the other one ?

I'm wondering in particular if the revision_date shouldn't be taken into 
account when all other criteria have been considered, that is take the entry 
with the most recent revision date ? or perhaps better, taking the 
transformation which has the highest number in parenthesis in its name, since 
it is likely to have been added after previous ones ?

It happens that PROJ 6 will select EPSG:1609 (it arbitrarily fallbacks to 
ASCII ordering on the transformation name, and using the entry with the 
'lowest' sorting order), whereas older PROJ versions selected EPSG:15929 (in a 
different arbitrary choice). This question is triggered by https://github.com/
OSGeo/PROJ/issues/1580

Hum, perhaps answering myself, I see the following entry in epsg_change:

OGRFeature(epsg_change):1998
  change_id (Real) = 2017.055
  report_date (Date) = 2017/12/14
  date_closed (Date) = 2018/01/19
  reporter (String) = Melita Kennedy; for ESRI Belgium
  request (String) = Add Belgium transformation data
  tables_affected (String) = Coordinate Operation
  codes_affected (String) = 15929
  change_comment (String) = 
  action (String) = Added transformation 8369. For transformation 15929, 
amended version, scope and accuracy from 0.5m to 1m to compensate for tectonic 
plate motion since definition.

So in the past EPSG:15929 had an acuracy of 0.5m, which probably explains why 
it was picked up by past PROJ versions. So for the WGS84 of past years, EPSG:
15929 was probably the best choice, but for WGS84 as of today, not that much.

Side question: what is the point of having also "EPSG:1610, BD72 to WGS 84 
(2), IGN-Bel 3, For applications to an accuracy of 5 metres" in the database, 
which also has the same area of use than 1609 and 15959 ? Are there use cases 
where a user would prefer using that one instead of 1609 or 15929 ?

Thanks,

Even

His answer

Hi Even

The problem is caused by EPSG 'copy' transformations policy. Although we have always recommended that a direct transformation be used in preference to a hub, we recognise the validity of the hub concept if it is done properly, and so (i) have documented it in our Guidance Note 7-1, and (ii) recognising that WGS 84 is the de facto choice as hub CRS, have made 'EPSG copies' to WGS 84 of CTs that are to another CRS (such as, but not limited to, ETRS89). There are two cases for these 'EPSG copies':

i)      The local CRS is defined as being coincident with ITRF at some chosen epoch. As WGS 84 can be assumed to be coincident with ITRF, then it may be assumed that the local CRS is coincident with WGS 84 at the chosen epoch. An example is ETRS89, equal to ITRF at epoch 1989.00. In this case EPSG will create a 'null' transformation (geocentric translations method, with tX=tY=tZ=0.0m). For example EPSG::1149.

ii)     A local CRS is related to ITRF or to another local CRS defined as being coincident with ITRF at some chosen epoch. An example id BD72 in Belgium, which has CTs to ETRS89. In this case, EPSG copies the BD72 to ETRS89 transformation as a BD72 to ETRS89 transformation. More on these below.

Now, in the data model all transformations have an accuracy attribute. It is there to recognise that transformations are empirically derived, will therefore inherit error from the coordinates through which they were derived, which in turn come from errors inherent in the observations measurements. A transformation accuracy of 0 is errorless. From memory, the highest transformation accuracy in the EPSG Dataset is about 45m, representing a very poor transformation.

When EPSG creates these 'EPSG copy' CTs, we increase the transformation accuracy to allow for the uncertainty in the assumption of equality of local CRS to WGS 84. Usually the biggest contribution to this is that the local CRS is static and WGS 84 is dynamic, so they drift apart by the tectonic plate movement since the epoch of coincidence. EPSG policy has been to assume neglecting tectonic plate movement of up to 1m as being acceptable, and assign to the EPSG copy CT an accuracy of 1m when the accuracy of the source CT is under 1m (examples below), or the source CT accuracy if it is over 1m. However last year we changed the accuracy of a GDA94 to WGS 84 EPSG null CT for Australia (EPSG::1150) from 1m to 3m to account for more tectonic plate movement than we had previously allowed for (so there are other examples to that you cite for Belgium of the accuracy being changed).

So in the Belgian case we have 6 transformations:

BD72 to ETRS89                                                                     BD72 to WGS 84

(no equivalent)                                                                         EPSG::1610, source IGN, accuracy 5m.

EPSG::1652, source IGN via EuroGeographics, accuracy 1m.    EPSG copy = EPSG::1609, source IOGP, accuracy made 1m as 1652 =< 1m

EPSG::15928, source IGN, accuracy 0.2m.                                EPSG copy = EPSG::15929, source IOGP, accuracy made 1m as 15928 =< 1m

EPSG::8369, source IGN, accuracy 0.01m.                                (no equivalent)

We do not make copies to WGS 84 of high accuracy CTs if there is already at least one CT to WGS 84 at the 1m level, so there is no EPSG copy of EPSG::8369 because EPSG::15929 has been copied from EPSG::15928.

The reason that EPSG::1609 and EPSG::15929 both have an accuracy of 1m is that they are both copies from CTs with accuracy of 1m or less.

This exposes a flaw in our transformation accuracy assignment policy. Based on their transformation accuracy alone, EPSG::15928 (0.2m) is clearly better than EPSG::1652 (1m) so one would assume that EPSG::15929 is preferable to EPSG::1609, but our assigned transformation accuracy does not expose that. The difference between the two CTs to ETRS89 is within the 1m tectonic movement uncertainty we have applied. A large number of copy CTs to WGS 84 have been allocated a transformation accuracy of exactly 1m. IOGP needs to review this.

Meanwhile, given the existing situation of transformation accuracy equality, what might you do? There are clues, but not easily machine-implementable because they are mainly in free-text remarks.

a)      All EPSG copy CTs will have remarks along the lines of (from 15929) "Parameter values from BD72 to ETRS89 (2) (code 15928) assuming that ETRS89 is equivalent to WGS 84 within the accuracy of the transformation". When you find two like this with equal accuracy, examining the transformation accuracy of the source CRSs will likely suggest a way forward (as in previous paragraph with EPSG::1652 and EPSG::15928).

        The information source (sic, not data source) of the transformation to WGS 84 being IOGP (or OGP for older records) is a likely indication that the CT is a copy and the originals could be used to differentiate.

b)      The transformation version field sometimes (but definitely not always) contains an indication of the transformation accuracy.

c)      The transformation method used might be considered an indication of transformation accuracy. A 7-parameter Helmert similarity transformation would be expected to be better than a 3-parameter geocentric translation similarity transformation (because it better models CRS offsets by including orientation and scale differences). And a grid interpolation would be expected to be better than a 7-parameter transformation (because it better models CRS offsets by including local variations). But a 10-parameter transformation will not be any more accurate than a 7-parameter one: its derivation will be better, but once derived in application there is no difference. Polynomial methods could be assumed to be better than similarity methods but not as good as grid interpolation. But all rather speculative, and likely to already be factored into the transformation accuracy value.

d)      It could be assumed that a later CT is more accurate than an earlier one, as in your question. Variant value can be used for determining that. Also date of entry (earliest change request value) is a possible source. Although EPSG codes are now assigned sequentially, this has not always been the case so again EPSG code is no indication of data entry. But there definitely are cases of transformations of lesser accuracy being loaded after higher accuracy, so I would not advocate this path using any or all of these indicators.

Of the above, (a) is by some way the best approach.

>    Side question: what is the point of having also "EPSG:1610, BD72 to WGS 84
>    (2), IGN-Bel 3, For applications to an accuracy of 5 metres" in the database,
>      which also has the same area of use than 1609 and 15959 ? Are there use cases
>     where a user would prefer using that one instead of 1609 or 15929 ?

It may be needed when the user is using applications that do not implement the methods used by higher accuracy transformations. For example, most hand-held GPS receivers only implement 3 translations (not Helmert 7-parameter and certainly not grid interpolation such as NTv2). So users will have to use that simple method. For Belgium their choice would be use EPSG::1610 or derive their own for their project area from one of the more accurate ones (especially if there is a grid interpolation method available, such as EPSG::8369, I usually use that to derive a 3-parameter transformation applicable to my project area).

Another use case is that the low accuracy transformation was used before the higher accuracy transformations became available, and for consistency rather than absolute accuracy the use wants to continue to use the same transformation as he began with.

Feedback requested

FYI, at its meeting in October the IOGP Geodesy Subcommittee will be reviewing its policy on these copy transformations. Internally we have raised the question about whether we should be doing it, and if doing it is ok then is our actual way of doing it appropriate. What probably should have been happening is that because WGS 84 is dynamic, CTs to it should have been using time-dependent methods. And conflating two different bits of data – the error induced by applying a transformation and an allowance for tectonic plate movement – cannot be considered good practice. So we need to make decisions about our future data population strategy, and once that is made decisions on whether or not the same criteria should be retrospectively applied to existing entries. Any input you or your colleagues have on this will be taken into consideration (but not necessarily adopted if conflicting with a bigger picture).

Regards,

Roger

@rouault
Copy link
Member

rouault commented Sep 4, 2019

Re-opening. Given Roger Lott's hints and following discussion, it seems that the EPSG dataset might be amended in a later version to tweak the accuracy of EPSG:15929 and EPSG:1609 so that the former is selected.

@rouault rouault reopened this Sep 4, 2019
@rouault rouault changed the title Incorrect data for EPSG: 31370 Incorrect data for EPSG:31370 (Belgian Lambert 72) / non-optimal transformation chosen for BD72 to WGS84 Sep 11, 2019
@stale
Copy link

stale bot commented Jan 8, 2020

This issue has been automatically marked as stale because it has not had recent activity. It will be closed if no further activity occurs. Thank you for your contributions.

@stale stale bot added the wontfix label Jan 8, 2020
@kbevers kbevers added the pinned Prevent stale-bot from closing an issue label Jan 8, 2020
@stale stale bot removed the wontfix label Jan 8, 2020
@SBylemans SBylemans mentioned this issue Sep 30, 2021
Merged
@rouault
Copy link
Member

rouault commented Feb 18, 2023

A "recent" change has favored,all other things equal, transformation names that are lexicographically bigger than other ones, so now BD72 to WGS 84 (3) / EPSG:15929 is the prefered one

@rouault rouault closed this as completed Feb 18, 2023
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
pinned Prevent stale-bot from closing an issue
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
4 participants
@rouault @Wim-De-Clercq @busstoptaktik @kbevers