Pyproj – Vertical Datum Transformation Techniques

pyprojvertical-transformation

I have been using Pyproj for some time now doing only simple projection and transformation. I am totally new to vertical datum and differences in GEOID among Coordinate systems.
I have been trying to do vertical transformation using pyproj but I feel like I am missing something.
For example I tried this : 

import pyproj
lat = 43.70012234
lng = -79.41629234
z = 100
wgs84 = pyproj.Proj("+init=EPSG:4326")
NAD83NDV88 = pyproj.Proj("+init=EPSG:5498")
results = pyproj.transform(wgs84, NAD83NDV88, lng, lat, z)

which returns : (43.70012234, -79.41629234, 100.0).

My logic was to take a WGS84 coordinate with an ellipsoid elevation and try to convert that into NAD83 NADV88 (epsg:5498). But it seems like it is not that simple.

I found some EPSG codes for transformation related to this question like : EPSG:3858 +> WGS 84 to EGM2008 height (1) but I haven't been able to figure out how to use those codes in Pyproj. I am missing some pieces of the puzzle!

So to summarize can this be done?

Best Answer

The accepted answer has a couple of flaws I'd like to address.

First, the objection to the use of NAD83 is misplaced. It is true that NAD83 is a horizontal datum, but the setup in the question uses EPSG:5498 -- a compound coordinate reference system (CRS) that combines both NAD83 (EPSG:4269) as a horizontal CRS and NAVD88 (EPSG:5703) as a vertical CRS.

However, this concern is valid when looking at the particular form of WGS84 being used -- EPSG:4326. Looking at the linked info page, we see that it uses the WGS84 datum, but a 2D coordinate system:

COORDINATE SYSTEM: Ellipsoidal 2D CS. Axes: latitude, longitude. Orientations: north, east. UoM: degree

Note there is no axis for height. This is the 2D geographic form of WGS84. If we instead use EPSG:4979 -- the 3D geographic form of WGS84 -- we see a difference:

COORDINATE SYSTEM: Ellipsoidal 3D CS. Axes: latitude, longitude, ellipsoidal height. Orientations: north, east, up. UoM: degree, degree, metre.

Now the conversion makes more sense, and can be summarized as:

Parameter Source Destination
Horizontal Datum WGS84 NAD83
Horizontal Coordinates geographic geographic
Vertical Datum WGS84 ellipsoid NAVD88
Vertical Coordinate height above ellipsoid height above geoid

Second, given that the question was scoped to pyproj, it would be best to come up with a more definitive answer in that environment rather than changing over to PROJ.

So it looks like you are getting out the same coordinates that you put in. This can be the result of not having the right support data in place, namely grid files for vertical datums. There are some diagnostic tests that can be done to understand what pyproj (and/or the underlying PROJ code) is trying to do in a transformation. However, it is easier to use if you switch from the PROJ-string way of working with pyproj, to the object-oriented way. So I will rewrite your code snippet as:

import pyproj
lat = 43.70012234
lng = -79.41629234
z = 100
wgs84 = pyproj.crs.CRS.from_epsg(4979)
nad83_navd88 = pyproj.crs.CRS.from_epsg(5498)

Now, rather than immediately doing a one-off transform, we will try to create a Transformer object that we can keep around and investigate:

tform = pyproj.transformer.Transformer.from_crs(crs_from=wgs84, crs_to=nad83_navd88)

If you look at the value of tform, and you have an issue with missing grid files, you'll see something like the following:

(At this point I will have to deviate from your example slightly because my install is doing the NAD83+NAVD88 conversion correctly. I will use WGS84+EGM84 instead. So you will see that in the outputs I give.)

<Other Coordinate Operation Transformer: noop>
Description: Inverse of Transformation from EGM84 height to WGS 84 (ballpark vertical transformation, without ellipsoid height to vertical height correction)
Area of Use:
- undefined

The important words here are noop (meaning "no operation", computer science jargon for a function that does nothing), and ballpark vertical transformation. If pyproj/PROJ doesn't have the necessary grid files for a vertical datum transformation, it will simply give you back the input you gave as a very naive approximation.

If you want to see exactly what file is missing, you can use a TransformerGroup object, following these instructions in the pyproj documentation. NOTE: some datums use more than one file to define their grid. This method won't be able to tell you which of multiple files is missing.

After you learn what's missing, you can refer to the Transformation Grids section to learn how to find the data directory and put the right files in it.

Also, to address this really quick:

I found some EPSG codes for transformation related to this question like : EPSG:3858 +> WGS 84 to EGM2008 height (1) but I haven't been able to figure out how to use those codes in Pyproj. I am missing some pieces of the puzzle!

Like you said, that EPSG code (3858) belongs to a transformation. This is what the pyproj Transformer object represents, and when you ask it for a Transformer, it is internally looking up the correct EPSG transformation definition. But as far as I know, there is no way to directly access a transformation using its EPSG code in pyproj -- you have to access it indirectly by telling pyproj which CRSes you want to convert.

Related Question