Transformation failed after PromoteTo3D

[aharondavid]

PROJ version: 9.3
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"]]]"

OGRSpatialReference srcOGRSpatialReference;
srcOGRSpatialReference.importFromWkt(srcWKT.c_str());
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]

The issue is that your source WKT is invalid. It lacks a required UNIT[] node in the VERT_CS[] node. You should check the return code of importFromWkt() against OGRERR_NONE

In #3926, I've submitted an enhancement to tolerate that in non-strict mode of WKT (which is used by GDAL), with a warning, and default to metre, as was done by the GDAL 2.4 WKT1 parser.

[aharondavid]

The original WKT is:
COMPD_CS["ENU (-77.410692720411:39.4145340892321) + EGM96 geoid height",PROJCS["ENU (-77.410692720411:39.4145340892321)",GEOGCS["GCS_WGS_1984",DATUM["D_unknown",SPHEROID["WGS84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Orthographic"],PARAMETER["Latitude_Of_Center",39.4145340892321],PARAMETER["Longitude_Of_Center",-77.410692720411],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]],VERT_CS["EGM96 geoid height",VERT_DATUM["EGM96 geoid",2005,AUTHORITY["EPSG","5171"],EXTENSION["PROJ4_GRIDS","egm96_15.gtx"]],AXIS["Up",UP],AUTHORITY["EPSG","5773"]]]

OGRSpatialReference.importFromWkt is failed on that WKT.

So, I did replace it by corrected WKT which OGRSpatialReference.importFromWkt is succeeded :

COMPD_CS["ENU (-77.410692720411:39.4145340892321) + EGM96 geoid height",PROJCS["ENU (-77.410692720411:39.4145340892321)",GEOGCS["GCS_WGS_1984",DATUM["D_unknown",SPHEROID["WGS84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Orthographic"],PARAMETER["Latitude_Of_Center",39.4145340892321],PARAMETER["Longitude_Of_Center",-77.410692720411],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]],VERT_CS["EGM96 geoid height",VERT_DATUM["EGM96 geoid",2005,AUTHORITY["EPSG","5171"],EXTENSION["PROJ4_GRIDS","egm96_15.gtx"]],AXIS["Up",UP],AUTHORITY["EPSG","5773"],UNIT["m",1.0]]]

Then I did OGRSpatialReference.exportToWkt and got that different WKT:

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"]]]

can you explain that behavior ?

[rouault]

Then I did OGRSpatialReference.exportToWkt and got that different WKT:

It is equivalent.

can you explain that behavior ?

PROJ uses an object model closer to WKT2 than WKT1. When you import from WKT1, it normalises to the object representation, and re-export from it.

[aharondavid]

but if the destination WKT is COMPD_CS the transformation is work without any error:

COMPD_CS["WGS84 Coordinate System + egm96 geoid height",
    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"]],
    VERT_CS["EGM96 height",
        VERT_DATUM["EGM96 geoid",2005,
            AUTHORITY["EPSG","5171"]],
        UNIT["metre",1,
            AUTHORITY["EPSG","9001"]],
        AXIS["Gravity-related height",UP],
        AUTHORITY["EPSG","5773"]]]

[rouault]

but if the destination WKT is COMPD_CS the transformation is work without any error:

yes, note that one has a UNIT[] node inside VERT_CS[].

[aharondavid]

But you said that the problem is in the source WKT, not in destination WKT

[rouault]

But you said that the problem is in the source WKT, not in destination WKT

Eh, the problem in your original reproducer was in your source WKT indeed. And then in #3925 (comment) you're going with another test case that doesn't exhibit the characteristics of your original report. I'm stopping this conversation here

[aharondavid]

Ok, Thanks