I've narrowed down a problem I met when processing some points having altitude values to the following:
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?)
|"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?