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

[Mikael Rittri]

> > 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
...
>
> Sorry about my ignorance: where does the command line above say you're
> converting to CH1903 (epsg:4149) rather than WGS84 ?

Right at the start, where it says proj instead of cs2cs.
The original proj utility cannot do any datum shift, only cs2cs can do that.
So when the proj command uses +init=epsg:21781, it can convert from longlat
to projected coordinates (the default behaviour), or vice versa (the -I for
inverse projection). 

> It could be I found my drift ? See this:
> 
> cs2cs +init=epsg:21781 +to +proj=longlat +ellps=bessel +towgs84=674.4,15.1,405.3,0,0,0,0 +no_defs -f "%.9f"
> 602030.680 191775.030
> 7.466225442     46.878408624 0.012045878

Strange, and I cannot reproduce it; I get identical results from proj and cs2cs:

>cs2cs +init=epsg:21781 +to +proj=longlat +ellps=bessel +towgs84=674.4,15.1,405.3,0,0,0,0 +no_defs -f "%.9f"
602030.680 191775.030
7.466225970     46.878408135 0.000000000

Or simpler:

>cs2cs +init=epsg:21781 +to +init=epsg:4149 -f "%.9f"
602030.680 191775.030
7.466225970     46.878408135 0.000000000

The 1.2 cm height difference you get, may be caused by your having different
+towgs84 parameters on your left and right hand side, but I don't understand
how that can be. (I am using OSGeo4W with Proj version 4.7.1, I think.)

> So do you confirm your multi-step would be as follows ?
> 
> (1) convert CH1903/LV03 to WGS84 using epsg:21781
> (2) apply NTv2 on the result of step 1 (get CH1903+)
> (3) convert CH1903+ back to CH1903+/LV95 using epsg:2056

Yes, more or less:

Step (1) is really conversion from CH1903/LV03 to CH1903 longlat.
(I think that's what you mean).

Then in step (2), we pretend that these longlat values are on the WGS84
ellipsoid, and we pretend that we use the NTv2 grid to get WGS84 longlat,
but we really get CH1903+ longlat. (That's how you already do it.)

And step (3) should be possible with the proj command or cs2cs.

Maybe it would be possible to combine step (1) and (2) to a single cs2cs
command, by replacing +proj=longlat by +proj=somerc etc. in the source CRS
that uses the NTv2 grid. But I am not sure what ellipsoid to specify in that
case. The +proj=somerc would need the Bessel ellipsoid to work properly, and
that would perhaps confuse the grid conversion. 

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: Wednesday, January 25, 2012 11:24 AM
To: Mikael Rittri
Cc: PROJ.4 and general Projections Discussions
Subject: Re: [Proj] Help with use of swisstopo NTv2 grid shift ?

On Wed, Jan 25, 2012 at 09:54:49AM +0000, Mikael Rittri wrote:

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

Good to know I'm on the right track !

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

Right. That's what I had.

> I have tried to reproduce your results, but my mismatches were ten times
> smaller, or at most 8 mm. 
...
> 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
...

Sorry about my ignorance: where does the commandline above say you're
converting to CH1903 (epsg:4149) rather than WGS84 ?

It could be I found my drift ? See this:

cs2cs +init=epsg:21781 +to +proj=longlat +ellps=bessel +towgs84=674.4,15.1,405.3,0,0,0,0 +no_defs -f "%.9f"
602030.680 191775.030
7.466225442     46.878408624 0.012045878

The destination projection above is epsg:4149

So do you confirm your multi-step would be as follows ?

 (1) convert CH1903/LV03 to WGS84 using epsg:21781
 (2) apply NTv2 on the result of step 1 (get CH1903+)
 (3) convert CH1903+ back to CH1903+/LV95 using epsg:2056

--strk;

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