1

I have a database with many coordinates that were transformed using PyProj 1.9.6 (with PROJ 4.9.3 in wheels). I want to ensure that the positioning remains consistent after ugprading to PyProj 3.0.1 (with PROJ 7.2.1 in wheels), and beyond. I'm pretty sure this requires picking the correct pyproj.Transformer within the pyproj.transformer.TransformerGroup. The default Transformer does not always replicate PyProj 1.9.6, presumably because of changes in PROJ's algorithm for picking transformations, availability of datum shift grids, etc. Is there a way to determine which transformation PyProj 1.9.6 would have picked from within PyProj 3, for consistency purposes?

Also, is Transformer.description the best way to canonically describe the transformation that was chosen? Doing Transformer.to_json_dict() and extracting the transformation step out of the pipeline sometimes yields an EPSG code, but not always. I just want something that I can save in a database that will continue to work, even with future versions of Py(Proj).


For example, let's unproject some EPSG 27700 x,y coordinates to lng,lat in PyProj 1.9.6.

>>> import pyproj
>>> pyproj.proj_version_str
'4.9.3'
>>> x, y = 533595.48, 180463.35
>>> lng, lat = pyproj.transform(pyproj.Proj(init='epsg:27700'), pyproj.Proj(init='epsg:4326'), x, y)
>>> lng, lat
(-0.07640596925708688, 51.50731699924586)

In PyProj 3.0.1, the default transformation Inverse of British National Grid + OSGB 1936 to WGS 84 (9) produces a point with a geodesic error of 1.81 meters. The second transformation Inverse of British National Grid + OSGB 1936 to WGS 84 (6) is the only with negligible error, so presumably that's what PyProj 1.9.6 had used. I could repeat this process for each of the 700+ different coordinate systems in this database, but I'm wondering if there's an easier and more precise way.

>>> import pyproj
>>> pyproj.proj_version_str
'7.2.1'
>>> pyproj.network.set_network_enabled(True)
>>> geod = pyproj.geod.Geod(ellps='WGS84')
>>> x, y = 533595.48, 180463.35
>>> lng, lat = -0.07640596925708688, 51.50731699924586 # from PyProj 1.x
>>> tg = pyproj.transformer.TransformerGroup(27700, 4326)
>>> for i, t in enumerate(tg.transformers):
...     lat2, lng2 = t.transform(x, y)
...     dist = geod.inv(lng, lat, lng2, lat2)[-1]
...     epsg = next(s for s in t.to_json_dict()['steps'] if s['type'] == 'Transformation').get('id', {}).get('code', '???')
...     print(f'Transformer {i} is {dist:.2f} meters off: EPSG {epsg}, {t.description}') 
... 
Transformer 0 is 1.81 meters off: EPSG 7710, Inverse of British National Grid + OSGB 1936 to WGS 84 (9)
Transformer 1 is 0.00 meters off: EPSG 1314, Inverse of British National Grid + OSGB 1936 to WGS 84 (6)
Transformer 2 is 6.16 meters off: EPSG 1195, Inverse of British National Grid + OSGB 1936 to WGS 84 (1)
Transformer 3 is 1.08 meters off: EPSG 1196, Inverse of British National Grid + OSGB 1936 to WGS 84 (2)
Transformer 4 is 1.74 meters off: EPSG 1197, Inverse of British National Grid + OSGB 1936 to WGS 84 (3)
Transformer 5 is 16.83 meters off: EPSG 1198, Inverse of British National Grid + OSGB 1936 to WGS 84 (4)
Transformer 6 is 4.43 meters off: EPSG 1199, Inverse of British National Grid + OSGB 1936 to WGS 84 (5)
Transformer 7 is 3.49 meters off: EPSG 5622, Inverse of British National Grid + OSGB 1936 to WGS 84 (8)
Transformer 8 is 125.57 meters off: EPSG ???, Inverse of British National Grid + Ballpark geographic offset from OSGB 1936 to WGS 84

1 Answer 1

1

Doing what you are doing with the transformer group would be what I would recommend.

is Transformer.description the best way to canonically describe the transformation that was chosen?

The PROJ pipeline string from the transformer group is your best bet at consistency. You can get it using Transformer.to_proj4() or Transformer.definition depending on your version of pyproj.

4
  • Thanks for clarifying, @snowman2. I may set up a microservice that has PyProj 1.9.6 in a separate virtual environment from PyProj 3 so I can verify that it's always picking the same transformation from the transformer group. Or at least a transformation with a negligible difference.
    – dmahr
    Jul 7, 2021 at 12:47
  • 1
    You can use Transformer.from_pipeline once you get the pipeline string.
    – snowman2
    Jul 7, 2021 at 12:52
  • Follow-up question: given two Transformer objects, is there any way to know if they are the inverse of one another? I feel like the nth transformer in the forward group should match the nth transformer in the inverse group, but I'm not positive.
    – dmahr
    Jul 20, 2021 at 10:51
  • 1
    There might be, but I am not sure how you would tell in a reliable way.
    – snowman2
    Jul 20, 2021 at 11:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.