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

Wed Oct 16 00:43:35 PDT 2024

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> de la part de Christian 
> DOULIAC via PROJ <proj at lists.osgeo.org>
> *Envoyé :* mercredi 2 octobre 2024 14:55:02
> *À :* Even Rouault; 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>
> *Envoyé :* mardi 1 octobre 2024 15:13:11
> *À :* Christian DOULIAC; 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 
>> <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/8fabab9f/attachment.htm>