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

[Sandro Santilli]

On Wed, Jan 25, 2012 at 04:45:13PM +0100, Sandro Santilli wrote:

> The +towgs84 values look like swapped !

Alright, so with the NEW definitions for the known projections a full
roundtrip from LV03 to LV95 using the shift grid gives these distances
in meters:

        name      |     st_distance
  ----------------+---------------------
   Zimmerwald     | 0.00692816117873887
   Chrischona     | 0.00428319953182237
   Pfaender       | 0.00016855795287734
   La Givrine     | 0.00555304264270726
   Monte Generoso | 0.00804385131339353
  (5 rows)

For the record, here's the PostGIS query used to perform the round-trip:

  SELECT name, st_distance(lv95, 
   postgis_transform_geometry(      -- CH1903+      to   CH1903+/LV95
     postgis_transform_geometry(    -- CH1903       to   CH1903+ 
       postgis_transform_geometry(  -- CH1903/LV03  to   CH1903
          lv03,
          '+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 +no_defs  <>', -- 21781 (good)
          '+proj=longlat +ellps=bessel +towgs84=674.4,15.1,405.3,0,0,0,0 +no_defs  <>', -- 4149 (good)
          4149
       ),
       '+proj=longlat +ellps=WGS84 +nadgrids=CHENYX06a.gsb +no_defs', -- grid shift
       '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs', -- ch1903+ pretending to be WGS84
       4150),
     '+proj=longlat +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +no_defs  <>', -- 4150 (good)
     '+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=2600000 +y_0=1200000 +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs  <>', -- 2056 good
     2056))
  FROM euref;

The "lv95" and "lv03" field in the "euref" table contain official numbers
found in in the PDF downloaded from the swisstopo website.

--strk; 

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