Adding EPSG codes to Pyproj

Question

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

Answer 1

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)

Answer 2

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