createOperations(): tune geod_crs_of_datum_ensemble <--> geod_crs_of_another_datum

[rouault]

This reverts largely #4736 ("make ETRS89 great again"), except the test cases it fixed and added. The fragile logic of #4736 is replaced by an hopefully less fragile one of commit 400d5c5

createOperations(): tune geod_crs_of_datum_ensemble <--> geod_crs_of_another_datum

typically ETRS89 to some national CRS (e.g. Pulkovo 42(58), etc.)

If the existing lookup paths have not tried to go through an
intermediate datum already (which can be the case if direct
transformations exists), then force trying with the members of the datum
ensemble.

e.g. ETRS89 to Pulkovo 42(58) has direct transformations (EPSG:15993
Romania 10m and EPSG:1644 Poland 1m), but there's also a ETRS89-ROU to
Pulkovo 42(58), EPSG:15994, for Romania 3m.

This has been tuned quite a bit to obtain the same runtime when running ctest as in master (earlier attempts resulted in terrible runtime in some cases)

This implements this suggestion of "ETRS89 datum change issues v3 2026-05-08":

The software logic for finding a transformation between two specified CRSs where one has a datum ensemble is: 
    a) list all transformations that are an exact match of source and target CRS, or have source and target transposed and the CT method is reversible [collectively (i) above], and 
    b) find additional transformations that may also be used [(ii) above] by using the datum ensemble.
    • Check whether the source or target CRS is based on a datum that has the conventional RS set (and therefore is a member of a datum ensemble);
    • Find transformations associated with CRSs based on members of the datum ensemble and the second of the two specified CRSs; 
    • Test whether that datum ensemble is the same as that for the ensemble-based specified CRS;
    • Test that the extent and accuracy of the transformation meet requirements.

A difficulty with using datum ensembles is that their use requires data mining. In some circumstances, using name and/or alias may produce results more quickly, but because of issues with uncertainty in data population and text matching these results will be less complete and reliable than from using the datum ensemble. 

I confirm the difficulty of the data mining!!!

[rouault]

@kbevers Note that commit 'createOperationsEnsembleCRSToOtherGeodCRS(): better deal with geographic 3D / geocentric CRS' has a consequence on 1 test case of interest for you

Now projinfo -s EPSG:7789 -t EPSG:4936 --area EPSG:1080 --summary --hide-ballpark (that is from ITRF2014 geocentric to ETRS89 geocentric, using Denmark area of use) returns

Candidate operations found: 2
    EPSG:10894, ITRF2014 to ETRS89-DNK (1), 0.006 m, Denmark - onshore and offshore., time-dependent operation
    NKG:ITRF2014_TO_DK, ITRF2014 to ETRS89(DK), 0.01 m, Denmark - onshore and offshore., time-dependent operation

Previously it only returned the NKG:ITRF2014_TO_DK operation.

[rouault]

@nich0lasbr0wn (hoping you're still on topics related to Australian coordinate transformation), also CC'ing @nyalldawson

This pull requests changes the way PROJ infers coordinate operations when a datum ensemble is involved in the source or target CRS. This is based on a ongoing discussion with EPSG, following recent changes in ETRS89. It enhances the traditional lookups that PROJ used to do by also considering transformations between a member of the datum ensemble (a realization of ETRS89 or WGS 84 typically) and the other CRS. A typical exemple in Europe would be to use "ETRS89-ROU to Pulkovo 42(58)" when the user asked for generic ETRS89 to "Pulkovo 42(58)" (the later being what users do since ETRS89-ROU is a recent addition. This PR is generic for all datum ensembles, not just ETRS89, and so WGS 84 'generic' (EPSG:4326/4978/4978) is also concerned. Which brings us to Australia...

For now I've blocklisted the GDA94 <--> WGS 84 and GDA2020 <--> WGS 84 cases as I'm unusure of what is desired. But if I remove the blocking logic,

projinfo GDA94 "WGS 84" --summary --spatial-test intersects --grid-check none
would now return:

Candidate operations found: 4
EPSG:1150, GDA94 to WGS 84 (1), 3.0 m, Australia including Lord Howe Island, Macquarie Island, Ashmore and Cartier Islands, Christmas Island, Cocos (Keeling) Islands, Norfolk Island. All onshore and offshore.
EPSG:9688, GDA94 to WGS 84 (2), 3.0 m, Australia including Lord Howe Island, Macquarie Island, Ashmore and Cartier Islands, Christmas Island, Cocos (Keeling) Islands, Norfolk Island. All onshore and offshore.
EPSG:9687, GDA94 to WGS 84 (G1762) (2), 0.25 m, Australia - Australian Capital Territory; New South Wales; Northern Territory; Queensland; South Australia; Tasmania; Western Australia; Victoria., at least one grid missing, time-dependent operation
DERIVED_FROM(EPSG):9689, GDA94 to WGS 84 (3), 3.0 m, Australia - Australian Capital Territory; New South Wales; Northern Territory; Queensland; South Australia; Tasmania; Western Australia; Victoria., at least one grid missing

The novelty (if the unblocking would be removed) is the addition of the "EPSG:9687, GDA94 to WGS 84 (G1762) (2)" operation that was previously not used because targetting specifically WGS 84 (G1762) is now considered. And as its accuracy is 0.25 compared to the other ones, it would be used (when the GDA94_GDA2020_conformal_and_distortion grid is available). Whereas currently the "GDA94 to WGS 84 (1)" no-op is selected.

Similarly, for GDA2020 to WGS 84 (but only for 3D operations),

$ bin/projinfo GDA2020 "WGS 84" --3d  --spatial-test intersects --grid-check none --summary

returns

Candidate operations found: 4
unknown id, Conversion from GDA2020 (geog3D) to GDA2020 (geocentric) + GDA2020 to WGS 84 (G1762) (1) + Conversion from WGS 84 (G1762) (geocentric) to WGS 84 (G1762) (geog3D) + Inverse of WGS 84 to WGS 84 (G1762), 2.2 m, Australia including Lord Howe Island, Macquarie Island, Ashmore and Cartier Islands, Christmas Island, Cocos (Keeling) Islands, Norfolk Island. All onshore and offshore., time-dependent operation
DERIVED_FROM(EPSG):8450, GDA2020 to WGS 84 (2), 3.0 m, Australia including Lord Howe Island, Macquarie Island, Ashmore and Cartier Islands, Christmas Island, Cocos (Keeling) Islands, Norfolk Island. All onshore and offshore.
INVERSE(DERIVED_FROM(EPSG)):9690, Inverse of WGS 84 to GDA2020 (3), 3.0 m, Australia including Lord Howe Island, Macquarie Island, Ashmore and Cartier Islands, Christmas Island, Cocos (Keeling) Islands, Norfolk Island. All onshore and offshore.
INVERSE(DERIVED_FROM(EPSG)):9691, Inverse of WGS 84 to GDA2020 (4), 3.0 m, Australia - Australian Capital Territory; New South Wales; Northern Territory; Queensland; South Australia; Tasmania; Western Australia; Victoria., at least one grid missing

The novelty (if the unblocking would be removed) would be the first operation "GDA2020 to WGS 84 (G1762) (1)" that would be used, whereas currently the null transformation "GDA2020 to WGS 84 (2)" is used.

Would such a change be desired ? I'd note that PROJ logic can use time-dependent operations, even when the input coordinate tuple has no epoch (in which case the time-dependent part is not applied). So for "GDA2020 to WGS 84 (G1762) (1) " that would still be essentially a no-op (just the tiny difference of semi-minor axis between the GRS80 and WGS 84 ellipsoids causing sub-millimetric shift)

[kbevers]

Note that commit 'createOperationsEnsembleCRSToOtherGeodCRS(): better deal with geographic 3D / geocentric CRS' has a consequence on 1 test case of interest for you

This is good, thanks for testing it. The two transformations returned are the same so probably I should mark the NKG ones deprecated.

[rouault]

The two transformations returned are the same

not stritcly the same PROJ string, but apparently equivalent ITRF2014<-->ETRF2014 Helmert using different centrals epochs. On a Copenhagen test point, both indeed match perfectly:

echo "55.6761 12.5683 0 2026" | bin/cs2cs ITRF2014 EPSG:7789 --3d | PROJ_NETWORK=ON PROJ_DATA=data  bin/cct -d 8 "ITRF2014 to ETRS89(DK)"
3518305.34403383  784390.17421352  5244191.46952714     2026.0000

echo "55.6761 12.5683 0 2026" | bin/cs2cs ITRF2014 EPSG:7789 --3d | PROJ_NETWORK=ON PROJ_DATA=data bin/cct -d 8 "ITRF2014 to ETRS89-DNK (1)"
3518305.34403383  784390.17421352  5244191.46952714     2026.0000

[RogerLott]

@kbevers, have you checked whether the commit works for Danish or Greenlandic vertical ensembles?

[kbevers]

@kbevers, have you checked whether the commit works for Danish or Greenlandic vertical ensembles?

@RogerLott I've poked around a bit and haven't found any issues compared to a previous PROJ version

[rouault]

@kbevers @jjimenezshaw I believe this can be merged independently whether EPSG reverts a number of the changes as discussed in today's meeting or keep them, right ?

[kbevers]

I would say so, yes

[jjimenezshaw]

@kbevers @jjimenezshaw I believe this can be merged independently whether EPSG reverts a number of the changes as discussed in today's meeting or keep them, right ?

Yes, that's my understanding