[PROJ] Correct syntax of vertical datum

[Even Rouault]

Martin,

> Is this WKT invalid or just using some obscure feature that PROJ does
> not implement? 

Unless there's some hidden WKT spec I'm not aware of, this WKT is invalid. You 
cannot put a GEOGCS as a child node of a VERT_CS. And as this is a sort of 
WKT1, even the GEOGCS to express EPSG:4937 (ETRS89 as a Geographic 3D CRS) is 
invalid, as WKT1 spec restricts to 2 axis.

> And if the WKT is wrong, how can I produce a VERT_CS
> clause for EPSG:4937 that would be valid, so that I can notify the
> authority to fix the metadata?

In theory, you cannot express the concept of a Projected CRS + an ellipsoidal 
height as a compoundCRS. Because an ellipsoidal height cannot define a 
VerticalCRS (you must associate the horizontal coordinates as well), contrary 
to gravity-based orthometric heights which can be standalone.

That said... The old WKT1 specs allow a vertical datum to have a 2002 code 
expressing an ellipsoidal height. And US Lidar users have use that to create 
such WKT as you can see in https://github.com/OSGeo/PROJ/issues/2228

So if you transpose that example to yours, that would be:

COMPD_CS["ETRS89 / UTM zone 34N (N-E) + Ellipsoid (Meters)",
    PROJCS["ETRS89 / UTM zone 34N (N-E)",
        GEOGCS["ETRS89",
            DATUM["European_Terrestrial_Reference_System_1989",
                SPHEROID["GRS 1980",6378137,298.257222101,
                    AUTHORITY["EPSG","7019"]],
                AUTHORITY["EPSG","6258"]],
            PRIMEM["Greenwich",0,
                AUTHORITY["EPSG","8901"]],
            UNIT["degree",0.0174532925199433,
                AUTHORITY["EPSG","9122"]],
            AUTHORITY["EPSG","4258"]],
        PROJECTION["Transverse_Mercator"],
        PARAMETER["latitude_of_origin",0],
        PARAMETER["central_meridian",21],
        PARAMETER["scale_factor",0.9996],
        PARAMETER["false_easting",500000],
        PARAMETER["false_northing",0],
        UNIT["metre",1,
            AUTHORITY["EPSG","9001"]],
        AUTHORITY["EPSG","3046"]],
    VERT_CS["Ellipsoid (Meters)",
        VERT_DATUM["Ellipsoid",2002],
        UNIT["metre",1.0,
            AUTHORITY["EPSG","9001"]],
        AXIS["Up",UP]]]

Note that this will only be recognized properly by PROJ >= 7.1.1. Earlier 
versions would ingest that, but will not understand that it is an ellipsoidal 
VertCRS, and doing coordinate transformations with them will lead to improper 
transformation of the vertical component

If you provide this WKT as input of projinfo ( >= 7.1.1 ), you'll note that 
the output as WKT2:2019 is

PROJCRS["ETRS89 / UTM zone 34N (N-E)",
    BASEGEOGCRS["ETRS89",
        DATUM["European Terrestrial Reference System 1989",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4937]],
    CONVERSION["UTM zone 34N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",21,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]],
        ID["EPSG",16034]],
    CS[Cartesian,3],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["ellipsoidal height (h)",up,
            ORDER[3],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
    REMARK["Promoted to 3D from EPSG:3046"]]

So upon ingestion, this is no longer a CompoundCRS, but a PROJCRS whose 
BASEGEOGCRS is EPSG:4937, and with a cartesian coordinate system with 3 axis. 
(The REMARK is something new in PROJ master)
This construct is allowed per WKT:2019
( http://docs.opengeospatial.org/is/18-010r7/18-010r7.html , look for "For a 
projected 3D CRS") , but this remains a quite exotic object.

You can also generate the above WKT with:
projinfo EPSG:3046 --3d

Even

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com