[PROJ] Topocentric WKT representation

[Erixen Cruz]

Hello,

I think there may be a bug with topocentric WKT representation in PROJ, but I want to confirm before making an issue in the repo.

We have data in a east-north-up coordinate system with origin on the surface of WGS84 ellipsoid. I believe the way to represent that in WKT is with a topocentric conversion, which is codified in EPSG:9836 https://epsg.org/coord-operation-method_9836/Geocentric-topocentric-conversions.html. PROJ implements this since version 8 https://proj.org/en/9.3/operations/conversions/topocentric.html.

I build latest PROJ master locally and use projinfo to get WKT for a CRS defined by topocentric conversion at origin -3982059.42, 3331314.88, 3692463.58 (EPSG:4978).

$ PROJ_DATA=~/Documents/OSGeo/PROJ/build/data/ ~/Documents/OSGeo/PROJ/build/bin/projinfo -o WKT2:2019 "+proj=topocentric +ellps=WGS84 +X_0=-3982059.42 +Y_0=3331314.88 +Z_0=3692463.58 +no_defs +type=crs"
PROJCRS["unknown",
    BASEGEOGCRS["unknown",
        DATUM["Unknown based on WGS 84 ellipsoid",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1],
                ID["EPSG",7030]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8901]]],
    CONVERSION["unknown",
        METHOD["Geocentric/topocentric conversions",
            ID["EPSG",9836]],
        PARAMETER["Geocentric X of topocentric origin",-3982059.42,
            LENGTHUNIT["metre",1],
            ID["EPSG",8837]],
        PARAMETER["Geocentric Y of topocentric origin",3331314.88,
            LENGTHUNIT["metre",1],
            ID["EPSG",8838]],
        PARAMETER["Geocentric Z of topocentric origin",3692463.58,
            LENGTHUNIT["metre",1],
            ID["EPSG",8839]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]

This is peculiar because the base CRS is geographic (4979?) instead of geodetic 4978. I feed this to cs2cs to see if I can go back from the origin 0, 0, 0 to 4978.

$ echo 0 0 0 | PROJ_DATA=~/Documents/OSGeo/PROJ/build/data/ ~/Documents/OSGeo/PROJ/build/bin/cs2cs "PROJCRS[\"unknown\",BASEGEOGCRS[\"unknown\",DATUM[\"Unknown based on WGS 84 ellipsoid\",ELLIPSOID[\"WGS 84\",6378137,298.257223563,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",7030]]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8901]]],CONVERSION[\"unknown\",METHOD[\"Geocentric/topocentric conversions\",ID[\"EPSG\",9836]],PARAMETER[\"Geocentric X of topocentric origin\",-3982059.42,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8837]],PARAMETER[\"Geocentric Y of topocentric origin\",3331314.88,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8838]],PARAMETER[\"Geocentric Z of topocentric origin\",3692463.58,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8839]]],CS[Cartesian,2],AXIS[\"(E)\",east,ORDER[1],LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]],AXIS[\"(N)\",north,ORDER[2],LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]]]" +to EPSG:4978
pipeline: Pipeline: Mismatched units between step 1 and 2
*     * inf

And the peculiarity reveals itself: the topocentric conversion goes to cartesian (geodetic), not geographic. So there is units mismatch meters and degrees/radians. When I replace the BASEGEOGCRS with BASEGEODCRS with 4978 definition, it works as expected:

$ echo 0 0 0 | PROJ_DATA=~/Documents/OSGeo/PROJ/build/data/ ~/Documents/OSGeo/PROJ/build/bin/cs2cs "PROJCRS[\"unknown\",BASEGEODCRS[\"WGS 84\",ENSEMBLE[\"World Geodetic System 1984 ensemble\",MEMBER[\"World Geodetic System 1984 (Transit)\"],MEMBER[\"World Geodetic System 1984 (G730)\"],MEMBER[\"World Geodetic System 1984 (G873)\"],MEMBER[\"World Geodetic System 1984 (G1150)\"],MEMBER[\"World Geodetic System 1984 (G1674)\"],MEMBER[\"World Geodetic System 1984 (G1762)\"],MEMBER[\"World Geodetic System 1984 (G2139)\"],ELLIPSOID[\"WGS 84\",6378137,298.257223563,LENGTHUNIT[\"metre\",1]],ENSEMBLEACCURACY[2.0]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433]],CS[Cartesian,3],AXIS[\"(X)\",geocentricX,ORDER[1],LENGTHUNIT[\"metre\",1]],AXIS[\"(Y)\",geocentricY,ORDER[2],LENGTHUNIT[\"metre\",1]],AXIS[\"(Z)\",geocentricZ,ORDER[3],LENGTHUNIT[\"metre\",1]],USAGE[SCOPE[\"Geodesy. Navigation and positioning using GPS satellite system.\"],AREA[\"World.\"],BBOX[-90,-180,90,180]],ID[\"EPSG\",4978]],CONVERSION[\"unknown\",METHOD[\"Geocentric/topocentric conversions\",ID[\"EPSG\",9836]],PARAMETER[\"Geocentric X of topocentric origin\",-3982059.42,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8837]],PARAMETER[\"Geocentric Y of topocentric origin\",3331314.88,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8838]],PARAMETER[\"Geocentric Z of topocentric origin\",3692463.58,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8839]]],CS[Cartesian,2],AXIS[\"(E)\",east,ORDER[1],LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]],AXIS[\"(N)\",north,ORDER[2],LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]]]" +to EPSG:4978
-3982059.42 3331314.88 3692463.58

Am I misusing projinfo and/or the topocentric conversion? Or should PROJ be returning a BASEGEODCRS instead of BASEGEOGCRS in this case? Thanks.

Sincerely, Erixen
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20231018/32b95fd6/attachment-0001.htm>