Pyproj – Transform EPSG:4979 to Given WKT (~EPSG:6654)

nad83pyprojutmvertical-transformationwgs84

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?

Best Answer

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]

Related Solutions

How to Project EPSG:4326 Data on a 2D Map

OpenLayers uses the term 'EPSG:4326' to mean the Plate Caree projection. Referring to 'WGS84' and EPSG:4326 as a projection has been common for so long that it is a source of confusion. This short-hand has been going on since before Google and OpenLayers came on to the scene. For instance, ESRI have been fudging the terms for as long as I can remember. OpenLayers does not reproject to EPSG:900913 on the fly unless you tell it to, which you have to if you want to mix data with Google base maps because the Google API only uses 900913 (which was invented by Google - the numbers being vaguely reminiscent of the word 'google').

PostGIS – Is Geometry(4326) an Alias for Geography Datatype?

tl;dr: No. But see the update at the end.

GEOMETRY and GEOGRAPHY are different PostgreSQL TYPES; exposed as higher level (SQL), user-defined composite types to the psql environment, but implemented as lower level (C) base types.

These base types are defined on the C level with e.g. their own typemod constraints, TOAST support, data I/O functions, operator classes and casting behaviour. On top of that, both types have their own indexing mechanics.

On the SQL level, both types (optionally) accept a TEXT and a NUMERIC parameter, which translate to geometry type, and spatial reference. Additionally, you can index both with the same command (USING GIST(<geom>)), although internally those indexes differ, and CAST between them almost seamlessly, with some caveats as seen in the example below.

In every way, these are different types with fundamentally different mechanics, despite their deliberately similar handling.

Under the hood, the core difference is the mathematical representation and calculation of and with the given coordinates, in two words: planar vs. spheroidal.

The GEOMETRY type is bound by definition to a 2D planar reference, and operates with Cartesian math regardless of the nature of the CRS; distance between two points is a straight line as e.g. per Pythagoras, even if your geometries are defined in a geodetic reference (e.g. EPSG:4326). The GEOMETRY type is meant to work best for projected CRS.

The GEOGRAPHY type is inherently assuming that the given coordinates are geodetic coordinates, thus operating solely on great circle math (spherical) and/or it's more complicated application on spheroidal/ellipsoidal bodies. The GEOGRAPHY type is meant to work with geographic/geodetic RS (Longitudes/Latitudes) only.

Neither does PostGIS "the same 'Geography calculations'", nor detects the nature of a given SRS to do so.

The WGS84 ellipsoid (EPSG:7030) is allegedly among the most usable representation of our planet on a global scale, and is widely supported in transformations to and from projections and other ellipsoids. Also, it is and has ever since been the base reference of most of our positioning systems.

It is de facto a global standard for spatial data storage and its interchangeability.

In PostGIS, the actual default SRS is no SRS, or 0 (which is bad to have...).

Now, as a small example of what PostGIS does with wrong SRS definitions; you can do this (note that the coordinates should represent degrees...)

SELECT ST_AsText(ST_SetSRID('POINT(1000 1000)'::GEOMETRY, 4326));

returning a valid GEOMETRY

st_astext       |
----------------|
POINT(1000 1000)|

while this

SELECT ST_AsText(ST_SetSRID('POINT(1000 1000)'::GEOGRAPHY, 4326));

returns a valid GEOGRAPHY

st_astext     |
--------------|
POINT(-80 -80)|

along with the note

Coordinate values were coerced into range [-180 -90, 180 90] for GEOGRAPHY

Btw. coercing here means sth. like

Lat: -90 + (1000 % 90) = -80
Lon: -180 + (1000 % 180) = -80

so it effectively walks around the globe until 1000 degrees are traversed.

Yet the GEOGRAPHY type does not care if you do

SELECT ST_AsText(ST_SetSRID('POINT(1000 1000)'::GEOGRAPHY, 3857));

except for the above note of coordinate range coercion, while it actually should (if it could; we're assigning a SRS on the created type, so it would be the function assigning it that needed to) inform you about the nature of the SRS (and, since we're at it, IMHO also the uselessness of EPSG:3857...).

However, a CAST to GEOGRAPHY will indeed check the SRS and denies the cast for non-geodetic SRS, e.g.

SELECT ST_SetSRID('POINT(0 0)'::GEOMETRY, 3857)::GEOGRAPHY;

will error out with

ERROR: Only lon/lat coordinate systems are supported in geography.

In no way does PostGIS any fail-saving for you in these matters, except when using a CAST to GEOGRAPHY.

Update as per comment request:

Because PostGIS defines the cast from GEOMETRY to GEOGRAPHY (and only this way) AS IMPLICIT, PostgreSQL is free to cast any GEOMETRY to GEOGRAPHY whenever it is necessary, or possible, to solve ambuiguity. This affects parameters to functions in the following ways, where PostgreSQL will CAST implicitly:

if a distinct signature is unambiguously invoked that accepts only GEOGRAPHY types; an example would be the ST_Area signatures where, when explicitly used with a value for use_spheroid, PG casts a GEOMETRY to GEOGRAPHY
an ambiguous signature for both types is used (i.e. a type overloaded function), and only one parameter is passed in as GEOGRAPHY, the other will get cast; an example is ST_DWithin where, when one parameter gets passed as GEOMETRY and the other as GEOGRAPHY, the former will get cast to GEOGRAPHY

One might say, since only the CAST to GEOGRAPHY is IMPLICIT, it has inherent precedence.

In these cases, since the CAST is used, PostGIS will inform you (and deny the execution) if the GEOMETRY in question does not use a geodetic SRS!

The secret ingredient to ambiguity is an implicit CAST.

Best Answer

Related Solutions

How to Project EPSG:4326 Data on a 2D Map

PostGIS – Is Geometry(4326) an Alias for Geography Datatype?

Related Question