2

I am trying to correct ellipsoid heights to orthometric height, transforming from WGS84 to EGM2008 using the EPSG code 3855. I have tried with pyroj, as in this question:

import pyproj
EGM2008height = pyproj.Proj("EPSG:3855") 
wgs84_ll = pyproj.Proj("EPSG:4326")
print(pyproj.transform(wgs84_ll, EGM2008height, 34.68016909181368, 38.31245226053967,100))

which returns the error:

(Internal Proj Error: proj_create: unrecognized format / unknown name)

I understand that I need the egm08_25.gtx file, which I have downloaded but I don't know how to use it so with pyproj. Is there a way to add EPSG codes to pyproj?

I can get it to work by looping over my data and using gdaltransform by calling the command:

cmd = 'echo ' + '"' + str(X[i]) + ' ' + str(Y[i]) + ' ' + str(Z[i]) + '"' + ' | gdaltransform -s_srs "+proj=longlat +datum=WGS84 +no_def" -t_srs "+proj=longlat +datum=WGS84 +no_defs +geoidgrids=/Usegm08_25.gtx"'

But this is a slow process and I think using pyproj will be more efficient.

There are several questions on this already, but none of them comprehensively show how to use the .gtx file or add EPSG codes into pyproj.

EDIT.....

Related Qs where I think the answe requires reviewing:

Vertical Datum transformation using Pyproj

Converting Ellipsoidal height to Orthometric height in Python

Converting EPSG:2263 to WGS84 using Python Pyproj

1
  • 2
    What are links to those previous questions? It will be useful to identify them so that they and/or their answer(s) can be reviewed.
    – PolyGeo
    Jul 3 '20 at 18:39
4

With Pyproj 2.x, it's better to use the Transformer class, as below :

from pyproj import Transformer

transformer = Transformer.from_crs("epsg:3855", "epsg:4326")
print(
    transformer.transform(34.68016909181368, 38.31245226053967, 100)
)

and the output is :

(34.68016909181368, 38.31245226053967, -100.0)
4
  • Awesome, thanks! But surprised this was not so clear in previous questions on this... Jul 3 '20 at 18:27
  • Note: -100 isn't the correct answer. Charles Karney's online geoid calculator gives a difference of 35.2833 m (don't ask me the sign) assuming 34 is the latitude.
    – mkennedy
    Jul 12 '20 at 2:27
  • @mkennedy, it looks like you got the lat and lon backwards. that geoid height calculator estimates 22.2315m with the 34 as latitude. your point is still valid though. the Epsg transformer didn't seem to do anything Jul 18 '20 at 4:36
  • the Unavco Geoid Height Calculator reports -23.107m for these coordinates assuming 0m GPS elevation. The numbers are close in magnitude but sign is a bit unclear Jul 18 '20 at 4:36

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.