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

[Paul Kelly]

Hello Brad

On Mon, 4 Sep 2006, Brad Douglas wrote:

> On Sun, 2006-08-27 at 00:00 +0100, Paul Kelly wrote:
>> 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).
>
> This is definitely desirable.  I would almost consider not being able to
> project Z data a bug.

Do you mean it was desirable the way the default behaviour was or the way 
it is now? If the former, have you any examples? Just out of interest. I 
really couldn't think of any but I am neither a geographer nor a 
geodisist...

>> 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 would reject projections involving different geoids.  Otherwise, Z
> data may be corrupted.  At a minimum, a big fat disclaimer, but I would
> prefer checking geoids and error if needed.

But the problem is there is no way of telling what the z data is relative 
to; it might even mean nothing height-related at all. At present there is 
nowhere that this information is stored. I suggested in an earlier e-mail 
that perhaps it could be stored as part of the metadata with each vector 
map, but it is a big job I think.

>> 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.
>
> Does your patch lacks output of Z data (I see it does calculate values)?

No; the z data passes straight through untouched. pj_do_transform() 
modifies only the x and y co-ordinates (it has pointers to these passed to 
it).

> Other than the comments above, I'd like to see it applied.  Great job!

Already done a little while ago :)

Paul