PostGIS ST_Transform seems to treat the vertical datum but not the equivalent pyproj Transformer, even so they both rely on Proj

Question

I've narrowed down a problem I met when processing some points having altitude values to the following:

It seems that PostGIS and pyproj do not treat the vertical datum in the same way but both are actually based on Proj.

Detailed explanation

Here's a simple reproducible PostgresSQL query and its response:

SELECT 
  ST_AsText(
    ST_Transform(
      ST_SetSRID(
        ST_MakePoint(8.37909, 47.01987, 1000),
      4326),
    2056)
  ) AS geom;

Here, the altitude of 1000m seems to be evaluated as an ellipsoidal height (relative to the WGS84 ellipsoid I guess?)
Hence =>

geom (text)
"POINT Z (2671499.8913080636 1208075.1135782269 951.4265527641401)"

And here is the same transformation in Python using pyproj:

import numpy as np
from shapely import ops as sp_ops
from shapely.geometry import Point
import pyproj
from pyproj import transform, Proj, Transformer

project = Transformer.from_crs('EPSG:4326', 'EPSG:2056', always_xy=True)
point = Point(8.37909, 47.01987, 1000)
geom = sp_ops.transform(project.transform, point)

np.set_printoptions(suppress=True)
print(np.array(geom.coords))

which results in:

[[ 2671499.8913080636   1208075.1135782297    1000.  ]]

As you can see here, the altitude stays the same after the transformation. I also tried using the EPSG code for the true 3D WGS84, namely EPSG:4979 according to: Vertical Datum transformation using Pyproj and https://epsg.org/crs_4979/WGS-84.html?sessionkey=00zgez8kt5 but the result is exactly the same.

At this stage, you can notice that PostGIS' result is actually "accurate" according the local geoid-ellipsoid separation which is roughly 48.21m at this position, so 1000 (ellipsoidal_height) - 48.21 (separation) ≅ 951.79 (orthometric_height).

So, why does the altitude stay the same when using pyproj and not when using PostGIS when 1) both actually rely on Proj itself and 2) the EPSG codes I used are the same in both cases?
May it be that PostGIS is a little to lenient?

It is said on the ST_Transform documentation that:

Requires PostGIS be compiled with Proj support. Use PostGIS_Full_Version to confirm you have proj support compiled in.

and when I actually check for its full version, it's aligned with that statement:

SELECT PostGIS_Full_Version();

"POSTGIS=""3.1.0 5e2af69"" [EXTENSION]     
  PGSQL=""130""    
  GEOS=""3.7.1-CAPI-1.11.1 27a5e771""
  PROJ=""Rel. 5.2.0, September 15th, 2018""
  LIBXML=""2.9.4""
  LIBJSON=""0.12.1""
  LIBPROTOBUF=""1.3.1""
  WAGYU=""0.5.0 (Internal)""    
TOPOLOGY"

On the other hand, my pyproj version is:

>>> pyproj.__version__
'3.0.0.post1'

But I cannot figure out any reference to altitude on vertical datum in the pyproj documentation https://pyproj4.github.io/pyproj/stable/api/transformer.html

I would be very grateful if someone had a clear explanation as to why things are currently the way they are. And how could I properly handle the case in Python, i.e. the same way PostGIS does?

Answer 1

One thing to note is the version of PROJ that is being used. PostGIS in your version is using PROJ 5.2 and the version of pyproj you are using is likely using PROJ 7.2+.

There are quite a few things that have changed since then. If you go back to PROJ 5.2 with pyproj 1.9.6 you get:

>> import pyproj
>>> pyproj.__version__
'1.9.6'
>>> pyproj.proj_version_str
'0.5.20'
>>> from pyproj import Proj, transform
>>> transform(Proj(init="epsg:4326"), Proj(init="epsg:2056"), 8.37909, 47.01987, 1000)
(2671499.8913080636, 1208075.1135782269, 951.4265527641401)

However, with the latest PROJ/pyproj you get:

>>> pyproj.__proj_version__
'8.0.0'
>>> pyproj.__version__
'3.0.1'
>> project = pyproj.Transformer.from_crs('EPSG:4979', 'EPSG:2056', always_xy=True)
>>> project.transform(8.37909, 47.01987, 1000)
(2671499.8913080636, 1208075.1135782297, 1000.0)
>>> project
<Concatenated Operation Transformer: pipeline>
Description: axis order change (2D) + Inverse of CH1903+ to WGS 84 (1) + Swiss Oblique Mercator 1995
Area of Use:
- name: Liechtenstein; Switzerland.
- bounds: (5.96, 45.82, 10.49, 47.81)

I suspect it has something to do with the fact that your destination CRS is 2D:

>>> pyproj.CRS("EPSG:2056")
<Projected CRS: EPSG:2056>
Name: CH1903+ / LV95
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- name: Liechtenstein; Switzerland.
- bounds: (5.96, 45.82, 10.49, 47.81)
Coordinate Operation:
- name: Swiss Oblique Mercator 1995
- method: Hotine Oblique Mercator (variant B)
Datum: CH1903+
- Ellipsoid: Bessel 1841
- Prime Meridian: Greenwich

The results from PROJ 8 are the same:

echo 47.01987 8.37909 1000 | cs2cs -d 8 'epsg:4326' 'epsg:2056'
2671499.89130806    1208075.11357823 1000.00000000
$ cs2cs
Rel. 8.0.0, March 1st, 2021

Sounds like something to bring up here: https://github.com/osgeo/PROJ/

UPDATE:

PROJ Information: https://github.com/OSGeo/PROJ/issues/2609

pyproj 3.1 update: https://github.com/pyproj4/pyproj/pull/808

    >>> from pyproj import CRS, Transformer
    >>> transformer_3d = Transformer.from_crs(
    ...     CRS("EPSG:4326").to_3d(),
    ...     CRS("EPSG:2056").to_3d(),
    ...     always_xy=True,
    ...)
    >>> transformer_3d.transform(8.37909, 47.01987, 1000)
    (2671499.8913080636, 1208075.1135782297, 951.4265527743846)