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

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

(attachments)

v.proj.diff (1.35 KB)

I agree with your proposal.

On 8/26/06, Paul Kelly <paul-grass@stjohnspoint.co.uk> 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).

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


grass-dev mailing list
grass-dev@grass.itc.it
http://grass.itc.it/mailman/listinfo/grass-dev

ā€“
David Finlayson

David Finlayson wrote:

I agree with your proposal.

+1

Maciek

Hi,

me too... would be committed to CVS?

Martin

2006/8/27, David Finlayson <david.p.finlayson@gmail.com>:

I agree with your proposal.

On 8/26/06, Paul Kelly <paul-grass@stjohnspoint.co.uk> 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).

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

_______________________________________________
grass-dev mailing list
grass-dev@grass.itc.it
http://grass.itc.it/mailman/listinfo/grass-dev

--
David Finlayson
_______________________________________________
grass-dev mailing list
grass-dev@grass.itc.it
http://grass.itc.it/mailman/listinfo/grass-dev

--
Martin Landa <landa.martin@gmail.com> * http://gama.fsv.cvut.cz/~landa *

OK I have committed the change now.

In the longer term I think this is an area we could start to think about vertical datum support for. Helena has mentioned it on the list before but I could never see an easy way to associate a particular property of a raster or vector map with it being elevation data. But it could certainly be done on a per-map basis for 3-D vectors.
Just a vague thought for the future.

Paul

On Wed, 30 Aug 2006, Martin Landa wrote:

Hi,

me too... would be committed to CVS?

Martin

2006/8/27, David Finlayson <david.p.finlayson@gmail.com>:

I agree with your proposal.

On 8/26/06, Paul Kelly <paul-grass@stjohnspoint.co.uk> 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).

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

_______________________________________________
grass-dev mailing list
grass-dev@grass.itc.it
http://grass.itc.it/mailman/listinfo/grass-dev

--
David Finlayson
_______________________________________________
grass-dev mailing list
grass-dev@grass.itc.it
http://grass.itc.it/mailman/listinfo/grass-dev

--
Martin Landa <landa.martin@gmail.com> * http://gama.fsv.cvut.cz/~landa *

Hello Paul,

do you want me to merge that to the 6.2-release branch?

Markus

On Wed, Aug 30, 2006 at 04:27:38PM +0100, Paul Kelly wrote:

OK I have committed the change now.

In the longer term I think this is an area we could start to think about
vertical datum support for. Helena has mentioned it on the list before but
I could never see an easy way to associate a particular property of a
raster or vector map with it being elevation data. But it could certainly
be done on a per-map basis for 3-D vectors.
Just a vague thought for the future.

Paul

On Wed, 30 Aug 2006, Martin Landa wrote:

>Hi,
>
>me too... would be committed to CVS?
>
>Martin
>
>2006/8/27, David Finlayson <david.p.finlayson@gmail.com>:
>>I agree with your proposal.
>>
>>
>>On 8/26/06, Paul Kelly <paul-grass@stjohnspoint.co.uk> 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).
>>
>>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
>>
>>_______________________________________________
>>grass-dev mailing list
>>grass-dev@grass.itc.it
>>http://grass.itc.it/mailman/listinfo/grass-dev
>>
>>
>>
>>
>>
>>--
>>David Finlayson
>>_______________________________________________
>>grass-dev mailing list
>>grass-dev@grass.itc.it
>>http://grass.itc.it/mailman/listinfo/grass-dev
>>
>>
>
>
>--
>Martin Landa <landa.martin@gmail.com> * http://gama.fsv.cvut.cz/~landa *
>

_______________________________________________
grass-dev mailing list
grass-dev@grass.itc.it
http://grass.itc.it/mailman/listinfo/grass-dev

--
Markus Neteler <neteler itc it> http://mpa.itc.it/markus/
ITC-irst - Centro per la Ricerca Scientifica e Tecnologica
MPBA - Predictive Models for Biol. & Environ. Data Analysis
Via Sommarive, 18 - 38050 Povo (Trento), Italy

Ha ha I was hoping I wouldn't have to make that decision! Well it does change functionality (in certain cases), but 6.2 is a major change in version number so perhaps that is OK. And anyway it is a bugfix. So I would say yes.

Paul

On Wed, 30 Aug 2006, Markus Neteler wrote:

Hello Paul,

do you want me to merge that to the 6.2-release branch?

Markus

On Wed, Aug 30, 2006 at 04:27:38PM +0100, Paul Kelly wrote:

OK I have committed the change now.

In the longer term I think this is an area we could start to think about
vertical datum support for. Helena has mentioned it on the list before but
I could never see an easy way to associate a particular property of a
raster or vector map with it being elevation data. But it could certainly
be done on a per-map basis for 3-D vectors.
Just a vague thought for the future.

Paul

On Wed, 30 Aug 2006, Martin Landa wrote:

Hi,

me too... would be committed to CVS?

Martin

2006/8/27, David Finlayson <david.p.finlayson@gmail.com>:

I agree with your proposal.

On 8/26/06, Paul Kelly <paul-grass@stjohnspoint.co.uk> 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).

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

_______________________________________________
grass-dev mailing list
grass-dev@grass.itc.it
http://grass.itc.it/mailman/listinfo/grass-dev

--
David Finlayson
_______________________________________________
grass-dev mailing list
grass-dev@grass.itc.it
http://grass.itc.it/mailman/listinfo/grass-dev

--
Martin Landa <landa.martin@gmail.com> * http://gama.fsv.cvut.cz/~landa *

_______________________________________________
grass-dev mailing list
grass-dev@grass.itc.it
http://grass.itc.it/mailman/listinfo/grass-dev

--
Markus Neteler <neteler itc it> http://mpa.itc.it/markus/
ITC-irst - Centro per la Ricerca Scientifica e Tecnologica
MPBA - Predictive Models for Biol. & Environ. Data Analysis
Via Sommarive, 18 - 38050 Povo (Trento), Italy

_______________________________________________
grass-dev mailing list
grass-dev@grass.itc.it
http://grass.itc.it/mailman/listinfo/grass-dev

Paul,

I tried reprojecting a 3D points vector (3 arcsec strm) with v.proj,
from latlong WGS84 to UTM 34, and I see *no difference* in either x,y
or z coordinate, no matter whether the new implemented -z flag is used
or not. Is this OK?

Maciek

Did you change datum?

On 9/3/06, Maciej Sieczka <tutey@o2.pl> wrote:

Paul,

I tried reprojecting a 3D points vector (3 arcsec strm) with v.proj,
from latlong WGS84 to UTM 34, and I see no difference in either x,y
or z coordinate, no matter whether the new implemented -z flag is used
or not. Is this OK?

Maciek


grass-dev mailing list
grass-dev@grass.itc.it
http://grass.itc.it/mailman/listinfo/grass-dev

ā€“
David Finlayson

David Finlayson wrote:

Did you change datum?

Dylan,

No. I reprojected from latlon on WGS84 to UTM zone 34, which is also
based on WGS84 datum.

Maciek

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.

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.

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)?

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

--
Brad Douglas <rez touchofmadness com> KB8UYR
Address: 37.493,-121.924 / WGS84 National Map Corps #TNMC-3785

Hello Maciek
That is very strange. Can you post the output from v.proj please to see if there are any clues?

Paul

On Sun, 3 Sep 2006, Maciej Sieczka wrote:

Paul,

I tried reprojecting a 3D points vector (3 arcsec strm) with v.proj,
from latlong WGS84 to UTM 34, and I see *no difference* in either x,y
or z coordinate, no matter whether the new implemented -z flag is used
or not. Is this OK?

Maciek

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 :slight_smile:

Paul

Paul Kelly wrote:

That is very strange. Can you post the output from v.proj please to see
if there are any clues?

Hi Paul!

Attached are 2 sample Grass locations.

zak_ll is latlong on WGS84 and contains an input 'test' 3D points
vector (created with 'v.random -z output=test n=10 zmin=0.0 zmax=3000.0').

zak_utm34 is UTM34 on WGS84 and contains 2 v.proj outputs:
'test_reproj_Z_no' and 'test_reproj_Z_yes', created by the following
commands, respectively:

$ v.proj input=test location=zak_ll mapset=test output=test_reproj_Z_no

$ v.proj -z input=test location=zak_ll mapset=test output=test_reproj_Z_yes

(You'll find a full v.proj output in the attached text file.)

For both output files I have uploaded their x,y,z coordinates into a
table, and I find there is *no difference* in either coordinate,
although the 1st vector was created using the v.proj with -z switch
while the 2nd without this switch. Looks like '-z' doesn't work?

Maciek

(attachments)

vproj_output.txt (1.5 KB)
zak_ll.tar.bz2 (2.11 KB)
zak_utm34.tar.bz2 (2.88 KB)

Maciej Sieczka wrote:

Paul Kelly wrote:

That is very strange. Can you post the output from v.proj please to see
if there are any clues?

Hi Paul!

Attached are 2 sample Grass locations.

snip

Hi Paul!

Where do we stand with this?

Maciek

On Tue, 5 Dec 2006, Maciej Sieczka wrote:

Paul Kelly wrote:

Hello Maciek
We're talking about your example at
http://grass.itc.it/pipermail/grass5/2006-September/025668.html, right?

Right.

In that case both your co-ordinate systems are based on the same
ellipsoid (WGS84) so it is normal and expected for the ellipsoidal
height to stay the same.

Sorry, but I'm quite a dummy as to ellipsoids vs height. So I should
expect a difference whether using the '-z' switch only when
reprojecting between 2 different ellipsoids, correct?

Yes. It is really not a very useful feature at all, quite confusing actually, and that's why I thought it was important to make it optional. Ellipsoidal height is not very useful because you can't use it to work out e.g. if a ball is going to roll down a hill. In theory the ellipsoidal height of the bottom of a hill could be higher than the top of the hill - what is more important in a practical situation is the height above the geoid, which is what most vertical datums use. And that is totally separate from the 2-D re-projection that the [rv].proj modules do. We've talked about better vertical datum support in GRASS but it is so complicated and would just open up a whole can of worms that I don't want to get into right now...

I'm not clear on what output you were expecting?

Well, you wrote yourself "That is very strange. Can you post the output
from v.proj please to see if there are any clues?". So I did, and
waited for your response.

Yes well I didn't quite pick up on that then either (the ellipsoids being the same). It is quite hard to get your head around it and understand it and I only saw that was the problem when looking at it afresh!

Paul

Paul,

Summarising:

Could you please confirm/deny the 2 statemements below - so that I'm
sure I'm not screwing/missunderstanding something:

1. The -z switch is not the default; but it *was* - before you made it
and option few months ago.

2. One should use it only when his data is 3D, he is sure the Z coord
is an an ellipsoidal height, and he is transforming between 2 different
ellipsoids.

3. Even if the Z coord of my vector is an ellipsoidal height, using -z
is not supposed to make any change when transformimg between eg.
UTM/WGS84 and ll/WGS84; *it doesn't do any harm though*.

I know I'm slow on this, but please bear with me :). I don't want to
mess something up in future with my data or advise something stupid to
others.

Cheers,
Maciek

On Tue, 5 Dec 2006, Maciej Sieczka wrote:

Paul,

Summarising:

Could you please confirm/deny the 2 statemements below - so that I'm
sure I'm not screwing/missunderstanding something:

1. The -z switch is not the default; but it *was* - before you made it
and option few months ago.

Yes. There used to be no way to avoid that effect. Now you have to explicitly enable it.

2. One should use it only when his data is 3D, he is sure the Z coord
is an an ellipsoidal height, and he is transforming between 2 different
ellipsoids.

No - even then you should not use it. In fact I can't think of any reason why anybody would want to use it under normal circumstances, and IMHO it was a very serious bug that it was the default before. GPS receivers often output WGS84 ellipsoidal height. In general you want to keep the height in that datum rather than transfer it to another ellipsoid. Generally different ellipsoids are used only for 2-D projection and the height above the ellipsoid is not considered - a separate vertical datum, often based on sea level or something like that is used for height measurements.

An example: in Northern Ireland the Ordnance Survey NI (national mapping agency) publishes a set of gridded values that represent the average height differential between WGS84 ellipsoidal height and Mean Sea Level Belfast (the common local vertical datum) at each point. This can thus be used to convert GPS measurements to the local geoid-based MSL Belfast vertical datum. Ellipsoidal height above the Modified Airy ellipsoid (used for Irish Grid 2-D positions) is not relevant at all.

3. Even if the Z coord of my vector is an ellipsoidal height, using -z
is not supposed to make any change when transformimg between eg.
UTM/WGS84 and ll/WGS84; *it doesn't do any harm though*.

No it shouldn't do any harm, although this is a matter/question for PROJ.4 and not for GRASS - GRASS just uses the PROJ pj_transform() function and the transformation happens in there.

I know I'm slow on this, but please bear with me :). I don't want to
mess something up in future with my data or advise something stupid to
others.

I understand - this is an important issue and does not seem to be discussed very widely so it is good that we are doing so here.

Paul

Paul Kelly wrote:

this is an important issue and does not seem to be
discussed very widely so it is good that we are doing so here.

Cheers Paul!

Continuing this - can you please look at my post on a related issue
from several months ago and let me know if I'm writing something wrong
there?:

http://www.nabble.com/-GRASSLIST%3A10786--Re%3A-SRTM-data-transformation-tf2068135.html#a5693757

Thanks!

Maciek