[Proj] Help with use of swisstopo NTv2 grid shift ?

[Sandro Santilli]

On Tue, Jan 24, 2012 at 12:32:40PM +0000, Mikael Rittri wrote:
> 
> Sandro Santilli wrote:
> 
> > It would surely help if I had an official set of points in CH1903
> > and CH1903+ to test my conversions against...
> 
> http://www.swisstopo.admin.ch/internet/swisstopo/en/home/topics/survey/sys/refsys/switzerland.parsysrelated1.37696.downloadList.97912.DownloadFile.tmp/swissprojectionen.pdf
> 
> There are five positions that are expressed in various Swiss systems.
> Unfortunately, CH1303 LongLat is not one of them, but LV03 is included,
> which I think is what EPSG calls "CH1903 / LV03" (EPSG:21781). So it is
> possible to find out the CH1903 LongLat coordinates via Proj.4, like this
> 
...
> 
> Which you can compare with the corresponding CH1903+ LongLat from Chapter 7:
> 
> 7 27 58.416328	46 52 42.269284 
> 7 40 10.574820	47 34  6.404965
> 9 47  8.465989	47 31  0.092644 
> 6  6  9.983811	46 27 19.272743 
> 9  1 20.606368	45 55 49.707052 
> 
> Hope this helps,

Thanks Mikael. It does partially help.

My early test with use of the gridshift give closer numbers
to the reference, but still a bit far.

These are distances from CH1903+ LongLat in Chapter 7 and:

  "proj": EPSG:21781 -> EPSG:4150 
  "trip": EPSG:21781 -> EPSG:4149 -> gridshift 

      name      |         proj         |         trip
----------------+----------------------+----------------------
 Zimmerwald     | 7.87246233515599e-07 | 7.60831787764113e-07
 Chrischona     | 9.71503452898624e-06 | 7.72176457507915e-07
 Pfaender       |  7.4382455747246e-06 | 7.18561592587409e-07
 La Givrine     | 1.05401600040938e-05 | 6.58608124095734e-07
 Monte Generoso | 1.36557256758873e-05 | 7.58098360996886e-07

Could EPSG:21781 -> EPSG:4149 have introduced that ~7.5e-07 units drift ?
Or am I doing it still wrong ...

The exact steps I took are as follow:

 (1) EPSG:21781 -> EPSG:4149 using standard EPSG entries:
     FROM (EPSG:21781): +proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +x_0=6 +y_0=2 +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs
     TO (EPSG:4149): +proj=longlat +ellps=bessel +towgs84=674.4,15.1,405.3,0,0,0,0 +no_defs

 (2) gridshift using custom entries:
     FROM: +proj=longlat +ellps=WGS84 +nadgrids=CHENYX06a.gsb +no_defs
     TO (EPSG:4150 pretending to be EPSG:4326): +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs

Am I doing step (2) correctly ? It's basically pretending the input is
just a gridshift away from the target (which is really EPSG:4150 in disguise)
   

--strk; 

  ()   Free GIS & Flash consultant/developer
  /\   http://strk.keybit.net/services.html