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?