[PROJ] Help needed for a cs2cs pipeline with shift grids
Even Rouault
even.rouault at spatialys.com
Tue Feb 20 10:05:55 PST 2024
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" <even.rouault at spatialys.com>
> *À: *"Nicolas Bellaiche" <nicolas.bellaiche at ign.fr>, "proj"
> <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 :
>
> 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
> 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.
>
--
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/20240220/d49b8983/attachment-0001.htm>
More information about the PROJ
mailing list