# Converting geographic coordinates into tangent-plane Cartesian with GDAL

I'm working on reading a variety of geographic raster data (primarily elevation data, land clutter type, and land clutter height) to project them into a local-tangent-plane Cartesian system (east-north-up) for running radio propagation calculations. The geodata comes in a variety of file formats and coordinate systems, all of which GDAL reads in just fine, giving me an array with the banded data and an WKT string describing the coordinates.

It looks like I can create a GDAL CoordinateTransformation object that will project from the geodata coordinates into the Cartesian system given the WKT strings for the two systems.

Is this the right approach? If so, how I do create the WKT string for the tangent-plane system?

I'm currently prototyping in Python so that I can visualize data with MatPlotLib, but will eventually do the coding in C/C++ or Fortran.

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.

• Thanks—this clarified much for me. Looks a possible next step would be deriving the WKT string for an arbitrary topocentric CRS to cover the area of interest. I see that someone asked this question a while back: gis.stackexchange.com/questions/241094/… Feb 7, 2020 at 21:36
• @AndreasYankopolus, you are welcome. I have added an answer to that question. The problem there is a mixture of Derived Engineering CRS from a Base Geodetic CRS. The Topocentric conversion methods (for Geocentric and Geographic 3D) feature request ticket in PROJ was closed by inactivity: github.com/OSGeo/PROJ/issues/500. Feb 7, 2020 at 22:44
• About the WKT string for a custom Topocentric CRS, the conversion parameters must be customized (wich customize all except the Base CRS). The parameters are just latitude, longitude and ellipsoidal height of the topocentric point. Feb 8, 2020 at 5:44
• It seems it has been recently revived by the PROJ dev team; see: github.com/OSGeo/PROJ/pull/2444 Dec 5, 2020 at 13:43