Transforming EPSG:4979 to given WKT (~EPSG:6654)

Question

I am setting up a GIS to transform points defined in lon/lat/z (referenced to WGS84, or in EPSG:4979) to an input .laz/.las file coordinate system. I've handled reading the wkt from laz/las files with pdal just fine. The problem is creating the transform via pyproj and maintaining projection complexity in 3D.

My pyproj version is 3.2.1

Suppose the input WKT translates precisely to EPSG:6654, or "NAD83(CSRS) / UTM Zone 11N + CGVD2013 height". I transform a point p = [-117,50,0] to the new coordinate system:

import numpy as np
import pyproj as proj
from pyproj import CRS

p = np.array([-117, 50, 0]) # lon, lat, z
print(p)

crs_from1 = CRS('EPSG:4326').to_3d() # 4326 -> 4979
crs_to1 = CRS('EPSG:6654').to_3d() # 6654 -> 6654, same
tr1 = proj.Transformer.from_crs(crs_from1, crs_to1, always_xy=True)
vec1 = tr1.transform(p[0],p[1],p[2])

print(vec1)

Prints out:

[-117   50    0]
(500001.03747013933, 5538629.990289256, 0.0)

This is odd because z is exactly the same, despite point p being defined relative to the WGS84 ellipsoid and crs_to1 is defined relative to the GRS80 ellipsoid. Removing CGVD2013, I compare to a related transformation, EPSG:4979 -> EPSG:2955 "NAD83(CSRS) / UTM Zone 11N":

import numpy as np
import pyproj as proj
from pyproj import CRS

p = np.array([-117, 50, 0]) # lon, lat, z
print(p)

crs_from2 = CRS('EPSG:4326').to_3d() # 4326 -> 4979
crs_to2 = CRS('EPSG:2955').to_3d() # 2955 -> None
tr2 = proj.Transformer.from_crs(crs_from2, crs_to2, always_xy=True)
vec2 = tr2.transform(p[0],p[1],p[2])

print(vec2)

Prints out:

[-117   50    0]
(500001.03747013933, 5538629.990289256, 0.4121707445010543)

In the second method, the only change is to remove CGVD2013.

If pyproj is not made aware that the CRS is "3D", it does not transform the z-axis coordinate. Is this what is happening for tr1 (even though I manually specified it)? I think I'm doing this wrong because it appears that the CGVD2013 vertical datum exactly undoes projection errors between WGS84 and NAD83. Does anyone know how to transform between EPSG:4979 and EPSG:6654?

Answer 1

It is likely that you are missing the transform grids. https://pyproj4.github.io/pyproj/stable/transformation_grids.html

You can either download the grids or enable the network. https://pyproj4.github.io/pyproj/stable/api/network.html

>>> import pyproj
>>> pyproj.network.is_network_enabled()
True

Also, I tried it out with PROJ 8.2.0 which is included in the wheels and it does not work. With a PROJ 9+ it seems to work. So, a patch was likely included for that version.

from pyproj import CRS, Transformer
crs_from1 = CRS("EPSG:4979")
crs_to1 = CRS("EPSG:6654")
trans1 = Transformer.from_crs(crs_from1, crs_to1, always_xy=True)

points = [-117, 50, 0]

print(trans1.transform(*points))

Output:

(500000.0, 5538630.702744345, 15.12400007247922)

Version Information:

>>> import pyproj
>>> pyproj.show_versions()
pyproj info:
    pyproj: 3.3.1
      PROJ: 9.0.0
  data dir: ~/miniconda/envs/midas/share/proj
user_data_dir: ~/.local/share/proj
PROJ DATA (recommended version): 1.9
PROJ Database: 1.2
EPSG Database: v10.054 [2022-02-13]
ESRI Database: ArcMap 12.9 [2022-02-18]
IGNF Database: 3.1.0 [2019-05-24]