2

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

2
  • Possible to give a specific example?
    – snowman2
    Commented Feb 27, 2020 at 22:05
  • @snowman2 sure, see the edit
    – mikewatt
    Commented Feb 27, 2020 at 22:43

1 Answer 1

3

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.

3
  • Thanks for the detailed explanation and confirming I wasn't missing some other straightforward way. I am filtering the warning which feels a bit dirty, but given the alternative I think I'll stick with it for now. IMO it would make sense to actually enforce the mutual exclusion of projparams and **kwargs in the Proj() class as I can't be the first person to be confused by it silently ignoring arguments, but maybe I'm again missing the reasoning. You sound like a maintainer so might as well throw it out there ;)
    – mikewatt
    Commented Feb 28, 2020 at 2:04
  • Not a bad idea. Mind making an issue on GitHub?
    – snowman2
    Commented Feb 28, 2020 at 2:26
  • Sure, you got it
    – mikewatt
    Commented Feb 28, 2020 at 2:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.