Finding 2D CRS from 3D to consistently build a compound horizontal + vertical CRS

[rhugonnet]

Hi @snowman2 and community!

I'm having issues consistently building a compound horizontal + vertical CRS in an automated manner for any DEM, which involves converting a 3D CRS to 2D. The horizontal CRS is known from the metadata, while the vertical CRS can sometimes be user-input based.

Why the issue arose
3D CRS seem to now be required to transform the z coordinate (since pyproj 3.3/PROJ 8? I didn't run into issues before and couldn't find detailed notes on that).
For instance, converting z from CRS("EPSG:4326") does not work anymore, but it does from CRS("EPSG:4979") which is also CRS("EPSG:4326").to_3d().

If I understood the several linked issues that I read (e.g., #630, https://gis.stackexchange.com/questions/436401/transforming-epsg4979-to-given-wkt-epsg6654), this is now the expected behaviour?

Issue
Having 3D required for z transformation is fine by me. However, when trying to build consistently a CompoundCRS from a VerticalCRS or a BoundCRS, 3D CRS are not accepted as part of the compounds:

CRS("EPSG:4979+5773")

pyproj.exceptions.CRSError: Invalid projection: EPSG:4979+5773: (Internal Proj Error: proj_create: components of the compound CRS do not belong to one of the allowed combinations of http://docs.opengeospatial.org/as/18-005r5/18-005r5.html#34)

while 2D work:

CRS("EPSG:4326+5773")

I need to be able to find the 2D version of a CRS (e.g., CRS("EPSG:4326")) from its 3D version (e.g., CRS("EPSG:4979")) in an automated manner, to build a CompoundCRS from any 2D version of a CRS and any VerticalCRS or BoundCRS.

But I could not find a CRS.to_2d() function or equivalent. Is there any that I didn't find? Or maybe an obvious solution that I missed? 😅
Thanks for the help! :)

[snowman2]

this is now the expected behaviour?

Yes, you have to explicitly specify a 2D or 3D CRS.

But I could not find a CRS.to_2d() function or equivalent. Is there any that I didn't find? Or maybe an obvious solution that I missed?

Normally, I would suggest using CRS.sub_crs_list, but that doesn't exist for the 3D Geographic CRS. I don't think there s a to_2d function at the moment.

@rouault do you have any suggestions?

[rhugonnet]

Amazing, thanks! Do you need a hand wrapping that into pyproj? I could give it a try right now.

[snowman2]

That would be great!

[snowman2]

Related: OSGeo/PROJ#3119

[rhugonnet]

For the record (if someone sees this discussion): we wrapped CRS.to_2d() in pyproj in #1267.