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

[phaarnes]

@jjimenezshaw , @snowman2

Hi,

Sorry for my late response. I understand that the datum transformation is defined between geographic CRS, in this case, EPSG:4230 and EPSG:4326. However, I believe it should be possible to create a BoundCRS object as described above, and Pyproj/PROJ should then reference the base CRS of the projected CRS. The TransformerGroup class, for instance, builds all the necessary steps, including both conversion and datum transformation. A similar functionality in the BoundCRS class would not only be expected but also preferable.

Regardless, Pyproj should at least throw an error in the case above rather than selecting a different transformation than provided when a projected CRS is used as input. If a user defines a specific transformation but Pyproj selects a different one without any warning or error message, that is not desirable behavior. I find this behavior problematic, as it does not support the geodetic integrity of the library and should be labeled as a bug.

The BoundCRS class is a powerful and great class, allowing users to easily tie CRSs and transformations together. It would be a significant improvement if it could exhibit similar behavior to TransformerGroup by breaking it down down into multiple steps/operations and handling both conversions and transformations seamlessly.

[snowman2]

@phaarnes those are interesting suggestions.

Pyproj/PROJ should then reference the base CRS of the projected CRS.

I hesitate to do something like that. That could cause confusion for users learning about building CRS. Maybe throw an an error, as you suggested, that would tell users to pass in a geographic CRS.

[phaarnes]

Yes, I got your point on that one. Personally, I would prefer that the BoundCRS class had similar functionality for building a complete pipeline, including all the necessary steps to get from the source CRS to the target CRS, where the datum transformation is explicitly defined by the user. This would make things much easier for the user, as the transformation could be done in a single step from the user's perspective, rather than requiring three separate steps:

1. Convert from projected to geographic (EPSG:23031 → EPSG:4230)
2. Perform the datum transformation ED50 → WGS84 using the EPSG:1612 transformation
3. Convert from geographic to projected (EPSG:4326 → EPSG:32631)

Furthermore, I believe it should throw an error and not produce any results if there is an issue with the user input. It should not automatically select the best available transformation without informing the user. This makes me somewhat reluctant to use this class, as I feel I cannot fully trust that it is actually using the transformation I provided.