Support AUTH:nnnnn+ellipsoidal syntax in osgeo::proj::io::createFromUserInput

[hernando]

I would like to pitch an idea for extending the behavior of osgeo::proj::io::createFromUserInput that I would offer to implement if it gets accepted.
Quite often I find myself having to write things like cs2cs "$(projinfo -q -o WKT2:2019 --single-line --3d EPSG:32632)" "$(projinfo -q -o WKT2:2019 --single-line --3d EPSG:2056)" -d 12 in order to get the vertical axis transformed between CRSs with different ellipsoids. I'm aware that this doesn't always make sense (e.g. specially with OSGB36), but in the previous example it's perfectly valid as far as I know.
Since there's already an option for 3D promotion in projinfo, I believe it would be handy to make the same functionality available through the CRS identifiers accepted by createFromUserInput, so the example above would become: cs2cs EPSG:32632+ellipsoidal EPSG:2056+ellipsoidal.

The open question to me is how to handle extensions to 3D of projected CRSs in units other than meter. I know that the promoteTo3Dmember function always adds the z axis in meter regardless of the unit of the other axes. I'm fine with this behaviour, however some users may find it surprising when writing something like EPSG:2234+ellipsoidal, where EPSG:2234 is a projected CRS using ftUS.

[hernando]

Sorry for being pushy, but I would like to rely on this behavior for something else than just a handy shortcut in the command line, so if it's an absolute no go, I would like to know to find another solution to my problem.

[kbevers]

Seems to me that implementing the --3d option in cs2cs would solve your problem.

With projinfo we can get a pipeline with almost the same syntax as would be passed to cs2cs:

$ projinfo -s EPSG:32632 -t EPSG:2056 -o PROJ --3d
Candidate operations found: 1
-------------------------------------
Operation No. 1:

unknown id, Inverse of UTM zone 32N + Inverse of CH1903+ to WGS 84 (1) + Swiss Oblique Mercator 1995, 1 m, Liechtenstein; Switzerland.

PROJ string:
+proj=pipeline
  +step +inv +proj=utm +zone=32 +ellps=WGS84
  +step +proj=cart +ellps=WGS84
  +step +proj=helmert +x=-674.374 +y=-15.056 +z=-405.346
  +step +inv +proj=cart +ellps=bessel
  +step +proj=somerc +lat_0=46.9524055555556 +lon_0=7.43958333333333 +k_0=1
        +x_0=2600000 +y_0=1200000 +ellps=bessel

which is equivalent to your cs2cs call:

$ projinfo -s "$(projinfo -q -o WKT2:2019 --single-line --3d EPSG:32632)" -t "$(projinfo -q -o WKT2:2019 --single-line --3d EPSG:2056)" -o PROJ
Candidate operations found: 1
-------------------------------------
Operation No. 1:

unknown id, Inverse of UTM zone 32N + Inverse of CH1903+ to WGS 84 (1) + Swiss Oblique Mercator 1995, 1 m, Liechtenstein; Switzerland.

PROJ string:
+proj=pipeline
  +step +inv +proj=utm +zone=32 +ellps=WGS84
  +step +proj=cart +ellps=WGS84
  +step +proj=helmert +x=-674.374 +y=-15.056 +z=-405.346
  +step +inv +proj=cart +ellps=bessel
  +step +proj=somerc +lat_0=46.9524055555556 +lon_0=7.43958333333333 +k_0=1
        +x_0=2600000 +y_0=1200000 +ellps=bessel

I believe the logic from projinfo can be reused in cs2cs. Exercise left for the reader :-)

[hernando]

Thanks for the reply.
In fact, I was having the secret hope of being able to extend the semantics of osgeo::proj::io::createFromUserInput by motivating it through a use case in cs2cs. From your answer I infer that don't want to fiddle with the behavior of that function.

The use case I have in mind is labeling data always with a 3D CRS. For that I would like to express 2D CRSs promoted to 3D with a syntax accepted by createFomUserInput. The only current solution is taking the WKT string of the 3D promoted system, but then, the authority code is lost, which I would like to keep, hence the suggestion for Auth:code+ellipsoidal.
Having a vertical CRS that expresses ellipsoidal heights and making a compound would be another alternative, but as far as I understand from the WKT2 standard, ellipsoidal height axis are not allowed in vertical CRS (I guess because then the CRS datum is not fully defined).

[kbevers]

but then, the authority code is lost,

Yeah, cause now you've made it into a thing the authority didn't intend you to do :-) To me it sounds like this is a problem that needs to be solved upstream by registering proper 3D versions of the relevant systems.

[busstoptaktik]

this is a problem that needs to be solved upstream by registering proper 3D versions of the relevant systems.

... or, @hernando - perhaps just promote, and live with, the fact that transformations are more fundamental than CRS: There's always a multitude of highways between any two CRS, but if you want to specify a specific one, then refer to the transformation specifier, rather than trying to force your software system to guess what you mean with nothing but two CRS as hints. The latter is a nicety for low accuracy work - the former is a necessity for high accuracy work

[hernando]

I see that I'll have to look for another solution.

Yeah, cause now you've made it into a thing the authority didn't intend you to do

I get what you mean by that (for example, promoting OSGB36 doesn't work as expected as I mentioned above). However, I don't think the absence of a 3D system always implies it's invalidity, but I think it's a consequence of the EPSG data model not handling the combinatorial explosion very well, so many valid variants are omitted for pragmatic reasons.
Specially, I think it's too optimistic to expect local geodetic authorities to register a 3D variant for all the already existing 2D CRSs that lack it, specially projected ones (the same way I don't expect anybody to submit all UTM projections for all the realizations of WGS84 and ETRS89 ensembles).
Taking PPK in drone imaging as an example, a third party application could allow users in the US to use a known state plan as output CRS, while keeping the vertical axis as ellipsoidal. I don't know if this is common, but I don't see anything fundamentally wrong with it either.

perhaps just promote, and live with, the fact that transformations are more fundamental than CRS

That's on our radar, but much further down the road. Also, we've noticed that many users have troubles understanding what are the choices that they have to make, so while some power users will benefit from specifying transformations rather that picking codes, most of the users would still prefer to just pick a CRS from a list.
Given this, I would like to exchange 3D data between applications in the same ecosystem without implicit assumptions about the vertical component, which enforces using a explicit 3D descriptions of the CRS. And at the same time, I would like is to be able to display the code that was chosen by the user despite storing a 3D CRS internally.

[hernando]

Just in case it was missed, let me repeat that I'm offering to implement this extension (the same way I would have implemented --3d in cs2cs if it was solving all my problems).

[kbevers]

However, I don't think the absence of a 3D system always implies it's invalidity, but I think it's a consequence of the EPSG data model not handling the combinatorial explosion very well, so many valid variants are omitted for pragmatic reasons.
Specially, I think it's too optimistic to expect local geodetic authorities to register a 3D variant for all the already existing 2D CRSs that lack it, specially projected ones

I certainly agree that the EPSG database and the standards it is based on is limiting. @busstoptaktik and others are actively working towards fixing that. In the meantime I don't think it is a good idea to make our own extensions to the standards. That'll just make things even more messy. So for now it is up to those of us who work in geodetic authorities to register usable CRSs and transformations. And if we are lucky some day we may have standars that allow us to do things smarter.

Just in case it was missed, let me repeat that I'm offering to implement this extension (the same way I would have implemented --3d in cs2cs if it was solving all my problems).

That much was clear. I appreciate the initiative. I am not convinced yet that this is a good idea though. I'd like to hear @rouault's opinion on this too.

[rouault]

The issue with EPSG:XXXX+ellipsoidal is that it is not obvious which unit will be used for the vertical axis (same issue as with --3d).

There's a less sexy (ok let's be honest: completely user unfriendly), but more precise way, of creating Projected 3D CRS using a not very well known syntax of OGC CRS URNs ( OGC 07-092r2: para 7.5.4 "URN combined references for projected or derived CRSs"), which is already implemented (using a PROJ:ENh custom coordinate system that was included in the database for that purpose):

$ bin/projinfo urn:ogc:def:crs,crs:EPSG::4979,cs:PROJ::ENh,coordinateOperation:EPSG::16031
PROJ.4 string:
+proj=utm +zone=31 +datum=WGS84 +units=m +no_defs +type=crs

WKT2:2019 string:
PROJCRS["WGS 84 (3D) / UTM zone 31N",
    BASEGEOGCRS["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]],
        ID["EPSG",4979]],
    CONVERSION["UTM zone 31N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",3,
            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",16031]],
    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]]]]

[hernando]

The issue with EPSG:XXXX+ellipsoidal is that it is not obvious which unit will be used for the vertical axis (same issue as with --3d).

Yes, that's my concern as well. If it wasn't for this detail, I would have directly submitted a PR.

There's a less sexy (ok let's be honest: completely user unfriendly), but more precise way, of creating Projected 3D CRS using a not very well known syntax of OGC CRS URNs ( OGC 07-092r2: para 7.5.4 "URN combined references for projected or derived CRSs"), which is already implemented (using a PROJ:ENh custom coordinate system that was included in the database for that purpose):

$ bin/projinfo urn:ogc:def:crs,crs:EPSG::4979,cs:PROJ::ENh,coordinateOperation:EPSG::16031

This is quite interesting, where can I find the list of valid coodinate operations, directly in the database?
As for PROJ:ENh, it seems that it would only work for heights in meters and right-handed CRSs (i.e. always Easting, Northing), For other cases, I would need to add additional custom CS definitions to proj.db, right?

(edited, I wrote wrongly right-handed as Northing-Easting)

[rouault]

where can I find the list of valid coodinate operations, directly in the database?

this is the conversion table, which is referenced by entries in projected_crs.

As for PROJ:ENh, it seems that it would only work for heights in meters and right-handed CRSs (i.e. always Northing, Easting), For other cases, I would need to add additional custom CS definitions to proj.db, right?

yes, in the coordinate_system and axis tables

[hernando]

I think this ticket can be closed. @rouault, thanks for adding --3d to cs2cs.