Bug in Transformer class when using a BoundCRS object (incorrect results)

[phaarnes]

Example Code showing the issue where the transformation part of the BoundCRS object is ignored

from pyproj import CRS,Transformer
from pyproj.crs import BoundCRS
import numpy as np
np.set_printoptions(suppress=True)


wgs_84_proj_true_values = (572653.567, 7234083.246) # easting, northing


# Transform from ED 50 Projected to WGS 84 Projected (Pyproj pick the EPSG:1139 transformation for some reason)
src_crs_wkt2 ='PROJCRS["ED50 / UTM zone 31N",BASEGEOGCRS["ED50",DATUM["European Datum 1950",ELLIPSOID["International 1924",6378388,297,LENGTHUNIT["metre",1,ID["EPSG",9001]],ID["EPSG",7022]],ID["EPSG",6230]],ID["EPSG",4230]],CONVERSION["UTM zone 31N",METHOD["Transverse Mercator",ID["EPSG",9807]],PARAMETER["Latitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9102]],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",3,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9102]],ID["EPSG",8802]],PARAMETER["Scale factor at natural origin",0.9996,SCALEUNIT["unity",1,ID["EPSG",9201]],ID["EPSG",8805]],PARAMETER["False easting",500000,LENGTHUNIT["metre",1,ID["EPSG",9001]],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1,ID["EPSG",9001]],ID["EPSG",8807]],ID["EPSG",16031]],CS[Cartesian,2,ID["EPSG",4400]],AXIS["Easting (E)",east],AXIS["Northing (N)",north],LENGTHUNIT["metre",1,ID["EPSG",9001]],ID["EPSG",23031]]'
trg_crs_wkt2 = 'PROJCRS["WGS 84 / UTM zone 31N",BASEGEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble", MEMBER["World Geodetic System 1984 (Transit)", ID["EPSG",1166]], MEMBER["World Geodetic System 1984 (G730)", ID["EPSG",1152]], MEMBER["World Geodetic System 1984 (G873)", ID["EPSG",1153]], MEMBER["World Geodetic System 1984 (G1150)", ID["EPSG",1154]], MEMBER["World Geodetic System 1984 (G1674)", ID["EPSG",1155]], MEMBER["World Geodetic System 1984 (G1762)", ID["EPSG",1156]], MEMBER["World Geodetic System 1984 (G2139)", ID["EPSG",1309]], ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1,ID["EPSG",9001]],ID["EPSG",7030]], ENSEMBLEACCURACY[2],ID["EPSG",6326]],ID["EPSG",4326]],CONVERSION["UTM zone 31N",METHOD["Transverse Mercator",ID["EPSG",9807]],PARAMETER["Latitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9102]],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",3,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9102]],ID["EPSG",8802]],PARAMETER["Scale factor at natural origin",0.9996,SCALEUNIT["unity",1,ID["EPSG",9201]],ID["EPSG",8805]],PARAMETER["False easting",500000,LENGTHUNIT["metre",1,ID["EPSG",9001]],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1,ID["EPSG",9001]],ID["EPSG",8807]],ID["EPSG",16031]],CS[Cartesian,2,ID["EPSG",4400]],AXIS["Easting (E)",east],AXIS["Northing (N)",north],LENGTHUNIT["metre",1,ID["EPSG",9001]],ID["EPSG",32631]]'
transformation_wkt2 = 'COORDINATEOPERATION["ED50 to WGS 84 (23)",VERSION["EPSG-Nor N62 2001"],SOURCECRS[GEOGCRS["ED50",DATUM["European Datum 1950",ELLIPSOID["International 1924",6378388,297,LENGTHUNIT["metre",1,ID["EPSG",9001]],ID["EPSG",7022]],ID["EPSG",6230]],CS[ellipsoidal,2,ID["EPSG",6422]],AXIS["Geodetic latitude (Lat)",north],AXIS["Geodetic longitude (Lon)",east],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9102]],ID["EPSG",4230]]],TARGETCRS[GEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble", MEMBER["World Geodetic System 1984 (Transit)", ID["EPSG",1166]], MEMBER["World Geodetic System 1984 (G730)", ID["EPSG",1152]], MEMBER["World Geodetic System 1984 (G873)", ID["EPSG",1153]], MEMBER["World Geodetic System 1984 (G1150)", ID["EPSG",1154]], MEMBER["World Geodetic System 1984 (G1674)", ID["EPSG",1155]], MEMBER["World Geodetic System 1984 (G1762)", ID["EPSG",1156]], MEMBER["World Geodetic System 1984 (G2139)", ID["EPSG",1309]], ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1,ID["EPSG",9001]],ID["EPSG",7030]], ENSEMBLEACCURACY[2],ID["EPSG",6326]],CS[ellipsoidal,2,ID["EPSG",6422]],AXIS["Geodetic latitude (Lat)",north],AXIS["Geodetic longitude (Lon)",east],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9102]],ID["EPSG",4326]]],METHOD["Position Vector transformation (geog2D domain)",ID["EPSG",9606]],PARAMETER["X-axis translation",-116.641,LENGTHUNIT["metre",1,ID["EPSG",9001]],ID["EPSG",8605]],PARAMETER["Y-axis translation",-56.931,LENGTHUNIT["metre",1,ID["EPSG",9001]],ID["EPSG",8606]],PARAMETER["Z-axis translation",-110.559,LENGTHUNIT["metre",1,ID["EPSG",9001]],ID["EPSG",8607]],PARAMETER["X-axis rotation",0.893,ANGLEUNIT["arc-second",0.0000048481368111,ID["EPSG",9104]],ID["EPSG",8608]],PARAMETER["Y-axis rotation",0.921,ANGLEUNIT["arc-second",0.0000048481368111,ID["EPSG",9104]],ID["EPSG",8609]],PARAMETER["Z-axis rotation",-0.917,ANGLEUNIT["arc-second",0.0000048481368111,ID["EPSG",9104]],ID["EPSG",8610]],PARAMETER["Scale difference",-3.52,SCALEUNIT["parts per million",0.000001,ID["EPSG",9202]],ID["EPSG",8611]],OPERATIONACCURACY[1],ID["EPSG",1612]]'
ed50_proj = (572742.772, 7234299.211)


bound_crs = BoundCRS(
        source_crs = CRS.from_wkt(src_crs_wkt2),
        target_crs = CRS.from_wkt(trg_crs_wkt2),
        transformation = Transformer.from_pipeline(transformation_wkt2)
        )


# Create a transformer object for transforming coordinates
transformer = Transformer.from_crs(bound_crs, bound_crs.target_crs, always_xy=True)

# Transform the coordinates to WGS84
wgs84_proj_coords = transformer.transform(*ed50_proj)

diff_proj = np.array(wgs84_proj_coords) - np.array(wgs_84_proj_true_values)
print(diff_proj)


## Get the same as when doing
Transformer.from_crs("EPSG:23031", "EPSG:32631").transform(*ed50_proj)

Problem description

When the target CRS of a BoundCRS object is a Projected CRS, the Transformer class ignores the transformations provided in the BoundCRS. It just picks the "best available" transformation as when you run for example:

Transformer.from_crs("EPSG:23031", "EPSG:32631").transform(*ed50_proj)

However, when if the target CRS is a geographic CRS, the transformation part of the BoundCRS object is used. Here is a example code showing correct results if target CRS is geographic:

from pyproj import CRS,Transformer
from pyproj.crs import BoundCRS
import numpy as np
np.set_printoptions(suppress=True)


wgs_84_geog_true_values = (4.5537099576907, 65.22193148979875) # lon, lat

#%% Transform from ED 50 Projected to WGS 84 Geographic
src_crs_wkt2 ='PROJCRS["ED50 / UTM zone 31N",BASEGEOGCRS["ED50",DATUM["European Datum 1950",ELLIPSOID["International 1924",6378388,297,LENGTHUNIT["metre",1,ID["EPSG",9001]],ID["EPSG",7022]],ID["EPSG",6230]],ID["EPSG",4230]],CONVERSION["UTM zone 31N",METHOD["Transverse Mercator",ID["EPSG",9807]],PARAMETER["Latitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9102]],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",3,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9102]],ID["EPSG",8802]],PARAMETER["Scale factor at natural origin",0.9996,SCALEUNIT["unity",1,ID["EPSG",9201]],ID["EPSG",8805]],PARAMETER["False easting",500000,LENGTHUNIT["metre",1,ID["EPSG",9001]],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1,ID["EPSG",9001]],ID["EPSG",8807]],ID["EPSG",16031]],CS[Cartesian,2,ID["EPSG",4400]],AXIS["Easting (E)",east],AXIS["Northing (N)",north],LENGTHUNIT["metre",1,ID["EPSG",9001]],ID["EPSG",23031]]'
trg_crs_wkt2 = 'GEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble", MEMBER["World Geodetic System 1984 (Transit)", ID["EPSG",1166]], MEMBER["World Geodetic System 1984 (G730)", ID["EPSG",1152]], MEMBER["World Geodetic System 1984 (G873)", ID["EPSG",1153]], MEMBER["World Geodetic System 1984 (G1150)", ID["EPSG",1154]], MEMBER["World Geodetic System 1984 (G1674)", ID["EPSG",1155]], MEMBER["World Geodetic System 1984 (G1762)", ID["EPSG",1156]], MEMBER["World Geodetic System 1984 (G2139)", ID["EPSG",1309]], ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1,ID["EPSG",9001]],ID["EPSG",7030]], ENSEMBLEACCURACY[2],ID["EPSG",6326]],CS[ellipsoidal,2,ID["EPSG",6422]],AXIS["Geodetic latitude (Lat)",north],AXIS["Geodetic longitude (Lon)",east],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9102]],ID["EPSG",4326]]'
transformation_wkt2 = 'COORDINATEOPERATION["ED50 to WGS 84 (23)",VERSION["EPSG-Nor N62 2001"],SOURCECRS[GEOGCRS["ED50",DATUM["European Datum 1950",ELLIPSOID["International 1924",6378388,297,LENGTHUNIT["metre",1,ID["EPSG",9001]],ID["EPSG",7022]],ID["EPSG",6230]],CS[ellipsoidal,2,ID["EPSG",6422]],AXIS["Geodetic latitude (Lat)",north],AXIS["Geodetic longitude (Lon)",east],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9102]],ID["EPSG",4230]]],TARGETCRS[GEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble", MEMBER["World Geodetic System 1984 (Transit)", ID["EPSG",1166]], MEMBER["World Geodetic System 1984 (G730)", ID["EPSG",1152]], MEMBER["World Geodetic System 1984 (G873)", ID["EPSG",1153]], MEMBER["World Geodetic System 1984 (G1150)", ID["EPSG",1154]], MEMBER["World Geodetic System 1984 (G1674)", ID["EPSG",1155]], MEMBER["World Geodetic System 1984 (G1762)", ID["EPSG",1156]], MEMBER["World Geodetic System 1984 (G2139)", ID["EPSG",1309]], ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1,ID["EPSG",9001]],ID["EPSG",7030]], ENSEMBLEACCURACY[2],ID["EPSG",6326]],CS[ellipsoidal,2,ID["EPSG",6422]],AXIS["Geodetic latitude (Lat)",north],AXIS["Geodetic longitude (Lon)",east],ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9102]],ID["EPSG",4326]]],METHOD["Position Vector transformation (geog2D domain)",ID["EPSG",9606]],PARAMETER["X-axis translation",-116.641,LENGTHUNIT["metre",1,ID["EPSG",9001]],ID["EPSG",8605]],PARAMETER["Y-axis translation",-56.931,LENGTHUNIT["metre",1,ID["EPSG",9001]],ID["EPSG",8606]],PARAMETER["Z-axis translation",-110.559,LENGTHUNIT["metre",1,ID["EPSG",9001]],ID["EPSG",8607]],PARAMETER["X-axis rotation",0.893,ANGLEUNIT["arc-second",0.0000048481368111,ID["EPSG",9104]],ID["EPSG",8608]],PARAMETER["Y-axis rotation",0.921,ANGLEUNIT["arc-second",0.0000048481368111,ID["EPSG",9104]],ID["EPSG",8609]],PARAMETER["Z-axis rotation",-0.917,ANGLEUNIT["arc-second",0.0000048481368111,ID["EPSG",9104]],ID["EPSG",8610]],PARAMETER["Scale difference",-3.52,SCALEUNIT["parts per million",0.000001,ID["EPSG",9202]],ID["EPSG",8611]],OPERATIONACCURACY[1],ID["EPSG",1612]]'
ed50_proj = (572742.772, 7234299.211)

bound_crs = BoundCRS(
        source_crs = CRS.from_wkt(src_crs_wkt2),
        target_crs = CRS.from_wkt(trg_crs_wkt2),
        transformation = Transformer.from_pipeline(transformation_wkt2)
        )


# Create a transformer object for transforming coordinates
transformer = Transformer.from_crs(bound_crs, bound_crs.target_crs, always_xy=True)

# Transform the coordinates to WGS84
wgs84_geog_coords = transformer.transform(*ed50_proj)

diff_geog = np.array(wgs84_geog_coords) - np.array(wgs_84_geog_true_values)
print(diff_geog)

Expected Output

wgs_84_proj_true_values = (572653.567, 7234083.246) # easting, northing

Environment Information

pyproj 3.6.1
python 3.11.10

Installation method

pip wheel

[jjimenezshaw]

The target CRS must be geographic. In your case WGS 84 not WGS 84 / UTM zone 31N