4

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

3
  • 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, 2020 at 18:39
  • I am struggling to reproduce your solution. For me the height 100 remains unconverted. Did you use pyproj version 3.3.1 installed with conda or pip?
    – Anni
    Aug 4, 2023 at 8:44
  • Do you have proj-data installed?
    – GeoMonkey
    Aug 4, 2023 at 17:32

2 Answers 2

5
  1. You want to make the "to" CRS a combination: EPSG4326+3855. It still uses the WGS84 ellipsoid for horizontal positions, but it has a new vertical datum (3855).

  2. You'll wanna set the always_xy=True option for pyproj's Transformer class. Then you can always pass in (lon, lat), even when a specific transformation calls for the reverse convention.

from pyproj import Transformer
t = Transformer.from_crs("epsg:4326", "epsg:4326+3855", always_xy=True)
t.transform(34.68016909181368, 38.31245226053967,100)
# (34.68016909181368, 38.31245226053967, 64.7096561282928)

Now the difference is about 35.2 meters, as the commenter pointed out that this site says https://www.unavco.org/software/geodetic-utilities/geoid-height-calculator/geoid-height-calculator.html

(I'm assuming that you have longitude = 34.6, and latitude = 38.31)

I've got pyproj version 3.3.1

2
  • 1
    Well done! There are a bunch of SO posts trying to do this conversion, and your answer shows a way to do it. But I found I needed to set crs_from to epsg 4979 to get a valid conversion under pyproj 3.6.1 (proj 9.3.0). Nov 14, 2023 at 15:16
  • 1
    See gis.stackexchange.com/a/470274/6247 for how to download the necessary EGM 2008 grid. Nov 14, 2023 at 16:53
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)
6
  • Awesome, thanks! But surprised this was not so clear in previous questions on this...
    – GeoMonkey
    Jul 3, 2020 at 18:27
  • 1
    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, 2020 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, 2020 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, 2020 at 4:36
  • 1
    This doesn't work in more recent versions of pyproj. The EPSG:3855 is only the vertical datum, so you can't convert a horizontal (EPSG:4326) to purely a vertical. This code now gives the error: ProjError: Error creating Transformer from CRS.: (Internal Proj Error: proj_crs_get_geodetic_crs: CRS has no geodetic CRS) Aug 4, 2022 at 18:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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