2

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)

2

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. – Davma Aug 8 '18 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. – AndreJ Aug 8 '18 at 7:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.