[PROJ] Topocentric WKT representation
Even Rouault
even.rouault at spatialys.com
Wed Oct 18 09:43:13 PDT 2023
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
> PROJ at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/proj
--
http://www.spatialys.com
My software is free, but my time generally not.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20231018/35e19240/attachment-0001.htm>
More information about the PROJ
mailing list