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

Tue Feb 27 10:32:44 PST 2024

Just a follow up, 

I've added a sql file to the sql files that build the proj.db sqlite database with these 2 lines: 

-- NTFP+NGF-IGN69 height 

INSERT INTO compound_crs VALUES('RDI','NTFP.IGN69','NTF geographiques Paris (gr) + NGF-IGN69 height','','IGNF','NTFP','EPSG','5720',0); 

INSERT INTO usage VALUES('RDI','USAGE_COMPOUND_CRS_NTFP.IGN69','compound_crs','RDI','NTFP.IGN69','PROJ','EXTENT_UNKNOWN','PROJ','SCOPE_UNKNOWN'); 

and recompiled the library. It works perfectly. 

echo 0 56.000000 0 | ./cs2cs -f "%.3f" --3d +from IGNF:NTFP +to IGNF:ETRS89LCC --> 3474788.212 2655358.229 42.989 

echo 0 56.000000 0 | ./cs2cs -f "%.3f" --3d +from RDI:NTFP.IGN69 +to IGNF:ETRS89LC --> 3474789.997 2655359.291 43.644 

I will see if someone in the Institute could clarify these codes and update them in a future release of PROJ. 

Thanks for your help, 

N. 

De: "Even Rouault" <even.rouault at spatialys.com> 
À: "Nicolas Bellaiche" <nicolas.bellaiche at ign.fr> 
Cc: "proj" <PROJ at lists.osgeo.org> 
Envoyé: Mardi 20 Février 2024 19:05:55 
Objet: Re: [PROJ] Help needed for a cs2cs pipeline with shift grids 

Le 20/02/2024 à 18:53, Nicolas Bellaiche a écrit : 

Hi Even 

Thanks for your kind and developed answer. 
I've tried with the EPSG equivalents: ( [ https://www.sandre.eaufrance.fr/jeu-de-donnees/projection-des-coordonn%C3%A9es?lang=fr | https://www.sandre.eaufrance.fr/jeu-de-donnees/projection-des-coordonnées?lang=fr ] ) 

echo 0 56.000000 0 | ./cs2cs +from EPSG:4807 +to EPSG:3034 

But the result doesn't make any sense: 

-1640975.39 10099142.17 0.00 

Is there something with the units that is not correct maybe? 

No, just the usual axis order issue as EPSG geographic CRS are lat, long ordered, whereas IGNF uses the opposition order. So: 

$ echo 56 0 | cs2cs EPSG:4807 EPSG:3034 
2655358.23 3474788.21 0.00 

BQ_BEGIN

This works well: 

echo 0 56.000000 0 | PROJ_NETWORK=ON ./cs2cs -d 3 "NTF geographiques Paris (gr) + NGF-IGN69 height" IGNF:ETRS89LCC --3d 

and I have a few questions about it: 

1) the result remains unchanged if I remove or put PROJ_NETWORK to OFF. What does it do exactly? 

BQ_END
If you have already all the grids in the PROJ_DATA directory, PROJ_NETWORK will not do anything. This is just for people who don't have grids installed locally 

BQ_BEGIN

2) Does the string " NTF geographiques Paris (gr) + NGF-IGN69 height" correspond to an entry in the CRS database proj.db or is it interpreted as a compounded system because PROJ parses the string? 

BQ_END

Yes, there's logic in PROJ to split "A + B" into A and B, do simple CRS lookups for A and B in the database, and create a compound CRS from that. 
BQ_BEGIN

There is no simple code like EPSG or IGNF related to it? Can we create one easily? 

BQ_END

yes, cf [ https://github.com/OSGeo/PROJ/blob/master/data/sql/transformations_czechia_extra.sql#L86 | https://github.com/OSGeo/PROJ/blob/master/data/sql/transformations_czechia_extra.sql#L86 ] for a custom CompoundCRS 

You can also get the potential SQL statements (that can be simplified and improved) with: 

$ bin/projinfo IGNF:NTFP+EPSG:5720 -o SQL --output-id SOME_AUTH:MY_COMPOUND_CRS -q 
INSERT INTO geodetic_crs VALUES('SOME_AUTH','COMPONENT_MY_COMPOUND_CRS_1','NTF geographiques Paris (gr)','','geographic 2D','EPSG','6425','EPSG','6807',NULL,0); 
INSERT INTO scope VALUES('SOME_AUTH','SCOPE_geodetic_crs_COMPONENT_MY_COMPOUND_CRS_1','NATIONALE, HISTORIQUE',0); 
INSERT INTO extent VALUES('SOME_AUTH','EXTENT_geodetic_crs_COMPONENT_MY_COMPOUND_CRS_1','FRANCE METROPOLITAINE (CORSE COMPRISE)','FRANCE METROPOLITAINE (CORSE COMPRISE)',41,52,-5.5,10,0); 
INSERT INTO usage VALUES('SOME_AUTH','USAGE_GEODETIC_CRS_COMPONENT_MY_COMPOUND_CRS_1','geodetic_crs','SOME_AUTH','COMPONENT_MY_COMPOUND_CRS_1','SOME_AUTH','EXTENT_geodetic_crs_COMPONENT_MY_COMPOUND_CRS_1','SOME_AUTH','SCOPE_geodetic_crs_COMPONENT_MY_COMPOUND_CRS_1'); 
INSERT INTO compound_crs VALUES('SOME_AUTH','MY_COMPOUND_CRS','NTF geographiques Paris (gr) + NGF-IGN69 height','','SOME_AUTH','COMPONENT_MY_COMPOUND_CRS_1','EPSG','5720',0); 
INSERT INTO usage VALUES('SOME_AUTH','USAGE_COMPOUND_CRS_MY_COMPOUND_CRS','compound_crs','SOME_AUTH','MY_COMPOUND_CRS','PROJ','EXTENT_UNKNOWN','PROJ','SCOPE_UNKNOWN'); 

This could be simplified by just taking the last 2 statements and modify 'SOME_AUTH','COMPONENT_MY_COMPOUND_CRS_1' to 'IGNF', 'NTFP' (I believe this SQL synthesis code only references EPSG objects, hence the redefinition of an equivalent of IGNF:NTFP) 

BQ_BEGIN

3) Where can I see where this crs"NTF geographiques Paris (gr) + NGF-IGN69 height" is defined in the sql commands used to build it? 

BQ_END
well, if you want all the lineage, this is a cascade of INSERT statement. You have to fetch the geographic and vertical CRS, their datum, coordinate systems, axis, units in the relevant tables. 

BQ_BEGIN

Nicolas Bellaiche 

De: "Even Rouault" [ mailto:even.rouault at spatialys.com | <even.rouault at spatialys.com> ] 
À: "Nicolas Bellaiche" [ mailto:nicolas.bellaiche at ign.fr | <nicolas.bellaiche at ign.fr> ] , "proj" [ mailto:PROJ at lists.osgeo.org | <PROJ at lists.osgeo.org> ] 
Envoyé: Dimanche 18 Février 2024 15:16:44 
Objet: Re: [PROJ] Help needed for a cs2cs pipeline with shift grids 

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 : 

BQ_BEGIN

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 [ mailto:PROJ at lists.osgeo.org | PROJ at lists.osgeo.org ] [ https://lists.osgeo.org/mailman/listinfo/proj | https://lists.osgeo.org/mailman/listinfo/proj ] 

BQ_END

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

BQ_END

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

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20240227/dd388efb/attachment-0001.htm>