Gravity related height from ellipsoidal height using pyproj

Question

I'm using the pyproj library (http://pyproj4.github.io/pyproj/stable/index.html) to get translate ellipsoidal heights into gravity related height.

---

Code:

from pyproj import CRS
from pyproj import Transformer

crs7416 = CRS.from_epsg(7416)  # dk1 dvr90
crs7417 = CRS.from_epsg(7417)  # dk2 dvr90
crs4326 = CRS.from_epsg(4326)  # ellipsoid 2d 1984
crs4979 = CRS.from_epsg(4979)  # ellipsoid 3d 1984
crs7912 = CRS.from_epsg(7912)  # ellipsoid 3d 1980
crs7789 = CRS.from_epsg(7789)  # geocentric reference

geocentrictrans = Transformer.from_crs(7789, 7417, always_xy=True)
newgeocentrictrans = geocentrictrans.transform(12.568337, 55.676098, 100)
print(newgeocentrictrans)

transformer = Transformer.from_crs(4979, 7417, always_xy=True)
newcoord = transformer.transform(12.568337, 55.676098, 100)
print(newcoord)

Output:

(9451442.134520395, 0.0, -6378079.9229382025)
(347093.2566329163, 6172711.4819890065, 100.0)

The transformation seem to work with the geocentric transformation but for some reason the height isn't transformed when i go from ellipsoidal- to gravity related height. As far as i can understand i am using the correct EPSG-references, but correct me if i'm wrong.

It seems like the pyproj library should be able to handle this but i'm missing some piece of the puzzle.

Answer 1

From projinfo application of the PROJ library:

c:\>projinfo -s EPSG:4979 -t EPSG:7417

returns:

Operation n1:

unknown id, Inverse of DVR90 height to ETRS89 + Inverse of ETRS89 to WGS 84 (1) + UTM zone 33N, unknown accuracy, Denmark - onshore, at least one grid missing

PROJ string:
+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +inv +proj=vgridshift +grids=dvr90.gtx +multiplier=1 +step +proj=utm +zone=33 +ellps=GRS80

[...]  

Grid dvr90.gtx needed but not found on the system. Can be obtained from the proj-datumgrid-europe package at https://download.osgeo.org/proj/proj-datumgrid-europe-latest.zip  

When I try it:

C:\>cct +proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +inv +proj=vgridshift +grids=dvr90.gtx +multiplier=1 +step +proj=utm +zone=33 +ellps=GRS80
proj_create: Error -38: failed to load datum shift filecct: Bad transformation arguments - (failed to load datum shift file)
    'cct -h' for help

I don't know the pyproj library, but the grid is supposed to be necessary there too.

---

A second operation is returned by projinfo:

Operation n2:  

unknown id, Inverse of Transformation from DVR90 height to ETRS89 (ballpark vertical transformation, without ellipsoid height to vertical height correction) + Inverse of ETRS89 to WGS 84 (1) + UTM zone 33N, unknown accuracy, Europe - ETRS89, has ballpark transformation

PROJ string:
+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=utm +zone=33 +ellps=GRS80

[...]

Which is already indicating that it will be done without vertical correction. When I run it:

C:\>cct +proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=utm +zone=33 +ellps=GRS80
55.676098 12.568337 100
  347093.2566   6172711.4820      100.0000           inf 

Where it is seen that it does not perform vertical correction. cct returns inf for the fourth dimension (time) not specified.