Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Transformation failed after PromoteTo3D despite of WKT is correct and OGRSpatialReference .importFromWkt return OGRERR_NONE #3927

Closed
aharondavid opened this issue Oct 19, 2023 · 1 comment · Fixed by #3929
Assignees
Labels

Comments

@aharondavid
Copy link

aharondavid commented Oct 19, 2023

PROJ version: 9.3
GDAL: 3.1.2
OS: WIndows 11

Code to reproduced:


std:: string srcWKT = "COMPD_CS["ENU (-77.410692720411:39.4145340892321) + EGM96 geoid height",PROJCS["ENU (-77.410692720411:39.4145340892321)",GEOGCS["WGS 84",DATUM["unknown",SPHEROID["WGS84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.0174532925199433]],PROJECTION["Orthographic"],PARAMETER["latitude_of_origin",39.4145340892321],PARAMETER["central_meridian",-77.410692720411],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]],VERT_CS["EGM96 geoid height",VERT_DATUM["EGM96 geoid",2005,EXTENSION["PROJ4_GRIDS","egm96_15.gtx"],AUTHORITY["EPSG","5171"]],UNIT["m",1],AXIS["Up",UP],AUTHORITY["EPSG","5773"]]]"

// WKT is correct and include UNIT["m",1] in VERT_CS in the right place
 
OGRSpatialReference srcOGRSpatialReference;
int error = srcOGRSpatialReference.importFromWkt(srcWKT.c_str());
// error is OGRERR_NONE
srcOGRSpatialReference.PromoteTo3D(NULL);




std:: string destWKT = "GEOGCS["WGS84 Coordinate System",DATUM["WGS_1984",SPHEROID["WGS 1984",6378137,298.257223563],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]"

OGRSpatialReference destOGRSpatialReference;
destOGRSpatialReference.importFromWkt(destWKT.c_str());
destOGRSpatialReference.PromoteTo3D(NULL);


OGRCoordinateTransformation *pResult = OGRCreateCoordinateTransformation(&srcOGRSpatialReference,&destOGRSpatialReference);

//pResult is null

Error Result:
Cannot find coordinate operations from COMPOUNDCRS["ENU (-77.410692720411:39.4145340892321) + EGM96 geoid height",PROJCRS["ENU (-77.410692720411:39.4145340892321)",BASEGEOGCRS["WGS 84",DATUM["unknown",ELLIPSOID["WGS84",6378137,298.257223563,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["Degree",0.0174532925199433]]],CONVERSION["unnamed",METHOD["Orthographic (Spherical)"],PARAMETER["Latitude of natural origin",39.4145340892321,ANGLEUNIT["Degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",-77.410692720411,ANGLEUNIT["Degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["False easting",0,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]],BOUNDCRS[SOURCECRS[VERTCRS["EGM96 geoid height",VDATUM["EGM96 geoid"],CS[vertical,1],AXIS["up",up,LENGTHUNIT["m",1]],ID["EPSG",5773]]],TARGETCRS[GEOGCRS["WGS 84",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,3],AXIS["latitude",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["longitude",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],AXIS["ellipsoidal height",up,ORDER[3],LENGTHUNIT["metre",1]],ID["EPSG",4979]]],ABRIDGEDTRANSFORMATION["EGM96 geoid height to WGS 84 ellipsoidal height",METHOD["GravityRelatedHeight to Geographic3D"],PARAMETERFILE["Geoid (height correction) model file","egm96_15.gtx",ID["EPSG",8666]]]]]' to BOUNDCRS[SOURCECRS[GEOGCRS["WGS84 Coordinate System",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 1984",6378137,298.257223563,LENGTHUNIT["metre",1]],ID["EPSG",6326]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,3],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],AXIS["ellipsoidal height (h)",up,ORDER[3],LENGTHUNIT["metre",1,ID["EPSG",9001]]],REMARK["Promoted to 3D from EPSG:4326"]]],TARGETCRS[GEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 (Transit)"],MEMBER["World Geodetic System 1984 (G730)"],MEMBER["World Geodetic System 1984 (G873)"],MEMBER["World Geodetic System 1984 (G1150)"],MEMBER["World Geodetic System 1984 (G1674)"],MEMBER["World Geodetic System 1984 (G1762)"],MEMBER["World Geodetic System 1984 (G2139)"],ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,3],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],AXIS["ellipsoidal height (h)",up,ORDER[3],LENGTHUNIT["metre",1]],USAGE[SCOPE["Geodesy. Navigation and positioning using GPS satellite system."],AREA["World."],BBOX[-90,-180,90,180]],ID["EPSG",4979]]],ABRIDGEDTRANSFORMATION["Transformation from WGS84 Coordinate System to WGS84",METHOD["Position Vector transformation (geog2D domain)",ID["EPSG",9606]],PARAMETER["X-axis translation",0,ID["EPSG",8605]],PARAMETER["Y-axis translation",0,ID["EPSG",8606]],PARAMETER["Z-axis translation",0,ID["EPSG",8607]],PARAMETER["X-axis rotation",0,ID["EPSG",8608]],PARAMETER["Y-axis rotation",0,ID["EPSG",8609]],PARAMETER["Z-axis rotation",0,ID["EPSG",8610]],PARAMETER["Scale difference",1,ID["EPSG",8611]]]]'

@rouault
Copy link
Member

rouault commented Oct 19, 2023

Works on master when using projinfo directly

$ bin/projinfo -s 'COMPD_CS["ENU (-77.410692720411:39.4145340892321) + EGM96 geoid height",PROJCS["ENU (-77.410692720411:39.4145340892321)",GEOGCS["WGS 84",DATUM["unknown",SPHEROID["WGS84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.0174532925199433]],PROJECTION["Orthographic"],PARAMETER["latitude_of_origin",39.4145340892321],PARAMETER["central_meridian",-77.410692720411],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]],VERT_CS["EGM96 geoid height",VERT_DATUM["EGM96 geoid",2005,EXTENSION["PROJ4_GRIDS","egm96_15.gtx"],AUTHORITY["EPSG","5171"]],UNIT["m",1],AXIS["Up",UP],AUTHORITY["EPSG","5773"]]]' -t 'GEOGCS["WGS84 Coordinate System",DATUM["WGS_1984",SPHEROID["WGS 1984",6378137,298.257223563],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]' --3d
Candidate operations found: 1
-------------------------------------
Operation No. 1:

unknown id, Inverse of unnamed + Ballpark geographic offset from WGS 84 to WGS 84 + EGM96 geoid height to WGS 84 ellipsoidal height, unknown accuracy, World, has ballpark transformation, at least one grid missing

PROJ string:
+proj=pipeline
  +step +inv +proj=ortho +lat_0=39.4145340892321 +lon_0=-77.410692720411 +x_0=0
        +y_0=0 +ellps=WGS84
  +step +proj=vgridshift +grids=egm96_15.gtx +multiplier=1
  +step +proj=unitconvert +xy_in=rad +xy_out=deg
  +step +proj=axisswap +order=2,1

But not when using the WKT2 strings used internally by GDAL CoordinateTransformation:

$ bin/projinfo -s 'COMPOUNDCRS["ENU (-77.410692720411:39.4145340892321) + EGM96 geoid height",PROJCRS["ENU (-77.410692720411:39.4145340892321)",BASEGEOGCRS["WGS 84",DATUM["unknown",ELLIPSOID["WGS84",6378137,298.257223563,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["Degree",0.0174532925199433]]],CONVERSION["unnamed",METHOD["Orthographic",ID["EPSG",9840]],PARAMETER["Latitude of natural origin",39.4145340892321,ANGLEUNIT["Degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",-77.410692720411,ANGLEUNIT["Degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["False easting",0,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["easting",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["northing",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]],BOUNDCRS[SOURCECRS[VERTCRS["EGM96 geoid height",VDATUM["EGM96 geoid"],CS[vertical,1],AXIS["up",up,LENGTHUNIT["m",1]],ID["EPSG",5773]]],TARGETCRS[GEOGCRS["WGS 84",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,3],AXIS["latitude",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["longitude",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],AXIS["ellipsoidal height",up,ORDER[3],LENGTHUNIT["metre",1]],ID["EPSG",4979]]],ABRIDGEDTRANSFORMATION["EGM96 geoid height to WGS 84 ellipsoidal height",METHOD["GravityRelatedHeight to Geographic3D"],PARAMETERFILE["Geoid (height correction) model file","egm96_15.gtx",ID["EPSG",8666]]]]]' -t 'BOUNDCRS[SOURCECRS[GEOGCRS["WGS84 Coordinate System",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 1984",6378137,298.257223563,LENGTHUNIT["metre",1]],ID["EPSG",6326]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,3],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],AXIS["ellipsoidal height (h)",up,ORDER[3],LENGTHUNIT["metre",1,ID["EPSG",9001]]],REMARK["Promoted to 3D from EPSG:4326"]]],TARGETCRS[GEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 (Transit)"],MEMBER["World Geodetic System 1984 (G730)"],MEMBER["World Geodetic System 1984 (G873)"],MEMBER["World Geodetic System 1984 (G1150)"],MEMBER["World Geodetic System 1984 (G1674)"],MEMBER["World Geodetic System 1984 (G1762)"],MEMBER["World Geodetic System 1984 (G2139)"],ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,3],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],AXIS["ellipsoidal height (h)",up,ORDER[3],LENGTHUNIT["metre",1]],USAGE[SCOPE["Geodesy. Navigation and positioning using GPS satellite system."],AREA["World."],BBOX[-90,-180,90,180]],ID["EPSG",4979]]],ABRIDGEDTRANSFORMATION["Transformation from WGS84 Coordinate System to WGS84",METHOD["Position Vector transformation (geog2D domain)",ID["EPSG",9606]],PARAMETER["X-axis translation",0,ID["EPSG",8605]],PARAMETER["Y-axis translation",0,ID["EPSG",8606]],PARAMETER["Z-axis translation",0,ID["EPSG",8607]],PARAMETER["X-axis rotation",0,ID["EPSG",8608]],PARAMETER["Y-axis rotation",0,ID["EPSG",8609]],PARAMETER["Z-axis rotation",0,ID["EPSG",8610]],PARAMETER["Scale difference",1,ID["EPSG",8611]]]]'
Candidate operations found: 0

Can be simplified to:

$  bin/projinfo -s 'COMPOUNDCRS["ENU (-77.410692720411:39.4145340892321) + EGM96 geoid height",PROJCRS["ENU (-77.410692720411:39.4145340892321)",BASEGEOGCRS["WGS 84",DATUM["unknown",ELLIPSOID["WGS84",6378137,298.257223563,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["Degree",0.0174532925199433]]],CONVERSION["unnamed",METHOD["Orthographic",ID["EPSG",9840]],PARAMETER["Latitude of natural origin",39.4145340892321,ANGLEUNIT["Degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",-77.410692720411,ANGLEUNIT["Degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["False easting",0,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["easting",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["northing",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]],BOUNDCRS[SOURCECRS[VERTCRS["EGM96 geoid height",VDATUM["EGM96 geoid"],CS[vertical,1],AXIS["up",up,LENGTHUNIT["m",1]],ID["EPSG",5773]]],TARGETCRS[GEOGCRS["WGS 84",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,3],AXIS["latitude",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["longitude",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],AXIS["ellipsoidal height",up,ORDER[3],LENGTHUNIT["metre",1]],ID["EPSG",4979]]],ABRIDGEDTRANSFORMATION["EGM96 geoid height to WGS 84 ellipsoidal height",METHOD["GravityRelatedHeight to Geographic3D"],PARAMETERFILE["Geoid (height correction) model file","egm96_15.gtx",ID["EPSG",8666]]]]]' -t EPSG:4979
Candidate operations found: 0

But using a GeogCRS as the horizontal CRS of the CompoundCRS works:

$ bin/projinfo -s 'COMPOUNDCRS["test",GEOGCRS["WGS 84",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["latitude",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["longitude",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]]],BOUNDCRS[SOURCECRS[VERTCRS["EGM96 geoid height",VDATUM["EGM96 geoid"],CS[vertical,1],AXIS["up",up,LENGTHUNIT["m",1]],ID["EPSG",5773]]],TARGETCRS[GEOGCRS["WGS 84",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,3],AXIS["latitude",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["longitude",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],AXIS["ellipsoidal height",up,ORDER[3],LENGTHUNIT["metre",1]],ID["EPSG",4979]]],ABRIDGEDTRANSFORMATION["EGM96 geoid height to WGS 84 ellipsoidal height",METHOD["GravityRelatedHeight to Geographic3D"],PARAMETERFILE["Geoid (height correction) model file","egm96_15.gtx",ID["EPSG",8666]]]]]' -t EPSG:4979
Candidate operations found: 1
-------------------------------------
Operation No. 1:

unknown id, EGM96 geoid height to WGS 84 ellipsoidal height, unknown accuracy, unknown domain of validity, at least one grid missing

PROJ string:
+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +xy_in=deg +xy_out=rad
  +step +proj=vgridshift +grids=egm96_15.gtx +multiplier=1
  +step +proj=unitconvert +xy_in=rad +xy_out=deg
  +step +proj=axisswap +order=2,1

@rouault rouault self-assigned this Oct 19, 2023
rouault added a commit to rouault/PROJ that referenced this issue Oct 19, 2023
…RS of unknown horizontal datum + boundCRS of vertical (fixes OSGeo#3927)
rouault added a commit that referenced this issue Oct 27, 2023
…RS of unknown horizontal datum + boundCRS of vertical (fixes #3927)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
2 participants