[gdal-dev] Compound CRS

[Javier Jimenez Shaw]

I don't think you need two level nesting derivations. You can use a
topocentric "normal" system and add the affine transformation as a derived
projected.
There are two types of topocentric : geographic_topocentric
and geocentric_topocentric.
https://github.com/OSGeo/PROJ/blob/master/test/unit/test_operation.cpp#L4484
https://github.com/OSGeo/PROJ/blob/master/test/unit/test_operation.cpp#L4544

On Fri, 24 Apr 2026 at 02:39, Even Rouault via gdal-dev <
gdal-dev at lists.osgeo.org> wrote:

> Norman,
>
> (ISO-19111 pedantically speaking your object is not a CompoundCRS, which
> would be a geographic/projected 2D + vertical CRS, but more a DerivedCRS)
>
> the error is "expected" given that it is generally not possible to infer a
> CRS definition from a PROJ pipeline
>
> 2 ways of addressing this:
>
> - the hard one: create a "real" CRS fighting with WKT/PROJSON (maybe just
> PROJJSON as I'm not sure WKT allows that 2-level nesting of derivation) of
> a DerivedGeodeticCRS (affine transform) of a DerivedGeodeticCRS (geographic
> to topocentric) of a geographic CRS... or something like that. I don't
> have  an example of that handy.  Some potential inspiration at
> https://github.com/qgis/QGIS/pull/65873/changes#diff-443e243a4b840a314841dfd4d6aae7f0981d64233badade7af0a769e29536982R3332
> and
> https://github.com/OSGeo/PROJ/blob/master/test/unit/test_operationfactory.cpp#L11174
>
> - or likely easier, write the full PROJ pipeline of the transformation,
> which should probably be just a matter of inverting your pipeline
> (reversing the steps with a +inv keyword on each) and use
> ct.SetOperation(proj_pipeline) as in
> https://github.com/OSGeo/gdal/blob/0bb1c88ff756f174577c6056af6d9452942a5a1a/autotest/osr/osr_ct.py#L359
>
> Even
> Le 24/04/2026 à 02:08, Norman Barker via gdal-dev a écrit :
>
> Hi,
>
> I am trying to use a compound CRS in GDAL. I have a CRS that is
> essentially topocentric but with an angle and skew applied to the
> tangential plane by two unit vectors uiax and uiay.
>
> I am doing this in Python and the code is as follows;
>
>             # planar angle between unit vector and northing in topocentric
>             theta =  np.arccos(np.clip(np.dot(uiay[:2], [0, 1]), -1.0, 1.0
> ))
>             # planar angle between unit vector and easting in topocentric
>             phi =  np.arccos(np.clip(np.dot(uiax[:2], [1, 0]), -1.0, 1.0))
>             # Build rotation matrix
>             rot = np.array([
>                 [np.cos(theta), -np.sin(theta), 0.0, 0.0],
>                 [np.sin(theta),  np.cos(theta), 0.0, 0.0],
>                 [          0.0,            0.0, 1.0, 0.0],
>                 [          0.0,            0.0, 0.0, 1.0],
>             ]
>             )
>             # Build shear/skew matrix
>             m = np.tan(phi)
>             skew = np.array([
>                 [1.0, 0.0, 0.0, 0.0],
>                 [  m, 1.0, 0.0, 0.0],
>                 [0.0, 0.0, 1.0, 0.0],
>                 [0.0, 0.0, 0.0, 1.0],
>             ]
>             )
>             # get affine transform
>             a = rot @ skew
>             wgs84 = osr.SpatialReference(epsg=4326)
>             wgs84.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
>             topo = osr.SpatialReference(f"+proj=pipeline "
>                                         f"+step +proj=topocentric +X_0={
> ref_ecf[0]} +Y_0={ref_ecf[1]} +Z_0={ref_ecf[2]} "
>                                         f"+step +proj=affine "
>                                         f"+s11={a[0,0]} +s12={a[0,1]}
> +s13={a[0,2]} "
>                                         f"+s21={a[1,0]} +s22={a[1,1]}
> +s23={a[1,2]} "
>                                         f"+s31={a[2,0]} +s32={a[2,1]}
> +s33={a[2,2]} "
>                                         )
>             ct = osr.CoordinateTransformation(topo, wgs84)
>
> I get the error message "ERROR 1: PROJ: proj_crs_get_coordinate_system:
> Object is not a SingleCRS"
>
> Which is true, it is not. I have grep'd the tests and have not seen any
> examples. Using `CoordinateTransformationOptions` applies once the
> conversion to lat/lon is complete.
>
> Thanks,
>
> Norman
>
> _______________________________________________
> gdal-dev mailing listgdal-dev at lists.osgeo.orghttps://lists.osgeo.org/mailman/listinfo/gdal-dev
>
> -- http://www.spatialys.com
> My software is free, but my time generally not.
> Highly recommend OxiGDAL if you want to live in the 21th century and cure Bixonimania
>
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20260424/52490e50/attachment.htm>