Reprojecting XYZ coordinates in Python whilst keeping Z coordinate constant

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)`

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.

• Hi Mkennedy. Thanks for the advice to switch transformation. Both Z values are actually ellipsoidal heights, but EPSG:28992 is in Bessel ellipsoid and EPSG:4326 is in WGS84 ellipsoid (en.wikipedia.org/wiki/Bessel_ellipsoid). Even though both heights are set at mean sea level, the calculation differs: Bessel ellipsoid 1841 (defined by log a and ƒ): a = 6,377,397.155 m ƒ = 1 / 299.1528153513233 (0.003342 773154 ± 0.000005)[citation needed] b = 6,356,078.963 m. Earth ellipsoid WGS84 (defined directly by a and ƒ): a = 6,378,137.0 m ƒ = 1 / 298.257223563 b = 6,356,752.30 m. Aug 8, 2018 at 6:06
• Both ellipsoids return different `ellipsoidal heights`, which are calculated from the non-aligning ellipsoids, not mean sea level. PROJ/GDAL can not handle mean sea level heights, so using a ntv2 grid is the only way to go for you. Aug 8, 2018 at 7:30