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
_______________________________________________ PROJ mailing list [email protected] https://lists.osgeo.org/mailman/listinfo/proj
