Add vertical crs option for projects

[nyalldawson]

(refs qgis/QGIS-Enhancement-Proposals#267)

If the project CRS is a compound CRS, then the vertical CRS for the project will be the vertical component of the main project CRS. Otherwise it will be the value explicitly set by the user.

Users can specify the vertical crs through the Project Properties, Elevation tab:

Eg. when the project CRS is set to a compound crs, the user will see:

When the project CRS is a horizontal CRS, then the user can select a specific vertical CRS for their project:

Currently, this is a "metadata" type property only, and has no impact on rendering or handling of features. (It IS exposed through new @project_vertical_crs_xxx variables, though)

[qgis-bot]

@nyalldawson
This pull request has been tagged as requiring documentation.

A documentation ticket will be opened at https://github.com/qgis/QGIS-Documentation when this PR is merged.

Please update the description (not the comments) with helpful description and screenshot to help the work from documentors.
Also, any commit having [needs-doc] or [Needs Documentation] in will see its message pushed to the issue, so please be as verbose as you can.

Thank you!

[qgis-bot]

@nyalldawson

This pull request has been tagged for the changelog.

• The description will be harvested so please provide a "nearly-ready" text for the final changelog
• If possible, add a nice illustration of the feature. Only the first one in the description will be harvested (GIF accepted as well)
• If you can, it's better to give credits to your sponsor, see below for different formats.

You can edit the description.

Format available for credits

• Funded by NAME
• Funded by URL
• Funded by NAME URL
• Sponsored by NAME
• Sponsored by URL
• Sponsored by NAME URL

Thank you!

[nyalldawson]

@rouault @wonder-sk @uclaros this relates to recent discussions

[github-actions]

🪟 Windows builds ready!

Windows builds of this PR are available for testing here. Debug symbols for this build are available here.

(Built from commit f2f6513)

[wonder-sk]

@nyalldawson one thing that I am wondering about is whether it is a good design choice to have "main CRS" and "vertical CRS" handled separately. It definitely makes sense when using 2D "main" CRS and adding vertical CRS to it. But what about Geographic3d or Geocentric CRS - in those cases, a separate vertical CRS does not seem to make sense...? (just checked, for EPSG:4978 or EPSG:4979 the verticalCrs() call returns invalid object)

Another issue is that any code dealing with elevations may get a bit complicated, as it may need to first consult project's CRS with crs(), and then based on the CRS type, it may also need to check verticalCrs().

Related to the separate storage of "main" and "vertical" CRS in projects is that later, probably every other place that deals with 3D CRS would probably need to adopt this kind of separate storage, just to make things consistent?

My preference would be to have the whole CRS stored in one place in the project, and have just setCrs() and crs() set/get - or are there some problems connected to that? We could still have verticalCrs() getter, but users would need to understand it can still return invalid CRS even for 3D CRS...

As for GUI, have you thought about having the vertical CRS configuration embedded in the main CRS selection widget? To me it feels that it belongs there, especially due to the complications from different types of CRS. And also we would not need to worry about dedicated vertical CRS selectors in places where they may be needed (e.g. map layer CRS override, CRS for data export...)

[rouault]

But what about Geographic3d or Geocentric CRS - in those cases, a separate vertical CRS does not seem to make sense...?

No, it doesn't. A vertical CRS can only be adjoined to a 2D geographic or projected CRS. A Vertical CRS is about gravity related height and is an "autonomous" object. For example the 3rd axis in a geographic 3D CRS with the ellipsoidal height is different. An ellipsoidal height has no physical meaning without being associated with the (x,y) parts of the coordinate. Whereas for a gravity related height, you know that if H1 > H2, whatever their horizontal coordinates, water will flow from 1 to 2 in the absence of obstacles.

My preference would be to have the whole CRS stored in one place in the project, and have just setCrs() and crs() set/get - or are there some problems connected to that?

A "whole CRS" with a horizontal and vertical CRS would then match the definition of a compound CRS, and would be the "ISO-19999 compliant" / clean way of dealing with the issue. But I suspect a lot of code isn't ready to deal with compound CRS. Typically most compound CRS do not have a well-known EPSG identifier, and thus the separation between both is a pragmatic way of dealing with the issue.
The alternative would be to audit the whole code base and replace a number of crs() calls by crs().demoteTo2d() where needed. demoteTo2d() would be a sort of generalized horizontalCrs() (there's such as PROJ function). For compound CRS, it returns the horizontal part. For geographic / projected 3D CRS, it returns the corresponding 2D CRS. And for geocentric CRS, it returns the source object unmodified (at least that's the PROJ implementation).

[nyalldawson]

@wonder-sk @rouault my thinking was that I'd add a "crs3d()" getter, which would always return a compound crs of the crs() +verticalCrs(), or just the crs() if no vertical crs is available. Then we'd ensure that in all cases where elevation is important (3d, elevation profile,... ) we create the transform using crs3d(), not crs()...

Does that answer your concerns?

[rouault]

a "crs3d()" getter

on QgsProject? sounds to me to be a reasonable way forward

[nyalldawson]

@rouault

on QgsProject?

Yes, that's correct. Then it's all opt in, and with a low risk of regressions/soft API breakage.

[wonder-sk]

my thinking was that I'd add a "crs3d()" getter, which would always return a compound crs of the crs() +verticalCrs(), or just the crs() if no vertical crs is available. Then we'd ensure that in all cases where elevation is important (3d, elevation profile,... ) we create the transform using crs3d(), not crs()...

Do you foresee problems if we only had crs() (without crs3d()), which could return either 2d or 3d CRS? That would be simpler to use by API users... My hope was that QgsCoordinateReferenceSystem abstracts the differences well, and for large majority of cases where the crs() API is used, it does not really matter whether a CRS is 2d or 3d, but I may be terribly wrong 😄

[nyalldawson]

@wonder-sk

Do you foresee problems if we only had crs() (without crs3d()), which could return either 2d or 3d CRS? That would be simpler to use by API users... My hope was that QgsCoordinateReferenceSystem abstracts the differences well, and for large majority of cases where the crs() API is used, it does not really matter whether a CRS is 2d or 3d, but I may be terribly wrong

Honestly, I don't know what the impact of that would be. It might be issue free, but it also might cause a bunch of regressions. Possibly issues could be:

• use of .authid() (or even worse, the .srsid() / .srid() ) methods on crs(). Those won't return any meaningful value for compound crs. So if somewhere in qgis code we are storing authid as a way of persisting crs then that will break. (And similar for plugins). We do have methods for writing / reading crs to xml which shouldn't have issues, but I'm not 100% sure if they are used everywhere in master for serialisation of crs. (And I'd be confident that plugins aren't correctly using these methods 🤣 )
• possible speed regressions in transforms -- are transforms involving compound crs slower then their purely horizontal counterparts? (@rouault ?). If so, then we risk unwanted performance impacts on purely 2d use cases, including the all-important 2d map rendering.
• are there issues when we use a compound crs for transforms and DON'T have valid z values in the transformed points?

There's just a lot of unknowns which make me prefer to take a safe incremental approach here... Maybe we could use an opt-in crs3d() initially, and then revisit for 4.0? Or even in a few releases time make crs() / crs3d() aliases for each other and deprecated crs3d()?

[rouault]

• use of .authid() (or even worse, the .srsid() / .srid() ) methods on crs(). Those won't return any meaningful value for compound crs

That's certainly a concern. Users could be easily concerned by having a "unknown CRS" (they are already for 2D CRS that fail to identify to a EPSG code). I'd also say based on my experience with GDAL drivers that a lot of them could have issues on the write site when presented with a CompoundCRS.

possible speed regressions in transforms -- are transforms involving compound crs slower then their purely horizontal counterparts?

There are several aspects. First you'll get a 3D capable transform only if the source and target CRS are 3D (compound or Geographic/Projected 3D). Then, due to complications / trickyness within the PROJ pipeline inference logic, the 2D part of a 3D pipeline between 3D CRS might not be strictly equal to the 2D part of what you'll get using the corresponding 2D CRS. A 3D pipeline might also be slower at execution due to using a geoid grid, but when you really want to do 3D transformations, that's actually something you want.

are there issues when we use a compound crs for transforms and DON'T have valid z values in the transformed points?

If you have a 3D pipeline, and you don't have a valid Z, then you still need to put a dummy (non HUGE_VAL) value for the z, like 0, to get sensible results for the X and Y parts

There's just a lot of unknowns which make me prefer to take a safe incremental approach here... Maybe we could use an opt-in crs3d() initially, and then revisit for 4.0? Or even in a few releases time make crs() / crs3d() aliases for each other and deprecated crs3d()?

Maybe fullCrs() / crsPossibly3d() or something like that would be better than crs3d()? since it might actually return just a 2D CRS

[wonder-sk]

@nyalldawson my main question was whether there was anything wrong with having just crs() conceptually - from the look of it, it is mainly that the current implementation in core and plugins may have issues with CRS handling if we suddenly allow compound CRS in crs(). Thanks for the clarifications, I understand that addressing all the potential issues may be a lot of work if we fully switched to crs() being possibly 3D crs. I am happy with the current plan to keep crs() as is, and have crs3d() (or an alternative naming from Even) in addition to crs() that stays unchanged.

[nyalldawson]

@rouault @wonder-sk

Maybe fullCrs() / crsPossibly3d() or something like that would be better than crs3d()? since it might actually return just a 2D CRS

I could live with fullCrs(), although I'd prefer a name which better indicates "use this crs if you need z handling". What about "crsFor3d()"? or "crs3dAware()"?

[wonder-sk]

I could live with fullCrs(), although I'd prefer a name which better indicates "use this crs if you need z handling". What about "crsFor3d()"? or "crs3dAware()"?

I am quite easy... The originally suggested crs3d() would work for me too 🙂 (and I would probably prefer it for brevity)

[github-actions]

The QGIS project highly values your contribution and would love to see this work merged! Unfortunately this PR has not had any activity in the last 14 days and is being automatically marked as "stale". If you think this pull request should be merged, please check

• that all unit tests are passing

• that all comments by reviewers have been addressed

• that there is enough information for reviewers, in particular

  • link to any issues which this pull request fixes

  • add a description of workflows which this pull request fixes

  • add screenshots if applicable

• that you have written unit tests where possible
  In case you should have any uncertainty, please leave a comment and we will be happy to help you proceed with this pull request.
  If there is no further activity on this pull request, it will be closed in a week.

[qgis-bot]

@nyalldawson
A documentation ticket has been opened at qgis/QGIS-Documentation#8999
It is your responsibility to visit this ticket and add as much detail as possible for the documentation team to correctly document this change.
Thank you!