0

I'm trying to calculate the orthometric height between ITRF97(or EGM96) and wgs84 (or grs80) in specific latitudes and longitudes.

For example: in longitude 57.111 and latitude -5.582 i used an online calculator and the result was 55.380 [m], which seems to be correct.

Is there a simple way to calculate this height? I am using Python for the calculation. Is there a module which converts between geoid and ellipsoid? If so how can I use it?

1
0

First, we should clarify some points about the question,

ITRF and WGS84 are Geodetic Reference Frames (a.k.a. "horizontal datums"), Which are calculated by different organizations, but are in general compatible.

The realization strategy of both systems differs. One key difference is that ITRF is recalculated from time to time, giving its "realizations" named ITRF, i.e.: ITRF97, ITFR2000, etc. Meanwhile, WGS84 is not usually renamed nor fully recalculated, instead it is periodically aligned to ITRF to ensure consistency with Geodetic mesurement. (see: this answer and its references in the comments).

It is probably not your use case, but there may be some differences in positions due to different epoch of materialization (the point in time for the coordinates which are published) of the reference frame, because tectonic dynamics must be considered when high precision is needed.

EGM96 is a geoid model, which is a model of the behavior of the gravity on the earth surface. As a derived product, equipotential surfaces of the earth can be computed from EGM96. More specifically, the equipotential surface of the earth which best fits to the mean level of the seas, is computed for its usage as a vertical reference surface, and orthometric heights are computed related to this surface. (see this question, and this other -disclaimer, accepted answer is mine-).

With this in mind, your needed conversion is probably from orthometric heights defined with respect to the EGM96 model to ellipsoidal heights with respect to the WGS84/ITRF97 reference frame(s) and its associated ellipsoid(s).

EGM96 is so common, that many libraries and online calculators have its data pre-loaded for computations. This answer should help for python

Edit: Epsg code and example

The EPSG code for EGM96 heights is EPSG:5773 (link here). If you look at the details, the definition includes a grid which contains the relevant conversion: EXTENSION["PROJ4_GRIDS","egm96_15.gtx"]]

By example, the cct tool, included on the PROJ tools can handle this (if you have the datum grids installed).

cct +ellps=WGS84 +proj=vgridshift +grids=egm96_15.gtx <<EOF
< coordinates here >
EOF

Here i have added a gist with an ipynb as example it uses the cct commandline tool behind the scenes from the python notebook, you can open it on colab and make tests.

On pyproj (the python bindings for PROJ) the Transformer class should do the trick: See the docs

Further information

For information and data on Geoid models, a full reference is given on this link

For information on the geodetic reference frames, ITRF publications are usually a good and detailed reference, see the ITRF97 documentation

5
  • thank you very much for the clarification that helped a lot. how can i fine to correct EPSG related to EGM96? in addition i cannot make a vertical shift from WGS84 to EGM96 (i tried several EPSG's i find online but the final result was the same as the original value). i assumed that after vertical shift the latitude and longitude stay the same as the original only the height changes. am i missing something? do you know the python command for such transformation?
    – I d
    Jun 24 at 6:32
  • i just wanted to add i'm using both pyproj and GDAL in python so i'm guessing there are built in functions for that in one of them, but i could not find any working explanation on how to use them for such vertical shift.
    – I d
    Jun 24 at 7:28
  • I've added an example and some extra information. Please accept the answer if it works. Maybe you should improve your question with an example of what you have tried with python or GDAL.
    – Javier JC
    Jun 24 at 13:50
  • i'm probably doing something wrong, but the example that was attached is using the old version of pyproj. i tried to use the new version, with WGS84EPSG:4326, EGM96EPSG:5773. my code was: transformer = pyproj.Transformer.from_crs( {"proj": 'latlong', "epsg":'epsg:4326'}, {"proj": 'latlong', "epsg":'epsg:5773'}, ) lonout,latout,altout = transformer.transform(57.11, -5.58, 0, radians=False) the height differences is 53[m] at that point, so i would expect the output to be 57.11,-5.58,53).
    – I d
    Jun 27 at 6:41
  • I'm not sure if there is a simpler method, or more accurate one but i found the following pygeodesy routine that works fine: stackoverflow.com/questions/22564727/…
    – I d
    Jun 27 at 9:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.