[gdal-dev] Compound CRS

[Even Rouault]

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 list
> gdal-dev at lists.osgeo.org
> https://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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20260424/4a3b49a3/attachment-0001.htm>