[Proj] Help with use of swisstopo NTv2 grid shift ?
Mikael Rittri
Mikael.Rittri at carmenta.com
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
More information about the Proj
mailing list