1

I use the last versions of R -4.0.3- and sf -0.9.6- (thus relying on PROJ version 6.3.1).

I want to transform a layer from one system of coordinates (EPSG:27572) to another (EPSG:2154). This can be done with (accurate) or without (less accurate) using a transformation grid file. Once I have copied the appropriate transformation grid file (ntf_r93.gsb for PROJ 6) in the proj directory of sf (located at "C:\Program Files\R\R-4.0.3\library\sf\proj" on my PC), st_transform automatically use it and gives the most accurate results:

st_transform(MyLayer,crs=2154)

My question is:

  • how to get rid of the transformation grid, coming back to the non-gridded transformation, without removing the file in the proj directory?
  • Is there any sf command, or any parameter in st_transform to do this simply?
  • Perhaps using the pipeline parameter?
  • How?
2

You can check the pipeline that is used for the conversion with projinfo https://proj.org/apps/projinfo.html.

projinfo -s epsg:27572 -t epsg:2154

PROJ string:
+proj=pipeline +step +inv +proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=0 +k_0=0.99987742 +x_0=600000 +y_0=2200000 +ellps=clrk80ign +pm=paris +step +proj=hgridshift +grids=ntf_r93.gsb +step +proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000 +y_0=6600000 +ellps=GRS80

Test with cct https://proj.org/apps/cct.html to prove if your thoughts about the pipeline are correct. First conversion includes the gridshift step +step +proj=hgridshift +grids=ntf_r93.gsb, the second is without it.

echo 557994.61 2101992.70 0|cct +proj=pipeline +step +inv +proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=0 +k_0=0.99987742 +x_0=600000 +y_0=2200000 +ellps=clrk80ign +pm=paris +step +proj=hgridshift +grids=ntf_r93.gsb +step +proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000 +y_0=6600000 +ellps=GRS80
  606610.4430   6535955.8606        0.0000           inf

echo 557994.61 2101992.70 0|cct +proj=pipeline +step +inv +proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=0 +k_0=0.99987742 +x_0=600000 +y_0=2200000 +ellps=clrk80ign +pm=paris +step +proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000 +y_0=6600000 +ellps=GRS80
  606665.9540   6535960.4163        0.0000           inf

The results are different so editing the pipeline should resolve your problem. Unfortunately I do not know R and how to feed the pipeline for st_transform.

8
  • To feed the pipeline for st_transform in R, we just have to provide the PROJ string to the pipeline parameter: st_transform(MyLayer, pipeline="...") However, the change you made in the PROJ string does not give the good result (the one which is obtained without the proj grid file). Nov 9 '20 at 15:25
  • 1
    The correct PROJ string to provide is "+proj=pipeline +step +inv +proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=0 +k_0=0.99987742 +x_0=600000 +y_0=2200000 +ellps=clrk80ign +pm=paris +step +proj=push +v_3 +step +proj=cart +ellps=clrk80ign +step +proj=helmert +x=-168 +y=-60 +z=320 +step +inv +proj=cart +ellps=GRS80 +step +proj=pop +v_3 +step +proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000 +y_0=6600000 +ellps=GRS80", much longer probably because the grid transformation has to be replaced by the Helmert transformation, in the middle of the string. Nov 9 '20 at 15:46
  • I was wondering whether a simpler solution exists, just to prevent the use of the grid by st_transform. Nov 9 '20 at 15:46
  • By default, when no grid file exists in the sf library, st_transform "knows" the correct transformation to apply. It should be simple to mimic this case. Nov 9 '20 at 15:48
  • So what would happen if you edit the gridshift file +grids=ntf_r93.gsb -> +grids=foo.gsb? But my Proj version finds just one candidate operation projinfo -s epsg:27572 -t epsg:2154 --spatial-test intersects Candidate operations found: 1. I do not know how it could change the transformation path to use Helmert. My Proj version is rather old Rel. 6.3.2, May 1st, 2020.
    – user30184
    Nov 9 '20 at 15:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.