Handling of compound CRS

[jorisvandenbossche]

One of the type of CRSses is a "compound CRS", which is a combination of one CRS for the horizontal coordinates and a vertical CRS for the z coordinate.

An example is EPSG:5972 (http://epsg.io/5972, taken from https://github.com/OSGeo/proj.4/blob/master/data/sql/compound_crs.sql):

In [9]: crs = pyproj.CRS.from_epsg(5972)

In [10]: crs     
Out[10]: 
<CRS: epsg:5972>
Name: ETRS89 / UTM zone 32 + NN2000 height
Ellipsoid:
- semi_major_metre: 6378137.00
- semi_minor_metre: 6356752.31
- inverse_flattening: 298.26
Area of Use:
- name: Norway - onshore - 6°E to 12°E
- bounds: (6.0, 57.93, 12.0, 65.76)
Prime Meridian:
- longitude: 0.0000
- unit_name: degree
- unit_conversion_factor: 0.01745329
Axis Info:
- UNDEFINED

In [11]: print(crs.to_wkt())  
COMPOUNDCRS["ETRS89 / UTM zone 32 + NN2000 height",
    PROJCRS["ETRS89 / UTM zone 32N",
        BASEGEOGCRS["ETRS89",
            DATUM["European Terrestrial Reference System 1989",
                ELLIPSOID["GRS 1980",6378137,298.257222101,
                    LENGTHUNIT["metre",1]]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["degree",0.0174532925199433]]],
        CONVERSION["UTM zone 32N",
            METHOD["Transverse Mercator",
                ID["EPSG",9807]],
            PARAMETER["Latitude of natural origin",0,
                ANGLEUNIT["degree",0.0174532925199433],
                ID["EPSG",8801]],
            PARAMETER["Longitude of natural origin",9,
                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]]],
        CS[Cartesian,2],
            AXIS["(E)",east,
                ORDER[1],
                LENGTHUNIT["metre",1]],
            AXIS["(N)",north,
                ORDER[2],
                LENGTHUNIT["metre",1]]],
    VERTCRS["NN2000 height",
        VDATUM["Norway Normal Null 2000"],
        CS[vertical,1],
            AXIS["gravity-related height (H)",up,
                LENGTHUNIT["metre",1]]],
    USAGE[
        SCOPE["unknown"],
        AREA["Norway - onshore - 6°E to 12°E"],
        BBOX[57.93,6,65.76,12]],
    ID["EPSG",5972]]

In [12]: crs.is_projected        
Out[12]: False

So this does not fully fit in the current CRS class:

• "Axis Info" in the repr is undefined, although we could get that by combining the axis info of the projected and vertical crs
• do we want to expose functionality to get the two CRSses? (accessible with proj_crs_get_sub_crs) Or is that too much detail for end users, and is it fine enough with having the current API?
• We might want to have is_projected return True if the horizontal CRS is projected?

(some similar questions can be raised for the Bound CRS)

[snowman2]

I think adding a sub_crs_list propery to return CRS instances based on proj_crs_get_sub_crs is warranted for this use case.

[jorisvandenbossche]

We might want to have is_projected return True if the horizontal CRS is projected?

This is also what the GDAL "is projected" function does, it special cases for CompoundCRS:

https://github.com/OSGeo/gdal/blob/60ab2831abf9446c629b3a417a16cb35e7f41e63/gdal/ogr/ogrspatialreference.cpp#L7553-L7581

[snowman2]

I am going to close this as the solution is merged into master. Thanks!