Good analysis. Fixed per https://github.com/OSGeo/PROJ/pull/3924
Le 18/10/2023 à 16:23, Erixen Cruz via PROJ a écrit :
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
--
http://www.spatialys.com
My software is free, but my time generally not.
_______________________________________________
PROJ mailing list
[email protected]
https://lists.osgeo.org/mailman/listinfo/proj
