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

Thu Jan 26 00:37:27 PST 2012

> I'm having your results with +init but not with the full-codes.
> [...]
> 
> Now, here's why !
> [...]

I'm glad you found it. It is clear that we have been running different versions of Proj.4.

> So, how could you get the same results ?
> Did you manually update your <4149> defn with a more recent string ?

I did get 

    C:\Windows\system32>cs2cs --version
    Rel. 4.7.1, 23 September 2009

but when I look at my C:\OSGeo4W\share\proj\epsg, I have the
decimeter-resolution datum shift for CH1903:

    +towgs84=674.4,15.1,405.3,0,0,0,0

And this was introduced in Changeset 1824, from 28 Feb 2010, which is
after the Proj 4.7.1 release, so it seems I have been running a trunk
version without knowing it.  I downloaded OSGeo4W in November 2011,
without doing anything special to get the trunk version of Proj.4, 
but maybe that's what you get by default. 

Before that Changeset, the datum shift for CH1903 had millimeter
resolution, (identical to the shift for CH1903+).

By the way, the EPSG Remarks for the decimeter-resolution shift are

    "These parameters are derived from CH1903+ to ETRS89 (code 1647) and are
     used as from CH1903 as an approximation which is within the accuracy of
     the distortions in the CH1903 network."

As I understand it, if you want to use a 3- or 7-parameter datum shift for
CH1903, you can use the same as the one for CH1903+, since the datums are
approximately in the same place. But the local distortions of CH1903 mean
that the accuracy is only about 1.5 m, so someone must have figured that 
the millimeter resolution made no sense, and rounded to decimeters.

> So, how could you get the same results ?

Well, when you use cs2cs to convert between CH1903 / LV03 and
CH1903 / longlat, you are not really doing any datum shift, since it
is the same datum on both sides. So, the only important thing is to 
have the same datum shift on both sides; then they will cancel each
other and have no net effect. Even +towgs84=999,999,999 should work, if
you use it on both sides. I think you can omit the +towgs84 parameters
on both sides, too, in this case, but I dislike doing that: I can
never remember the semantics of doing so. 

So, when you wrote 

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

the +init on the left side gave you the millimeter-resolution datum shift,
since you are running Proj 4.7.1, but on the right side you have explicitly
used the decimeter-resolution datum shift. Of course, that gives a
difference of 2.6 cm in X, 4.4 cm in Y and 4.6 cm in Z; I am not sure
what that would be horizontally, but I think it explains your mismatches
up to 8 cm.
     Well, that's what you had already figured out. I guess the lesson is
that it can be risky to use +init for one CRS definition and an explicit
definition for the other CRS.

> Now next step would be making such transformation available with a single > proj call.
> 
> Would it be possible to create an init file containing all "pretending"
> referenced set of working transforms ?definitions so to be able to
> use <CH1903+> as the pivot and <CH1903>, <CH1903/LV03>, <CH1903+/LV96>
> and maybe <WGS84> to have a self-referenced set of working transforms ?
> 
> Is that what IGNF file is all about ?

I don't think I understand these questions. 

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 4:45 PM
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 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