[GRASS-dev] v.proj transforming z co-ordinates

[Paul Kelly]

I suspected v.proj was doing this for years (since the introduction of the 
new vector format) but never got round to verifying the behaviour until 
yesterday. When re-projecting a 3-D vector file, if datum transformation 
parameters are specified using towgs84= notation (rather than nadgrids=), 
v.proj will assume the z co-ordinates are ellipsoidal heights and 
transform them accordingly.

I think there are probably very few circumstances in which this behaviour 
would be desirable. It is much more common for height data to be specified 
relative to a vertical datum, e.g. mean sea level, than relative to the 
ellipsoid (which is normally only used for horizontal distances).

GPS data is more likely to include ellipsoidal height, but again it is 
(IMHO) of questionable usefulness to transform it to another ellipsoid: 
the geoid models that exist to transform GPS data to the local vertical 
datum in regular use (OSGM02 for Great Britain and Northern Ireland is the 
one I'm familiar with) convert directly from WGS84 ellipsoidal height data 
to height above mean sea level.

I propose the attached patch to add a specific flag to v.proj to enable 
the current behaviour, otherwise the z co-ordinates are always left 
untouched. I thought it was worth explaining and discussing properly 
because it changes the current behaviour and will give different results 
in certain circumstances, but IMHO the current behaviour is wrong.

Any comments?

Paul
-------------- next part --------------
Index: main.c
===================================================================
RCS file: /grassrepository/grass6/vector/v.proj/main.c,v
retrieving revision 1.10
diff -u -r1.10 main.c
--- main.c	19 Jul 2006 02:51:38 -0000	1.10
+++ main.c	26 Aug 2006 22:58:13 -0000
@@ -33,6 +33,7 @@
     struct Map_info Out_Map;
     struct { 
 	struct Flag *list;               /* list files in source location */
+	struct Flag *transformz;
     } flag;
 
     G_gisinit (argv[0]);
@@ -70,6 +71,10 @@
     flag.list->key = 'l';
     flag.list->description = _("List vector files in input location and exit (a dummy value must be given for input)");
 
+    flag.transformz = G_define_flag();
+    flag.transformz->key = 'z';
+    flag.transformz->description = _("(3-D vectors only) Assume z co-ordinate is ellipsoidal height and transform appropriately");
+   
     if (G_parser (argc, argv)) exit (EXIT_FAILURE);
 		 
     /* start checking options and flags */
@@ -192,7 +197,8 @@
 
 	if (type == -1) G_fatal_error(_("Reading input dig file.")) ;
 	if ( type == -2) break;
-	if(pj_do_transform( Points->n_points, Points->x, Points->y, Points->z,
+	if(pj_do_transform( Points->n_points, Points->x, Points->y, 
+			    flag.transformz->answer? Points->z : NULL,
 		              &info_in,&info_out)<0) 
 	{ 
 	    fprintf(stderr, _("Error in pj_do_transform\n"));