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ées?lang=fr
<https://www.sandre.eaufrance.fr/jeu-de-donnees/projection-des-coordonnées?lang=fr>)
echo 056.000000 0| ./cs2cs +from EPSG:4807 +to EPSG:3034
But the result doesn't make any sense:
-1640975.3910099142.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
This works well:
echo 056.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?
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
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?
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.
There is no simple code like EPSG or IGNF related to it? Can we create
one easily?
yes, cf
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)
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?
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.
Nicolas Bellaiche
------------------------------------------------------------------------
*De: *"Even Rouault" <[email protected]>
*À: *"Nicolas Bellaiche" <[email protected]>, "proj"
<[email protected]>
*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 :
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 560 | ./cs2cs+ from IGNF:NTFP+to IGNF:ETRS89LCC
result : 3474788.212655358.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.9972655359.29143.642
Thanks a lot for your help,
Nicolas Bellaiche
_______________________________________________
PROJ mailing list
[email protected]
https://lists.osgeo.org/mailman/listinfo/proj
--
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.
_______________________________________________
PROJ mailing list
[email protected]
https://lists.osgeo.org/mailman/listinfo/proj
