Reprojecting XYZ coordinates in Python whilst keeping Z coordinate constant

Question

When doing a reprojection from a projected coordinate system (e.g. EPSG:28992) to a geographic coordinate system (e.g. EPSG:4326), the Z coordinate gets an offset of 45 meters. The Z coordinates of both coordinate systems are in meters and have 0 set at mean sea level. So when reprojecting XYZ points with a height of 0 meters, you expect them to have a height of 0 meters after reprojection. This is not the case in Geopandas (which uses PyProj), PyProj and GDAL. How can one reproject XYZ coordinates with Python whilst keeping Z coordinate constant?

import geopandas as gpd
from shapely.geometry import Point

geometry = [Point(192080.218800001, 326727.25, 0.0)]
crs = {'init': 'epsg:28992'}
pointRDnew = gpd.GeoDataFrame(index=[0], crs=crs, geometry=geometry)
pointWGS84 = pointRDnew.to_crs(epsg="4326")
print(pointWGS84.geometry[0])

POINT Z (5.914611042609124 50.92907898070797 45.24482111539692)

from pyproj import Proj, transform

x1, y1, z1 = 192080.218800001, 326727.25, 0.0
inProj = Proj(init='epsg:28992')
outProj = Proj(init='epsg:4326')
x2, y2, z2 = transform(inProj, outProj, x1, y1, z1)
print(x2, y2, z2)

5.914611042609124 50.92907898070797 45.24482111539692

from osgeo import ogr, osr

wktGeometry = 'POINT (192080.218800001 326727.25 0.0)'
epsgTransformation = epsgToWGS84Transformation(28992)
geometryFromWkt = ogr.CreateGeometryFromWkt(wktGeometry)
geometryFromWkt.Transform(epsgTransformation)
print(geometryFromWkt)

POINT (5.91461112637604 50.9290789739909 45.2474540863186)

Answer 1

The transformation that you're using is probably an "equation-based" method like geocentric translation, coordinate frame/position vector, or Molodensky-Badekas. The methods will convert the z values. However, the z value must be ellipsoidal heights, not gravity-related heights which is what it sounds like you have.

You might be able to switch the transformation to use an NTv2 file (rdtrans2008.gsb), EPSG::7000, which should ignore any z values.

Or, don't send the z values or overwrite the output z values with the input ones.

I don't use geopandas or gdal so I can't tell you what over workarounds exist.

Disclosure: I work at Esri and help maintain the EPSG registry.