[PROJ] [Questions] Vertical transformations advices for PROJ >= 9.1 version

Wed Oct 16 08:15:37 PDT 2024

Thanks for the answer, yeah i see I think that we might have too strict thresholds for our tests.

Thank you for your time

________________________________
De : Even Rouault <even.rouault at spatialys.com>
Envoyé : mercredi 16 octobre 2024 09:43:35
À : Christian DOULIAC; proj at lists.osgeo.org
Objet : Re: [PROJ] [Questions] Vertical transformations advices for PROJ >= 9.1 version

Unless you are interested in sub-millimetric precision, you shouldn't care about the difference between GRS80 and WGS84:

$ echo 4189881.02 146313.87 4790558.75 | bin/cct -d 4 +proj=pipeline  +step +inv +proj=cart +ellps=WGS84  +step +proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000 +y_0=6600000 +ellps=GRS80
  626830.1633   6878250.0291        0.0017           inf

$ echo 4189881.02 146313.87 4790558.75 | bin/cct -d 4 +proj=pipeline  +step +inv +proj=cart +ellps=GRS80  +step +proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000 +y_0=6600000 +ellps=GRS80
  626830.1633   6878250.0292        0.0017           inf

And if you care about sub-millimetric precision you shouldn't attempt at doing that transformation at all, since the "RGF93 v1 to WGS 84 (1)" null-transformation has an advertized accuracy of 1 metre.

Instead of using EPSG:4978, you'd probably want to use EPSG:4937 (ETRS89 geocentric) to avoid entirely the issues with WGS 84.

Le 16/10/2024 à 09:26, Christian DOULIAC a écrit :

Hello and sorry to spam this subject a bit.

I have created a workaround in my code by checking the "PJ_TYPE" and "proj_crs_promote_to_3D" the CRS if it is a "PJ_TYPE_GEOGRAPHIC_2D_CRS" or a "PJ_TYPE_PROJECTED_CRS".

It works pretty well except for one case :

Conversion from "EPSG:4978" to "EPSG:2154"
-In proj 9.2 without promotion the "proj_create_crs_to_crs_from_pj" function return me a transformation like this :
+proj=pipeline
+step +inv +proj=cart +ellps=WGS84
+step +proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000 +y_0=6600000 +ellps=GRS80
-In proj 9.2 with promotion it gave me this :
+proj=pipeline

+step +inv +proj=cart +ellps=GRS80

+step +proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000 +y_0=6600000 +ellps=GRS80

The only difference is the change from "+ellps=WGS84" to "+ellps=GRS80".

The 2nd transformation seems to be more accurate because "EPSG:2154" is based on the GRS80, is it true?
How can I recreate the past behavior with the C API? (Having "+ellps=WGS84"). Ideally I would like a code that is not specific for this case.

Thank you for taking the time to read me.

________________________________
De : PROJ <proj-bounces at lists.osgeo.org><mailto:proj-bounces at lists.osgeo.org> de la part de Christian DOULIAC via PROJ <proj at lists.osgeo.org><mailto:proj at lists.osgeo.org>
Envoyé : mercredi 2 octobre 2024 14:55:02
À : Even Rouault; proj at lists.osgeo.org<mailto:proj at lists.osgeo.org>
Objet : Re: [PROJ] [Questions] Vertical transformations advices for PROJ >= 9.1 version

Hello thanks for your quick answer,

Indeed adding --3d to cs2cs command give me the past behavior.
The commands that I wrote were only here for the example, to illustrate that I can recreate the behaviors.
My goal is to be more generic, I would like to be able to do CRS conversions without knownledges of them at first, using the PROJ C API.

I found some advices here https://github.com/OSGeo/PROJ/pull/3119#issuecomment-1069272183
and here https://github.com/OSGeo/PROJ/pull/3119#issuecomment-1069272183 .
What I understood from theses comments is that promoting CRS to 3D only make sense for PJ_TYPE_GEOGRAPHIC_2D_CRS and PJ_TYPE_PROJECTED_CRS types.
I would like to detect theses cases and promoting the CRS accordingly with the API.

My questions are multiple :
-first did I understood right?
-Are there others cases (or types) that I should think about?
-Is looking for the CRS type is enough or should I look for other information?
________________________________
De : Even Rouault <even.rouault at spatialys.com><mailto:even.rouault at spatialys.com>
Envoyé : mardi 1 octobre 2024 15:13:11
À : Christian DOULIAC; proj at lists.osgeo.org<mailto:proj at lists.osgeo.org>
Objet : Re: [PROJ] [Questions] Vertical transformations advices for PROJ >= 9.1 version

Le 01/10/2024 à 11:04, Christian DOULIAC via PROJ a écrit :

Hello, I would like some precision about CRS conversions.

I'm aware that there are some changes about vertical transformations with PROJ >= 9.1 as stated here : <https://github.com/OSGeo/PROJ/pull/3119> https://github.com/OSGeo/PROJ/pull/3119
Or here : https://proj.org/en/9.4/apps/cs2cs.html#cmdoption-cs2cs-3d "Starting with PROJ 9.1, both CRS need to be 3D for vertical transformation to possibly happen."

I tested it and I recreated the behavior of pre 9.1 and post 9.1 versions :

-PROJ 8.2.0 :
echo 0 0 1000 | cs2cs +proj=longlat +datum=WGS84 +geoidgrids=/proj-build/data/for_tests/egm96_15.tif +vunits=m +no_defs +type=crs +to +proj=longlat +datum=WGS84 +no_defs +type=crs -f "%.14f"
0.00000000000000 0.00000000000000 1017.16157913208008

-PROJ 9.2.0
echo 0 0 1000 | cs2cs +proj=longlat +datum=WGS84 +geoidgrids=/proj-build/data/for_tests/egm96_15.tif +vunits=m +no_defs +type=crs +to +proj=longlat +datum=WGS84 +no_defs +type=crs -f "%.14f"
0.00000000000000 0.00000000000000 1000.00000000000000

Add --3d to promote the target CRS to 3D.

Or use: cs2cs -f "%.14f" EPSG:4326+5773 EPSG:4979  (with lat, long order)

--
http://www.spatialys.com
My software is free, but my time generally not.

--
http://www.spatialys.com
My software is free, but my time generally not.
Butcher of all kinds of standards, open or closed formats. At the end, this is just about bytes.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20241016/a37eaa27/attachment.htm>