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

Wed Jan 25 04:54:49 EST 2012

Hello Sandro, 

I don't see that you do anything wrong with the 'trip' calculation. 

You say you get a mismatch of about 7e-07, which I assume is degrees,
so that's about 7 or 8 cm.  

I have tried to reproduce your results, but my mismatches were ten times
smaller, or at most 8 mm. 

One reason for a mismatch is that Switzerland, like many countries, has
its own official algorithm for datum shifts; the Swiss one is called
FINELTRA. (The reason is probably a conspiracy to create work opportunities
for GIS programmers.) The NTv2 file that the Swiss have published is
intended only as an approximation. The mismatches compared to FINELTRA can
be up to 10 or 20 cm, apparently, but is less than 2 cm in most of the
country. I got this info from 

http://www.swisstopo.admin.ch/internet/swisstopo/en/home/topics/survey/lv03-lv95/chenyx06/distortion_grids.html

and there is more info in the PDF notes, but I am not good at German. 

Anyway, here is how I made my computations:

First, convert the five test points from CH1903 / LV03 to CH1903 LongLat:

>proj +init=epsg:21781 -f "%.9f" -I

602030.680 191775.030
7.466225970     46.878408135
617306.300 268507.300
7.669595854     47.568440713
776668.105 265372.681
9.785678730     47.516696408
497313.292 145625.438
6.102781793     46.455347324
722758.810 87649.670
9.022387820     45.930487539

Note that I used -f "%.9f" to get higher precision in the output. 

Second, transform from CH1903 LongLat to CH1903+ LongLat, using the NTv2
file in the same way you did, but again using -f "%.9f":

>cs2cs +proj=longlat +ellps=WGS84 +nadgrids=CHENYX06a.gsb +no_defs +to +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs -f "%.9f"

7.466225970     46.878408135
7.466226678     46.878408104 0.000000000
7.669595854     47.568440713
7.669604076     47.568445851 0.000000000
9.785678730     47.516696408
9.785684998     47.516692402 0.000000000
6.102781793     46.455347324
6.102773347     46.455353520 0.000000000
9.022387820     45.930487539
9.022390666     45.930474254 0.000000000

Third, compare these positions with the official (FINELTRA-calculated)
CH1903+ LongLat values from Chapter 7. I did this on Ed Williams's
calculator http://williams.best.vwh.net/gccalc.htm , and got mismatches
of

7.0 mm
4.3 mm
0.1 mm
5.5 mm
8.1 mm

I don't understand how you could get ten times larger mismatches.
Maybe you have some rounding going on somewhere. 

However, there is something strange in your +somerc definition, where you
have

+x_0=6 +y_0=2

This should be 

+x_0=600000 +y_0=200000

for CH1903 / LV03 (EPSG:21781), and 

+x_0=2600000 +y_0=1200000

for CH1903+ / LV95 (EPSG:2056). But perhaps you just made a cut-and-paste
error in your text (otherwise the errors would have been hundreds of
kilometers).

Regards,

Mikael Rittri
Carmenta
Sweden
http://www.carmenta.com

-----Original Message-----
From: Sandro Santilli [mailto:sandro.santilli at gmail.com] On Behalf Of Sandro Santilli
Sent: Tuesday, January 24, 2012 6:12 PM
To: Mikael Rittri
Cc: PROJ.4 and general Projections Discussions
Subject: Re: [Proj] Help with use of swisstopo NTv2 grid shift ?

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