2

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.

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.

| improve this answer | |

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.