[gdal-dev] Compound CRS

[Norman Barker]

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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20260423/01d040b5/attachment.htm>