QGIS fails to understand geocentric cartesian coordinates

[kbevers]

Describe the bug

Adding a new layer from a text file with geocentric cartesian coordinates results in a layer without geometry. Attached is a file, GR96_fiducial_stations_1996.txt, with EPSG:4908 coordinates which QGIS in theory should be able to understand and transform to the project CRS. This fails.

How to Reproduce

0. Create a new project
1. Add greenland.gpkg (packaged here: greenland.zip ). Make sure that project CRS is set to EPSG:4909.
2. Add delimited text layer
3. Open the attached GR96_fiducial_stations_1996.txt
4. Put field 1 as x, 2 as y, 3 as z.
5. Set "geometry CRS" to EPSG:4908
6. Add the layer
7. Note that the coordinates are not loaded correctly. They all end up on a line where they in fact should end up spread out in and around Greenland. See the image below. The orange points are the ones QGIS has imported and transformed. The blue ones are the correct positions transformed using PROJ.

Note that the x-coordinate (longitude) seems to be represented correct, but that the y-coordinate (latitude) is set to 0. Worth mentioning is also that when zooming in on the orange points they will disappear and the identify features tool won't recognize the points (independent on zoom level).

QGIS and OS versions

QGIS version
3.10.2-A Coruña
QGIS code revision
d4cd3cfe5a
Compiled against Qt
5.11.2
Running against Qt
5.11.2
Compiled against GDAL/OGR
3.0.3
Running against GDAL/OGR
3.0.3
Compiled against GEOS
3.8.0-CAPI-1.13.1
Running against GEOS
3.8.0-CAPI-1.13.1 
Compiled against SQLite
3.29.0
Running against SQLite
3.29.0
PostgreSQL Client Version
11.5
SpatiaLite Version
4.3.0
QWT Version
6.1.3
QScintilla2 Version
2.10.8
Compiled against PROJ
6.3.0
Running against PROJ
Rel. 6.3.0, January 1st, 2020
OS Version
Windows 10 (10.0)
Active python plugins
ClusterPoints; 
DataPlotly; 
DigitizingTools; 
Kortforsyningen; 
mmqgis; 
profiletool; 
QuickWKT; 
tiss; 
VectorFieldCalc; 
db_manager; 
MetaSearch; 
processing

Additional context

Using a recent version of PROJ the correct transformation can be performed:

cs2cs -f %.4f EPSG:4908 +to EPSG:4909 GR96_fiducial_stations_1996.txt
78.9296 11.8651 78.4355  NYAL
52.1784 5.8096 96.8524  KOSG
69.6627 18.9383 132.4359  TROM
57.3953 11.9255 45.5531  ONSA
76.5373 -68.7880 54.9832  THU1
64.1388 -21.9555 93.0416  REYK
47.5952 -52.6777 152.8254  STJO
62.4809 -114.4807 180.8542  YELL
66.9874 -50.9448 229.8540  KELY
58.7591 -94.0887 -19.4488  CHUR
76.5370 -68.8250 36.0272  THU2
60.7153 -46.0478 110.3945  QAQ1

Since QGIS is using PROJ this should work. I suspect the problem has something to do with the PROJ.4 representation of geocentric cartesian CRS's which ends up being +proj=geocent ... which was not designed for that use-case.

QGIS seems to know which transformation to apply as can be seen in the image below from the "Select Datum Transformations" dialog.

which matches the output from projinfo:

 projinfo -s EPSG:4908 -t EPSG:4909
Candidate operations found: 1
-------------------------------------
Operation n┬░1:

unknown id, Conversion from GR96 (geocentric) to GR96 (geog3D), 0 m, unknown domain of validity

PROJ string:
+proj=pipeline +step +inv +proj=cart +ellps=GRS80 +step +proj=unitconvert +xy_in=rad +z_in=m +xy_out=deg +z_out=m +step +proj=axisswap +order=2,1

WKT2_2018 string:
CONVERSION["Conversion from GR96 (geocentric) to GR96 (geog3D)",
    METHOD["Geographic/geocentric conversions",
        ID["EPSG",9602]]]

[nirvn]

Can you try with latest QGIS 3.10.3? There was a fair amount of fixes that could resolve your issue.

[kbevers]

Can you try with latest QGIS 3.10.3? There was a fair amount of fixes that could resolve your issue.

Probably not. At least not immediately (coorporate PC, you know). I may be able to give it a go with 3.12.0. Hang on :-)

[gioman]

2\. (packaged here: [greenland.zip].(https://github.com/qgis/QGIS/files/4281693/greenland.zip))

@kbevers this attachment/zip is empty.

[kbevers]

@kbevers this attachment/zip is empty.

Whoops. Try this one: greenland.zip

[gioman]

4\. Open the attached [GR96_fiducial_stations_1996.txt](https://github.com/qgis/QGIS/files/4281622/GR96_fiducial_stations_1996.txt)

@kbevers also this is not really a CSV (no clear seprator, is not space and also not TAB), can you attach a real CSV? Thanks!

[kbevers]

In 3.12.0 I get this error:

No transform is available between EPSG:4908 - GR96 and EPSG:4909 - GR96.
No such file or directory

and the "Select Datum Transformations" dialog doesn't show any transformations between EPSG:4908 and EPSG:4909. At least QGIS 3.12 aren't making promises it can't keep but it would be nice to be able to have support for geocentric cartesian coordinates though.

[kbevers]

@kbevers also this is not really a CSV (no clear seprator, is not space and also not TAB), can you attach a real CSV? Thanks!

Sure: GR96_fiducial_stations_1996_csv.txt

although whitespace is the clear separator in the other file. I have not problems loading it correctly.

[gioman]

although whitespace is the clear separator in the other file. I have not problems loading it correctly.

not here. Anyway I also tested QGIS 3.10.3 and 3.12 and confirm the reported results. The problem seems to be the reprojection of the data because if the project CRS is the same as the point data CRS then is ok.

[kbevers]

The problem seems to be the reprojection of the data because if the project CRS is the same as the point data CRS then is ok.

Trying that I at least get something that looks kinda correct. North has moved quite a bit west though :)

It doesn't really make sense to put geocentric cartesian coordinates on a 2D canvas. Which is what is done in this case. Basically the x and y coordinates are used no questions asked, which results in weird behavior since the origin of those coordinates are the center of the earth. See https://en.wikipedia.org/wiki/ECEF for more info.

[Pedro-Murteira]

I can still replicate this behaviour on QGIS 3.16.15 and 3.22.2.

[StefanoVagelli]

Hello, I'm using QGIS 3.28.1 and trying to change CRS yo a shapefile from EPSG:25832 (geographical) to EPSG:6704 (geocentric), but the result is rotated 90 degrees, see image (below is after trasformation). There a workaround ? Thank you.