Including geoidgrids when initializing projection via EPSG

Question

It's not clear to me how you're intended to initialize a pyproj.Proj() object with both an EPSG code and geoidgrids.

For my specific use case with the grid file downloaded from here:

epsg = 3162  # NAD83(CSRS) / Ontario MNR Lambert
grid_path = CGG2013an83.gtx  # CGVD2013(CGG2013a) height --> NAD83(CSRS)

Using all keyword args works, but init=epsg:xxx triggers a deprecation warning:

proj = pyproj.Proj(init=f'EPSG:{epsg}', geoidgrids=grid_path)
>>>    FutureWarning: '+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method.
>>>      return _prepare_from_string(" ".join(pjargs))
proj.crs.is_vertical
>>> True

Using the suggested syntax appears to make it silently ignore geoidgrids (and in fact all other keyword arguments):

proj = pyproj.Proj(f'EPSG:{epsg}', geoidgrids=grid_path)
proj.crs.is_vertical
>>> False

I'm not seeing any relevant examples in the docs.

Using pyproj 2.4.2

Answer 1

Here is the reasoning for the deprecation warning: https://github.com/pyproj4/pyproj/issues/355

However, you can probably still safely use it for a while as it is likely not going to be removed for a while in PROJ: https://github.com/OSGeo/PROJ/issues/1895

But, if you would like to avoid the deprecation warning, you can either filter the warning out or your path may be a little more complicated.

https://pyproj4.github.io/pyproj/stable/build_crs.html

When you use the +init= syntax, it will build out a lot of stuff for you:

from pyproj import CRS

crs = CRS(init="epsg:3162", geoidgrids="CGG2013an83.gtx")
crs

<Compound CRS: +init=epsg:3162 +geoidgrids=CGG2013an83.gtx +type= ...>
Name: unknown
Axis Info [cartesian|None]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- undefined
Datum: Unknown based on GRS80 ellipsoid
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich
Sub CRS:
- unknown
- unknown

However, without the +init= syntax, you get only part of what you want:

projected_crs = CRS("EPSG:3162")
projected_crs

<Projected CRS: EPSG:3162>
Name: NAD83(CSRS) / Ontario MNR Lambert
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- name: Canada - Ontario
- bounds: (-95.16, 41.67, -74.35, 56.9)
Coordinate Operation:
- name: Ontario MNR Lambert
- method: Lambert Conic Conformal (2SP)
Datum: NAD83 Canadian Spatial Reference System
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich

So, this is how you manually build what PROJ built for you:

from pyproj import Proj
from pyproj.crs import BoundCRS, CompoundCRS,  GeographicCRS, VerticalCRS
from pyproj.crs.coordinate_system import Ellipsoidal3DCS
from pyproj.crs.enums import Ellipsoidal3DCSAxis

vertical_crs = VerticalCRS(name="unknown", datum='VDATUM["unknown"]')
geographic3d_crs = GeographicCRS(
    name="WGS 84",
    ellipsoidal_cs=Ellipsoidal3DCS(
        axis=Ellipsoidal3DCSAxis.LATITUDE_LONGITUDE_HEIGHT
    ),
)
bound_crs = BoundCRS(
    source_crs=vertical_crs,
    target_crs=geographic3d_crs,
    transformation={
      "$schema": "https://proj.org/schemas/v0.2/projjson.schema.json",
      "type": "Transformation",
      "name": "unknown to WGS84 ellipsoidal height",
      "source_crs": vertical_crs.to_json_dict(),
      "target_crs": geographic3d_crs.to_json_dict(),
      "method": {
        "name": "GravityRelatedHeight to Geographic3D"
      },
      "parameters": [
        {
          "name": "Geoid (height correction) model file",
          "value": "CGG2013an83.gtx",
          "id": {
            "authority": "EPSG",
            "code": 8666
          }
        }
      ]
    }
)

compound_crs = CompoundCRS(
    name="unknown",
    components=[projected_crs, bound_crs],
)
compound_crs

<Compound CRS: {"$schema": "https://proj.org/schemas/v0.2/projjso ...>
Name: unknown
Axis Info [cartesian|None]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- undefined
Datum: NAD83 Canadian Spatial Reference System
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich
Sub CRS:
- NAD83(CSRS) / Ontario MNR Lambert
- unknown

proj = Proj(compound_crs)
proj.crs.is_vertical

True

Hopefully this helps.