Converting geographic coordinates into tangent-plane Cartesian with GDAL

The tangent plane CRS is known as a Topocentric CRS.

The conversion formulae from Geocentric, and from Geographic coordinates to Topocentric coordinates are defined in the Section 4.1.2, of the IOGP document: Geomatics Guidance Note 7, part 2.

The EPSG registry has two example definitions of Topocentric CRS: EPSG:5819 and EPSG:5820.

For instance, we can see the EPSG:5819 WKT2:2019 definition with:

C:>projinfo EPSG:5819
PROJ.4 string:
Error when exporting to PROJ string: Unsupported conversion method: Geographic/topocentric conversions

WKT2:2019 string:
PROJCRS["EPSG topocentric example A",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4979]],
    CONVERSION["EPSG topocentric example A",
        METHOD["Geographic/topocentric conversions",
            ID["EPSG",9837]],
        PARAMETER["Latitude of topocentric origin",55,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8834]],
        PARAMETER["Longitude of topocentric origin",5,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8835]],
        PARAMETER["Ellipsoidal height of topocentric origin",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8836]]],
    CS[Cartesian,3],
        AXIS["topocentric East (U)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["topocentric North (V)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
        AXIS["topocentric height (W)",up,
            ORDER[3],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["unknown"],
        AREA["To be specified"],
        BBOX[-90,-180,90,180]],
    ID["EPSG",5819]]  

But you can see there also the problem: The conversion method to Topocentric coordinates is not listed in the available PROJ Transformations.

I don't know if it can be assumed as a 3D affine transformation. In that case all the parameters may be derived from the definition of the tangent plane.

We can see also the WKT representation just for the conversion from EPSG:4979 (Geographic 3D (lat, lon, h)) to EPSG:5819 Topocentric CRS, with:

C:>projinfo -s EPSG:4979 -t EPSG:5819
Candidate operations found: 1
-------------------------------------
Operation n┬░1:

EPSG:15594, EPSG topocentric example A, 0 m, Not specified

PROJ string:
Error when exporting to PROJ string: Unsupported conversion method: Geographic/topocentric conversions

WKT2:2019 string:
CONVERSION["EPSG topocentric example A",
    METHOD["Geographic/topocentric conversions",
        ID["EPSG",9837]],
    PARAMETER["Latitude of topocentric origin",55,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8834]],
    PARAMETER["Longitude of topocentric origin",5,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8835]],
    PARAMETER["Ellipsoidal height of topocentric origin",0,
        LENGTHUNIT["metre",1],
        ID["EPSG",8836]],
    ID["EPSG",15594]]

As before, there is not a pipeline string for the operation, so I think that GDAL will fail trying to transform the data to that CRS.

---

But, if you are able to work with the geodetic coordinates as lists or arrays, you can use GeographicLib to convert them to a local cartesian system.

From C++, you can use the Geocentric and LocalCartesian classes.

From Python I think that it is not implemented in the geographiclib package (I'm not sure), but you can call the Cart Convert utility program from the subprocess module.