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

GeoMonkey
  • 1,357
  • 11
  • 26
  • 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 03 '20 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 04 '23 at 08:44
  • Do you have proj-data installed? – GeoMonkey Aug 04 '23 at 17:32

2 Answers2

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

Scott Staniewicz
  • 230
  • 1
  • 4
  • 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). – Robert Calhoun Nov 14 '23 at 15:16
  • 1
    See https://gis.stackexchange.com/a/470274/6247 for how to download the necessary EGM 2008 grid. – Robert Calhoun Nov 14 '23 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)
J. Monticolo
  • 15,695
  • 1
  • 29
  • 64
  • Awesome, thanks! But surprised this was not so clear in previous questions on this... – GeoMonkey Jul 03 '20 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 '20 at 02: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 – Marc Compere Jul 18 '20 at 04: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 – Marc Compere Jul 18 '20 at 04: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) – Scott Staniewicz Aug 04 '22 at 18:12
  • @ScottStaniewicz: thanks for this useful comment and your solution. – J. Monticolo Aug 04 '22 at 18:45