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

[Sandro Santilli]

On Wed, Jan 25, 2012 at 11:14:18AM +0000, Mikael Rittri wrote:
> > > 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.

Got it, thanks.

> > 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

I'm having your results with +init but not with the full-codes.

  ~$ cs2cs +init=epsg:21781 +to +init=epsg:4149 -f "%.9f"
  602030.680 191775.030
  7.466225970     46.878408135 0.000000000
  
  ~$ 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

Now, here's why !

  ~$ egrep '<(4149|21781)>' /usr/share/proj/epsg
  <4149> +proj=longlat +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +no_defs  <>
  <21781> +proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=600000 +y_0=200000 +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +un

I'm with 4.7.1 as well. On a 64bit.

  ~$ cs2cs --version
  Rel. 4.7.1, 23 September 2009

So, how could you get the same results ?
Did you manually update your <4149> defn with a more recent string ?
I see proj4 from SVN trunk has these other ones:

  ~$ egrep '<(4149|21781)>' /usr/src/proj/proj/proj/nad/epsg
  <4149> +proj=longlat +ellps=bessel +towgs84=674.4,15.1,405.3,0,0,0,0 +no_defs  <>
  <21781> +proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=600000 +y_0=200000 +ellps=bessel +towgs84=674.4,15.1,405.3,0,0,0,0 +units=m

The +towgs84 values look like swapped !

--strk;

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