geometric2geopotential

[Huw Lewis]

There was some confusion about the implementation of the geodesy routine geometric2geopotential to convert between height scales expressed relative to the geoid when the routines are based on Somigliana's equation, relative to the WGS-84 ellipsoid.

Although the same routine can be used to do height conversions relative to the ellipsoid or geoid to good accuracy, it is necessary to clarify the documentation to avoid confusion for users.

Email discussions:

Kent: 

We understand that ROPP doesn't include any geoid calculation
and that one can use the above routine to calculate the 
geopotential height to a good accuracy. However, I think it is 
confusing to say that the routine does calculate the 
geopotential height - since it doesn't. One could try to say 
that it calculates the 'geopotential height wrt the WGS84' but 
this is still not correct since the WGS84 is not an equi-potential
surface for the geopotential.

In addition, be clear in the documentation (+source code) about 
what the function does and also therein note that it can in fact 
be used to obtain the geopotential height but in doing so one 
makes a small error; and then give the equations that back up this 
(and estimate how big the error is).

Dave:

Strictly, we know it should be 'geopotential (as a height) above
geoid' (nominally EGM96), but as I understand it, and for reasons
unexplained, the coefficients datum is the ellipsoid (WGS84). Since the
undulation is included in the calculation, the end result is a geoid
datum, but we are aware of the small error that arises because the
undulation value strictly ought to be in geopotential too. I recall
Chris working out the error in geopot.ht to be in the order of a few mm
at 10km which is 2-3 orders of magnitude less than the accuracy of the
geoid (undulation) itself (EGM-96 claims an accuracy ~1m).

Hans:

I think it is OK to treat the geopotential heights the way
it is currently done in ROPP. However, it should be stated
very clearly in the header of the Fortran source file that
the geopotential heights are above the ellipsoid and not
above the geoid. Now the header only says that the function
returns 'geopotential heights' which per (standard) definition
is above the geoid. Most users will probably get this wrong
if they only read the header and not the formulas actually
implemented.

Dave:

My memory was faulty when I said that the undulation is being applied;
this is not the case *within* the routine, so if used stand-alone,
undulation would need to be applied by the user. This isn't at all
obvious from the code comments. In use, whatever datum is used for the
input 'h' is preserved for the output of the function call.

Dave:

Recall that the BUFR specification, (pg.17) and ROPP netCDF (pg.58)
geometric heights are supposed to be wrt geoid, NOT ellipsoid. From
Mikes' previous (indirect) analysis, direct questioning and recent
experience (undulation bug!) we know the impact of not applying
undulation (or applying it twice), so we are now sure that all current
GPSRO data is provided as 'geometric heights wrt geoid' and meets the
specification. The UCAR atmPrf netCDF files contain documented 'mean sea
level altitudes' so can be assumed consistent too (e.g. the output from
ucar2ropp is consistent with the UCAR-encoded BUFR).

So converting these geometric/geoid heights to geopotential/geoid
heights with geometric2geopotential() preserves the [EGM96] geoid datum
- and hence the ROPP output is consistent with our definition.

The routine itself, in isolation, assumes the input (& output) to be
geometric/ellipsoid (since the gravity & effective radius use
coefficients on WGS84) - hence the header comments. But within ROPP, 'h'
is actually given wrt geoid. It's true that the routine isn't actually
going from geometric/ellipsoid to geopotential/geoid, since the
undulation has to be applied externally. The end-to-end transform would
need to be either:

1) geometric/ellipsoid (-undulation) -> geometric/geoid -> [routine] 
-> geopotential/geoid

or

2) geometric/ellipsoid -> [routine] -> geometric/geoid  (-undulation) 
-> geopotential/geoid

Ie include undulation before or after the geometric2geopotential()
routine. (1) is the better way (as undulation is a geometric/ellipsoid
quantity - and in practice this is exactly how ROPP works as the
undulation has been applied at source - though as previously noted, this
does result in a very small error due to the datum shift, not present in
(2) - which has it's own error unless the undulation is also pre-
converted to geopotential wrt geoid. I would expect that for that either
errors are both insignificant so it doesn't matter which way round it's
done. (Something TBC).

I agree we should clarify the header comments regarding the datum,
but I don't think it's necessary to change the name of the routine
('geopotential' can be any datum you like, as long is it's specified) .
We could consider modifying the calling interface, e.g. to have an
optional 'undulation' argument, which if present is subtracted from the
altitude before the geopotential transform. (Within ROPP this option
would not be needed). The routine could then be a true transform in the
agreed definition.

Summary: 
1) Yes, the comments could be improved in case someone wishes to use
this routine stand-alone (and in particular that the input datum is
preserved, and how to change from ellipsoid to geoid), and possibly to
upgrade the code to do the complete transform internally as a user-
option.

2) ROPP has no need for an explicit ellipsoid/geoid transformation
(undulation calculation) since (a) all geometric altitudes are provided
in geoid space and (b) a pre-calculated value is provided by all
suppliers should the reverse operation be needed. But if a user should
have need of such a calculation, we can supply PES code separately (as
we have for Axel). Or if there is enough demand, bundle it within
ropp_tools - though as we've been busy taking unused code _out_ of the
package, I don't favour this route!

3) When used in the intended context of ROPP, the routine gives the
required output (as long as the input meets the spec), so the outputs
from ROPP are correct.

> So, since there is no EGM96 (or undulation) in ROPP it is impossible
> to calculate the geopotential-height. 

We rely on the undulation provided in the input (notably in the thinner
when generating nominal refrac heights from IP values), but otherwise
assume the undulation is already applied at source and all refrac hight-
based coordinates are wrt geoid, so no explicit calculation is
necessary.

> Therefore, a routine
> 'geometric2geopotential' confuses me (and others...)
> since it can clearly not be calculating the geopotential (height).

It is, but not on the 'correct' datum if the input is not wrt geoid.

> Whether the routine 'geometric2geopotential' in fact can be
> used in a 'smart' way such that it outputs a quantity
> that is close to the geopotential height is another story
> (and I believe I do understand how to use the routine that way).

I think ROPP uses it that way, but out of context, I now understand the
confusion.

Dave



PS: Hans also comments:

> I think it is OK to treat the geopotential heights the way
> it is currently done in ROPP. However, it should be stated
> very clearly in the header of the Fortran source file that
> the geopotential heights are above the ellipsoid and not
> above the geoid. Now the header only says that the function
> returns 'geopotential heights' which per (standard) definition
> is above the geoid. Most users will probably get this wrong
> if they only read the header and not the formulas actually
> implemented.

I agree that if taken out of the ROPP context, geometric2geopotential()
routine isn't actually doing the end-to-end transformation, but whatever
datum goes in, comes out (though the comments suggest an expectation of
the datum being wrt ellipsoid).

To clarify exactly what's happening, we should certainly improve the
header comments. The routines in ROPP_UTILS are documented as being
utility routines for use in other ROPP modules, and are not intended to
be used stand-alone (though they could be) so they have not been so
comprehensively documented as routines that we do expect users to call
via the API (e.g. top-level I/O, FM & 1dvar interfaces).

This one is probably a case where the comments should be clearer.

Kent:

Dave,
thanks for the explanation! This was very instructive.

I just talked to Hans and with this history we now better
understand what confused us. We think it would be helpful
with a few clarifications to the text; to sum up our suggestions
(see also my def's below):

1. name of the routine: ok

2. clarify in the text in the routine:
it outputs the _difference_ in geopotential heights between
the point in question and the point above which the point
in questions' height is measured relative to:
i) for the input quantity h_g the output is _the_ geopotential height, z:
z=A[h_g]
ii) for the input quantity h the output A[h] is z-z_WGS where z_WGS is
the geopotential height of the given point on the ellipsoid

3. also put in a note in the routine (and possibly also at other 
places in ROPP):
when used with input h, and with z'=A[h]-u, the difference between z' and
z is very small (since z_WGS is very nearly equal to -u).

What do you think of these changes to the text?


Kent

Definitions:
h: height for some point P above WGS84
h_g: height for some point P above EGM96
u: undulation 
A: the algorithm/routine in question

[Huw Lewis]

Code headers updated and committed to trunk (see [2104]). The present implementation does not make reference to undulation, as this quantity is not required in most implementations within ROPP.

Ticket closed as fixed.