OSGeo / PROJ - TransformerGroup Ordinality

Question

My question is within the context of working with pyproj, the popular python package wrapping OSGeo's PROJ library, but it probably spills over into coordinate transformations as a whole and transformation grids.

Pyproj allows for constructing TransfomerGroups, which allow enumerating *each of the possible transformations between two coordinate reference systems.

When a TransformerGroup is created, PROJ is able to identify whether an available transformer is the best_available TransformerGroup. I would suspect that in each case, the only way that there would be a better transformation available, is whether the library could make use of a particular Transformation grid in that area.

Take the code below, I've set up two TransformerGroups, not downloaded any transformation grids. From the transformations below, in one case, we have better transformers available, and in another case, we do not.

# python                                                                                                              
Python 3.12.1 (tags/v3.12.1:2305ca5, Dec  7 2023, 22:03:25) [MSC v.1937 64 bit (AMD64)] on win32                      
Type "help", "copyright", "credits" or "license" for more information.                                                
>>> import pyproj                                                                                                     
>>> from pyproj.transformer import TransformerGroup                                                                   
>>> tg = TransformerGroup(4326, 26917)                                                                                
>>> tg                                                                                                                
<TransformerGroup: best_available=True>                                                                               
- transformers: 1                                                                                                     
- unavailable_operations: 14                                                                                          
>>> tg2 = TransformerGroup(4326, 32155)                                                                               
C:\Users\Michael.DiFelice\AppData\Local\Programs\Python\Python312\Lib\site-packages\pyproj\transformer.py:207: UserWar
ning: Best transformation is not available due to missing Grid(short_name=us_noaa_wyhpgn.tif, full_name=, package_name
=, url=https://cdn.proj.org/us_noaa_wyhpgn.tif, direct_download=True, open_license=True, available=False)             
  super().__init__(                                                                                                   
>>> tg2                                                                                                               
<TransformerGroup: best_available=False>                                                                              
- transformers: 1                                                                                                     
- unavailable_operations: 5                                                                                           
>>> [t.name for t in tg2.unavailable_operations]                                                                      
['Inverse of NAD83 to WGS 84 (44) + SPCS83 Wyoming East zone (meters)', 'Inverse of NAD83 to WGS 84 (19) + SPCS83 Wyom
ing East zone (meters)', 'Inverse of NAD83 to WGS 84 (36) + SPCS83 Wyoming East zone (meters)', 'Inverse of NAD83 to W
GS 84 (16) + SPCS83 Wyoming East zone (meters)', 'Inverse of NAD83 to WGS 84 (28) + SPCS83 Wyoming East zone (meters)'
]                                                                                                                     

How do we order transformations?

What makes one transformation better than another?

Having a best_available transformation seems to imply ordinality with transformations, but without it always being tied to some transformation grid, I can't seem to guess what would make one better than another.

I've done some digging into pyproj, and it seems as though the transformations are provided in an ordered list: https://github.com/pyproj4/pyproj/blob/7aab7d91e14230d846f74c418810652859d3c69d/pyproj/_transformer.pyx#L227-L258. best_available is just indicating whether the first transformation is available or not. It doesn't answer the question, however, as to how transformations are ordered.

Answer 1

Accuracy is not the whole story because for 4326=>26917 the best transformation is non-grid with accuracy=4.0 while all grid transformations are accuracy=1.5 | 2.0

Looks like the rules for sorting are as follow:

1. Operations where grids are all available go before other
2. Operations where grids are all known in our DB go before other
3. Operations with known accuracy go before those with unknown accuracy
4. Operations with larger non-zero area of use go before those with lower one
5. Operations with better accuracy go before those with worse one
6. The less intermediate steps, the better
7. Compare number of steps in PROJ pipeline, and prefer the ones with less operations.
8. "Ballpark geographic offset from NAD83(CSRS)v6 to NAD83(CSRS)" > "Ballpark geographic offset from ITRF2008 to NAD83(CSRS)"
9. Prefer shorter name
10. Something about NTF (Paris)

As of 9.5 the function that sorts the transformations is osgeo::proj::operation::SortFunction::compare(..) and this is where I came up with the criteria above.

Since area takes precedence over accuracy the non-grid Inverse of NAD83 to WGS 84 (1) + UTM zone 17N takes precedence over gridded Inverse of NAD83 to WGS 84 (6) + UTM zone 17N

name=Inverse of NAD83 to WGS 84 (1) + UTM zone 17N area=3.5449 accuracy=4 isPROJExportable=1 hasGrids=0 gridsAvailable=1 gridsKnown=1 stepCount=2 projStepCount=3 ballpark=0 vertBallpark=0 isNull=0
name=Inverse of NAD83 to WGS 84 (6) + UTM zone 17N area=0.334836 accuracy=1.5 isPROJExportable=1 hasGrids=1 gridsAvailable=1 gridsKnown=1 stepCount=2 projStepCount=4 ballpark=0 vertBallpark=0 isNull=0
name=Inverse of NAD83 to WGS 84 (18) + UTM zone 17N area=0.410755 accuracy=2 isPROJExportable=1 hasGrids=1 gridsAvailable=1 gridsKnown=1 stepCount=2 projStepCount=4 ballpark=0 vertBallpark=0 isNull=0
name=Inverse of NAD83 to WGS 84 (26) + UTM zone 17N area=0.15291 accuracy=2 isPROJExportable=1 hasGrids=1 gridsAvailable=1 gridsKnown=1 stepCount=2 projStepCount=4 ballpark=0 vertBallpark=0 isNull=0
name=Inverse of NAD83 to WGS 84 (31) + UTM zone 17N area=0.103198 accuracy=2 isPROJExportable=1 hasGrids=1 gridsAvailable=1 gridsKnown=1 stepCount=2 projStepCount=4 ballpark=0 vertBallpark=0 isNull=0
name=Inverse of NAD83 to WGS 84 (55) + UTM zone 17N area=0.236125 accuracy=2 isPROJExportable=1 hasGrids=1 gridsAvailable=1 gridsKnown=1 stepCount=2 projStepCount=4 ballpark=0 vertBallpark=0 isNull=0
name=Inverse of NAD83 to WGS 84 (40) + UTM zone 17N area=0.228083 accuracy=2 isPROJExportable=1 hasGrids=1 gridsAvailable=1 gridsKnown=1 stepCount=2 projStepCount=4 ballpark=0 vertBallpark=0 isNull=0
name=Inverse of NAD83 to WGS 84 (17) + UTM zone 17N area=0.22057 accuracy=2 isPROJExportable=1 hasGrids=1 gridsAvailable=1 gridsKnown=1 stepCount=2 projStepCount=4 ballpark=0 vertBallpark=0 isNull=0
name=Inverse of NAD83 to WGS 84 (22) + UTM zone 17N area=0.0751742 accuracy=2 isPROJExportable=1 hasGrids=1 gridsAvailable=1 gridsKnown=1 stepCount=2 projStepCount=4 ballpark=0 vertBallpark=0 isNull=0
name=Inverse of NAD83 to WGS 84 (51) + UTM zone 17N area=0.217178 accuracy=2 isPROJExportable=1 hasGrids=1 gridsAvailable=1 gridsKnown=1 stepCount=2 projStepCount=4 ballpark=0 vertBallpark=0 isNull=0
name=Inverse of NAD83 to WGS 84 (49) + UTM zone 17N area=0.182359 accuracy=2 isPROJExportable=1 hasGrids=1 gridsAvailable=1 gridsKnown=1 stepCount=2 projStepCount=4 ballpark=0 vertBallpark=0 isNull=0
name=Inverse of NAD83 to WGS 84 (56) + UTM zone 17N area=0.222232 accuracy=2 isPROJExportable=1 hasGrids=1 gridsAvailable=1 gridsKnown=1 stepCount=2 projStepCount=4 ballpark=0 vertBallpark=0 isNull=0
name=Inverse of NAD83 to WGS 84 (57) + UTM zone 17N area=0.0937974 accuracy=2 isPROJExportable=1 hasGrids=1 gridsAvailable=1 gridsKnown=1 stepCount=2 projStepCount=4 ballpark=0 vertBallpark=0 isNull=0
name=Inverse of NAD83 to WGS 84 (37) + UTM zone 17N area=0.0565287 accuracy=2 isPROJExportable=1 hasGrids=1 gridsAvailable=1 gridsKnown=1 stepCount=2 projStepCount=4 ballpark=0 vertBallpark=0 isNull=0
name=Inverse of NAD83 to WGS 84 (24) + UTM zone 17N area=0.0380414 accuracy=2 isPROJExportable=1 hasGrids=1 gridsAvailable=1 gridsKnown=1 stepCount=2 projStepCount=4 ballpark=0 vertBallpark=0 isNull=0
name=Ballpark geographic offset from WGS 84 to NAD83 + UTM zone 17N area=3.5449 accuracy=-1 isPROJExportable=1 hasGrids=0 gridsAvailable=1 gridsKnown=1 stepCount=2 projStepCount=3 ballpark=1 vertBallpark=0 isNull=0

Answer 2

The EPSG database has an attribute "accuracy" for the transformations. For example this one https://epsg.org/transformation_1107/Afgooye-to-WGS-84-1.html? has an accuracy of 44 m.

Proj stores the data from the EPSG database into an SQLite database "proj.db" that is located in the Proj data directory. I think that the transformations are stored into 3 tables: helmert_transformation_table, grid_transformation, and other_transformation. All those tables have the accuracy column. However, it seems that the value of the accuracy can be 0.0 or NULL. This is a screenshot about the helmert transformations table structure.