[Proj] Problems using grid shift files

[Mikael Rittri]

Dear Stefan,
This is not the answer you are looking for, but maybe it helps.

>From the source code in src/pj_datum_set.c, my impression
is that +nadgrids takes precedence over +towgs84, at least 
in Proj4 (perhaps not in OGRSpatialReference etc.).  

That would mean that in your definition of CRS 902056,
where you use both +towgs84=674.374,... and +nadgrids=@null, 
the +towgs84 would be ignored and +nadgrids=@null would be 
used, which is the same thing (as far as I understand) as 
saying that the longitude/latitude in CRS 902056 is the same 
as longitude/latitude in WGS84.  

But that's not your intention, is it?  
The datum of EPSG:2056 is CH1903+, and the datum of 
EPSG:21781 is CH1903 (without the +).  

I thought CH1903+ was some kind of improved version of CH1903,
with less distortion, although I don't know if your gridshift 
file is intended for CH1903 or for CH1903+.

Anyway, I don't the see the reason for using +nadgrids=@null
in your definition of CRS 902056.  I would have expected either
the +towgs84 parameter or the +nadgrids=chenyx06.gsb.

Best regards,

--
Mikael Rittri
Carmenta AB
Box 11354
SE-404 28 Göteborg
Visitors: Sankt Eriksgatan 5
SWEDEN
mikael.rittri at carmenta.com
www.carmenta.com 

________________________________

From: proj-bounces at lists.maptools.org [mailto:proj-bounces at lists.maptools.org] On Behalf Of Ziegler Stefan
Sent: den 19 januari 2009 17:00
To: proj
Subject: [Proj] Problems using grid shift files


Dear List

we have some problems using grid shift files. We are using a grid shift file which was developed by the Swiss Federal Office of Topography [1]. We just duplicated epsg 21781 (to 921781) and 2056 (to 902056) in the epsg file and added "+nadgrids=chenyx06.gsb" and "+nadgrids=@null":

21781: 
+proj=somerc +lat_0=46d57'08.66"N +lon_0=7d26'22.50"E +ellps=bessel +x_0=600000 +y_0=200000 +towgs84=674.374,15.056,405.346 +units=m +k_0=1 +no_defs 

2056:
+proj=somerc +lat_0=46d57'08.66"N +lon_0=7d26'22.50"E +ellps=bessel +x_0=2600000 +y_0=1200000 +towgs84=674.374,15.056,405.346 +units=m +k_0=1 +no_defs

921781: 
+proj=somerc +lat_0=46d57'08.66"N +lon_0=7d26'22.50"E +ellps=bessel +x_0=600000 +y_0=200000 +towgs84=674.374,15.056,405.346 +units=m +nadgrids=chenyx06.gsb +k_0=1 +no_defs 

902056:
+proj=somerc +lat_0=46d57'08.66"N +lon_0=7d26'22.50"E +ellps=bessel +x_0=2600000 +y_0=1200000 +towgs84=674.374,15.056,405.346 +units=m +nadgrids=@null +k_0=1 +no_defs

4326:
+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs

904326:
+proj=longlat +ellps=WGS84 +datum=WGS84 +nadgrids=@null

Here are the results:

21781 -> 2056: correct!
921781 -> 90256: correct (with grid shift file)!

but:
921781 -> 2056: wrong!
921781 -> 4326: wrong!
921781 -> 904326: wrong!

Is the grid shift only applied if the destination definiton has a "+nadgrids=@null"? But then why is the transformation to 904326 not correct? Has this something to do with this bug report [2]? Or I'm just mixing up some parameters? Same results with Proj 4.5 and 4.6.1. Any help is appreciated.

regards
Stefan

[1]: http://www.swisstopo.admin.ch/internet/swisstopo/en/home/topics/survey/lv03-lv95/chenyx06/distortion_grids.html
[2]: http://trac.osgeo.org/proj/ticket/22

Mit freundlichem Gruss
Stefan Ziegler
Leiter Aufsicht

Kanton Solothurn
Bau- und Justizdepartement
Amt für Geoinformation
Rötistrasse 4
4501 Solothurn
Telefon 032 627 75 96
Telefax 032 627 75 98
stefan.ziegler at bd.so.ch
http://www.so.ch