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