Bug in Pyproj Transformer class when using a BoundCRS object (incorrect results) #4517
Description
The issue was originally created in the Pyproj issue board, but it was recommended to move it to PROJ. Here is the link to the issue in Pyproj:
pyproj4/pyproj#1467
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)Comment
I understand that datum transformations are typically defined between geographic CRSs, specifically EPSG:4230 and EPSG:4326 in this case. However, it should be possible to create a BoundCRS object as described above, with Pyproj/PROJ referencing the base CRS of the projected CRS. The TransformerGroup class constructs all necessary steps, including both conversion and datum transformation. A similar functionality in the BoundCRS class would be expected and preferable.
Pyproj should at least throw an error in the scenario above rather than selecting a different transformation than the one provided when a projected CRS is used as input. If a user specifies a particular transformation and Pyproj selects a different one without any warning or error message, this is undesirable behavior. This undermines the geodetic integrity of the library and should be considered a bug.
From a user perspective, it would be beneficial if the BoundCRS class had functionality for building a complete pipeline, including all necessary steps to transition from the source CRS to the target CRS, with the datum transformation explicitly defined by the user. This approach would simplify the process for users, allowing transformations to be executed in a single step rather than requiring three separate operations:
- Convert from projected to geographic (EPSG:23031 → EPSG:4230)
- Perform the datum transformation ED50 → WGS84 using the EPSG:1612 transformation
- Convert from geographic to projected (EPSG:4326 → EPSG:32631)
Furthermore, the system 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.
The BoundCRS class is powerful and allows users to easily tie CRSs and transformations together. It would be a significant improvement if it could exhibit similar behavior to TransformerGroup by breaking down operations into multiple steps and handling both conversions and transformations seamlessly.
Expected Output
wgs_84_proj_true_values = (572653.567, 7234083.246) # easting, northingEnvironment Information
pyproj 3.6.1
python 3.11.10
Installation method
pip wheel