[PROJ] Help needed for a cs2cs pipeline with shift grids

Nicolas Bellaiche nicolas.bellaiche at ign.fr
Tue Feb 20 06:41:32 PST 2024 

@Gabor

your pipeline gives me this result with Proj-9.3.1:
echo 0  56.000000 0  | ./cs2cs -f "%.3f"  +proj=eqc +R=6366197.72367581 +to_meter=100000 +lon_0=-2.3372291667E +nadgrids=fr_ign_ntf_r93.tif +geoidgrids=RAF09.gtx +to +proj=lcc +lat_0=52 +lon_0=10 +lat_1=35 +lat_2=65 +x_0=4000000 +y_0=2800000 +ellps=GRS80

3474840.045	2655360.932 0.000

@Even

I'm studying this, hopefully I find how it works without sweating too much!

Nicolas


----- Mail original -----
De: "molnar via PROJ" <proj at lists.osgeo.org>
À: "proj" <PROJ at lists.osgeo.org>
Envoyé: Lundi 19 Février 2024 21:15:54
Objet: Re: [PROJ] Help needed for a cs2cs pipeline with shift grids

Dear Nicolas,

If you really need a cs2cs string, that converts from grad to Lambert 
projection, here it is:

cs2cs -f "%.3f"  +proj=eqc +R=6366197.72367581 +to_meter=100000 
+lon_0=-2.3372291667E +nadgrids=fr_ign_ntf_r93.tif +geoidgrids=RAF09.gtx 
+to +proj=lcc +lat_0=52 +lon_0=10 +lat_1=35 +lat_2=65 +x_0=4000000 
+y_0=2800000 +ellps=GRS80


This produces exactly the same result as your cct code!
It simulates the grad to degree transformation with an equidistant 
cylindrical projection. (My favorite trick :-)) Input data (that is in 
grad) handled internally as meter, and the first part of the cs2cs 
command with equidistant cylindrical projection converts it to degree. 
As a surprise for me, nadgrids and geoidgrids should be applied before 
the +to switch.
Best wishes,
Gabor



On 2024-02-18 14:16, Even Rouault via PROJ wrote:
> Nicolas,
> 
> I believe the Helmert transformation is used when using IGNF codes,
> because the IGNF database has only recorded the use of the ntf_r93
> grid for NTF to RGF93, but not for NTF to ETRS89 (by the way, just
> reiterating that IGN is more than welcome to take over on maintenance
> & updates of the IGNF part of the PROJ database). And when
> transforming between 2 IGNF codes, PROJ by default will only consider
> transformations in the IGNF domain. But if you use EPSG codes, then
> ntf_r93 is available to transforme between NTF and ETRS89.
> 
> To get a vertical transformation, you also need to use 3D CRS.
> 
> For example
> 
> projinfo -s "NTF geographiques Paris (gr) + NGF-IGN69 height" -t
> IGNF:ETRS89LCC --3d --spatial-test intersects
> 
> reports the following pipeline:
> 
> unknown id, axis order change (2D) + NTF (Paris) to NTF (1) + NTF to
> RGF93 v1 (1) + Inverse of RGF93 v1 to NGF-IGN69 height (1) + RGF93 v1
> to ETRS89 (1) + axis order change (geographic3D horizontal) + ETRS89
> LAMBERT CONFORMAL CONIC, 1.6 m, France - mainland onshore., at least
> one grid missing
> 
> PROJ string:
> +proj=pipeline
>   +step +proj=unitconvert +xy_in=grad +xy_out=rad
>   +step +inv +proj=longlat +ellps=clrk80ign +pm=paris
>   +step +proj=push +v_3
>   +step +proj=cart +ellps=clrk80ign
>   +step +proj=xyzgridshift +grids=fr_ign_gr3df97a.tif
> +grid_ref=output_crs
>         +ellps=GRS80
>   +step +inv +proj=cart +ellps=GRS80
>   +step +proj=pop +v_3
>   +step +proj=vgridshift +grids=fr_ign_RAF18.tif +multiplier=1
>   +step +proj=lcc +lat_0=52 +lon_0=10 +lat_1=35 +lat_2=65 +x_0=4000000
>         +y_0=2800000 +ellps=GRS80
> 
> Here, as there's a mix of objects from different authorities (IGNF and
> EPSG), then EPSG transformations are considered, hence you get the
> gr3df97a grid.
> 
> $ echo 0  56.000000 0  | PROJ_NETWORK=ON bin/cs2cs -d 3  "NTF
> geographiques Paris (gr) + NGF-IGN69 height" IGNF:ETRS89LCC --3d
> 3474789.997    2655359.291 43.644
> 
> which is the application of the above pipeline:
> 
> $ echo 0  56.000000 0  | PROJ_NETWORK=ON bin/cct -d 3 +proj=pipeline
> +step +proj=unitconvert +xy_in=grad +xy_out=rad +step +inv
> +proj=longlat +ellps=clrk80ign +pm=paris +step +proj=push +v_3 +step
> +proj=cart +ellps=clrk80ign +step +proj=xyzgridshift
> +grids=fr_ign_gr3df97a.tif +grid_ref=output_crs +ellps=GRS80 +step
> +inv +proj=cart +ellps=GRS80 +step +proj=pop +v_3 +step
> +proj=vgridshift +grids=fr_ign_RAF18.tif +multiplier=1 +step +proj=lcc
> +lat_0=52 +lon_0=10 +lat_1=35 +lat_2=65 +x_0=4000000 +y_0=2800000
> +ellps=GRS80
> 
> If you want to get exactly the below pipeline, you'll have to tweak
> proj.db to create a compound CRS for "NTF geographiques Paris (gr) +
> NGF-IGN69 height", and then create a custom record in
> other_transformation for the transformation between this compound CRS
> and (a 3D projected CRS derived from) IGNF:ETRS89LCC. Cf
> data/sql/other_transformation_custom.sql for potential inspiration
> 
> Even
> 
> Le 15/02/2024 à 19:46, Nicolas Bellaiche via PROJ a écrit :
> 
>> Hi there,
>> 
>> I try to make a conversion using shiftgrids from IGNF: NTFP to
>> IGNF:ETRS89LCC. By default, it seems that it uses the helmert
>> approximation between the 2 crs and i cannot figure out how to use
>> the grid fr_ign_ntf_r93.
>> 
>> So far i've found how to create a pipeline that works with cct, but
>> i'd like to have a code that represents the crs and be able to go to
>> different destination crs with the cs2cs application. Any idea?
>> 
>> (helmert approx + no geoid with the standard codes)
>> 
>> echo 0 56  0 | ./cs2cs  + from IGNF:NTFP  +to IGNF:ETRS89LCC
>> 
>> result  : 3474788.21 2655358.23 0.00
>> 
>> (correct pipeline)
>> 
>> echo 0  56.000000 0  | PROJ_DEBUG=3  ./cct    +proj=pipeline +step
>> +proj=unitconvert +xy_in=grad +xy_out=rad +step +inv +proj=longlat
>> +ellps=clrk80ign +pm=paris  +step +proj=hgridshift
>> +grids=fr_ign_ntf_r93.tif +step +inv  +proj=vgridshift
>> +grids=RAF09.gtx  +step +proj=lcc +lat_0=52 +lon_0=10 +lat_1=35
>> +lat_2=65 +x_0=4000000 +y_0=2800000 +ellps=GRS80
>> 
>> result : 3474789.9972   2655359.2908       43.6421           inf
>> 
>> reference from THE official IGN/SGN dataset:
>> 
>> 3474789.997 2655359.291 43.642
>> 
>> Thanks a lot for your help,
>> 
>> Nicolas Bellaiche
>> 
>> _______________________________________________
>> 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.
> _______________________________________________
> PROJ mailing list
> PROJ at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/proj
_______________________________________________
PROJ mailing list
PROJ at lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/proj

More information about the PROJ mailing list