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

molnar molnar at sas.elte.hu
Wed Feb 21 07:29:18 PST 2024 

@Nicolas

I get the same result as you, if the nadgrid and geoidgrid files are not 
selected.
Maybe, in your system, cs2cs does not know where to find grid files. So 
maybe setting this, or copying the grid files in the actual directory 
can help.
You can check it with:
  echo 0  56.000000 0  | ./cs2cs -f "%.3f"  +proj=eqc
  +R=6366197.72367581 +to_meter=100000 +lon_0=-2.3372291667E
  +nadgrids=@null +geoidgrids=@null +to +proj=lcc
  +lat_0=52 +lon_0=10 +lat_1=35 +lat_2=65 +x_0=4000000 +y_0=2800000
  +ellps=GRS80
That produces the same wrong result that you got, because no grid files 
are selected.
Please let me know, if it works! Thank you!

Gabor




On 2024-02-20 14:41, Nicolas Bellaiche wrote:
> @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