Accept Datum and Ellipsoid objects for CRS datum= and ellps= constructor params

[zackw]

Suppose you have a standard lon/lat CRS, e.g.

>>> wgs84_crs = pyproj.CRS("epsg:4326")

and you want to create a projected CRS with the same datum. The most obvious way to do that would be to use the 'datum' property of the first CRS as the 'datum' argument to the CRS constructor, but that doesn't work:

>>> aeqd_crs = pyproj.CRS(proj="aeqd", lon_0=-80, lat_0=40.5,
...                       datum=wgs84_crs.datum)
Traceback (most recent call last):
  File "<stdin>", line 2, in <module>
  File "/usr/lib/python3/dist-packages/pyproj/crs.py", line 304, in __init__
    super(CRS, self).__init__(projstring)
  File "pyproj/_crs.pyx", line 1312, in pyproj._crs._CRS.__init__
pyproj.exceptions.CRSError: Invalid projection:
    +proj=aeqd +lon_0=-80 +lat_0=40.5 +datum=World Geodetic System 1984 +type=crs:
    (Internal Proj Error: proj_create: Error -9: unknown elliptical parameter name)

Can you make that work, please? Also for the 'ellps' argument and Ellipsoids.

[zackw]

N.B. The only way I've found to do this, that does work, is

>>> aeqd_crs = pyproj.CRS(proj="aeqd", lon_0=-80, lat_0=40.5,
...                       datum=wgs84_crs.to_dict()['datum'])
>>> aeqd_crs
<Projected CRS: +proj=aeqd +lon_0=-80 +lat_0=40.5 +datum=WGS84 +ty ...>
[...]

but this seems fragile -- how do I know the dictionary will always have a 'datum' parameter? how do I know that that parameter captures all of the necessary information?

[snowman2]

In a future version, you will be able to do this. But, it is a feature currently in progress in the PROJ repo.

>>> from pyproj import CRS
>>> aeqd_crs = CRS(proj="aeqd", lon_0=-80, lat_0=40.5)
>>> print(aeqd_crs.to_json(pretty=True))
{
  "type": "ProjectedCRS",
  "name": "unknown",
  "base_crs": {
    "name": "unknown",
    "datum": {
      "type": "GeodeticReferenceFrame",
      "name": "World Geodetic System 1984",
      "ellipsoid": {
        "name": "WGS 84",
        "semi_major_axis": 6378137,
        "inverse_flattening": 298.257223563
      },
      "id": {
        "authority": "EPSG",
        "code": 6326
      }
    },
    "coordinate_system": {
      "subtype": "ellipsoidal",
      "axis": [
        {
          "name": "Longitude",
          "abbreviation": "lon",
          "direction": "east",
          "unit": "degree"
        },
        {
          "name": "Latitude",
          "abbreviation": "lat",
          "direction": "north",
          "unit": "degree"
        }
      ]
    }
  },
  "conversion": {
    "name": "unknown",
    "method": {
      "name": "Modified Azimuthal Equidistant",
      "id": {
        "authority": "EPSG",
        "code": 9832
      }
    },
    "parameters": [
      {
        "name": "Latitude of natural origin",
        "value": 40.5,
        "unit": "degree",
        "id": {
          "authority": "EPSG",
          "code": 8801
        }
      },
      {
        "name": "Longitude of natural origin",
        "value": -80,
        "unit": "degree",
        "id": {
          "authority": "EPSG",
          "code": 8802
        }
      },
      {
        "name": "False easting",
        "value": 0,
        "unit": "metre",
        "id": {
          "authority": "EPSG",
          "code": 8806
        }
      },
      {
        "name": "False northing",
        "value": 0,
        "unit": "metre",
        "id": {
          "authority": "EPSG",
          "code": 8807
        }
      }
    ]
  },
  "coordinate_system": {
    "subtype": "Cartesian",
    "axis": [
      {
        "name": "Easting",
        "abbreviation": "E",
        "direction": "east",
        "unit": "metre"
      },
      {
        "name": "Northing",
        "abbreviation": "N",
        "direction": "north",
        "unit": "metre"
      }
    ]
  }
}
>>> crs_dict = aeqd_crs.to_json_dict()
>>> wgs_84_crs = CRS("EPSG:4326")
>>> wgs_84_crs.datum.to_json_dict()
{'type': 'GeodeticReferenceFrame', 'name': 'World Geodetic System 1984', 'ellipsoid': {'name': 'WGS 84', 'semi_major_axis': 6378137, 'inverse_flattening': 298.257223563}, 'area': 'World', 'bbox': {'south_latitude': -90, 'west_longitude': -180, 'north_latitude': 90, 'east_longitude': 180}, 'id': {'authority': 'EPSG', 'code': 6326}}
>>> crs_dict["base_crs"]["datum"] = wgs_84_crs.datum.to_json_dict()
>>> final_crs = CRS(crs_dict)
>>> final_crs
<Projected CRS: {"type": "ProjectedCRS", "name": "unknown", "base_ ...>
Name: unknown
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- undefined
Coordinate Operation:
- name: unknown
- method: Modified Azimuthal Equidistant
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

[snowman2]

After some thinking, I might add a simpler interface like:

>>> from pyproj.crs import Ellipsoid, CRS
>>> wgs_84_crs = CRS("EPSG:4326")
>>> wgs_84_crs
<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

>>> ell = Ellipsoid.from_string("urn:ogc:def:ellipsoid:EPSG::7001")
>>> wgs_84_crs.from_self(ellipsoid=ell)
<Geographic 2D CRS: {"type": "GeographicCRS", "name": "WGS 84", "datum ...>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984
- Ellipsoid: Airy 1830
- Prime Meridian: Greenwich

[snowman2]

@djhoese, maybe something like this would be useful for you in transitioning from PROJ to the more verbose WKT/PROJ JSON?

[djhoese]

Yes, it looks like it. Accessing/using the information that I used to get from the PROJ dict will be the main transition for some of my packages. Thanks for letting me know.

[snowman2]

Sounds good 👍. I think being able to tweak parameters in the CRS would be a nice feature to have. But, with the added complexity, I haven't come up with a simpler interface yet. Maybe having preset CRS classes where you can modify the parameters similar to cartopy would be useful. 🤔

[djhoese]

Oh I guess I didn't realize this was about assigning new values to a CRS. Most of my users are used to things like:

from pyresample.geometry import AreaDefinition
area = AreaDefinition('my_area_name', {PROJ parameter dictionary}, rows, columns, area_extent)

Or writing the YAML equivalent of the above. Creating areas from other areas isn't done a lot by me, but I'm sure users would find it useful.

As a library other being able to do my_crs.datum and the equivalent for a lot of the other "high-level" parameters (proj, datum, ellps, maybe reference longitude, etc) would be nice.

[snowman2]

As a library other being able to do my_crs.datum and the equivalent for a lot of the other "high-level" parameters (proj, datum, ellps, maybe reference longitude, etc) would be nice.

For proj, that one is a bit trickier as the name could be in the coordinate_operation or somewhere else depending on the type of CRS.

The datum & ellps have the name properties that could be useful. Not a 1-1 match to the original PROJ string, but something you could use.

>>> from pyproj import CRS
>>> crs = CRS("ESRI:102719")
>>> crs.datum.name
'North American Datum 1983'
>>> crs.ellipsoid.name
'GRS 1980'
>>> crs.coordinate_operation.name
'NAD_1983_StatePlane_North_Carolina_FIPS_3200_Feet'
>>> crs.coordinate_operation.method_name
'Lambert Conic Conformal (2SP)'

For the reference longitude access, I would recommend looking into the params property of the coordinate_operation property on the CRS. See: https://gis.stackexchange.com/a/325290/144357

[snowman2]

Addressed in #509

[zackw]

This was supposed to have been fixed by #509, but it still doesn't work as of 3.1.0:

>>> from pyproj.crs import CRS, Datum
>>> CRS(proj="aeqd", datum="WGS84", lon_0=123, lat_0=1.3).to_wkt()
'PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ID["EPSG",6326]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",METHOD["Modified Azimuthal Equidistant",ID["EPSG",9832]],PARAMETER["Latitude of natural origin",1.3,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",123,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["False easting",0,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]'

>>> CRS(proj="aeqd", datum=Datum.from_name("WGS84"), lon_0=123, lat_0=1.3).to_wkt()
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/.../python3.9/site-packages/pyproj/crs/crs.py", line 313, in __init__
    self._local.crs = _CRS(self.srs)
  File "pyproj/_crs.pyx", line 2326, in pyproj._crs._CRS.__init__
pyproj.exceptions.CRSError: Invalid projection: +proj=aeqd +datum=World Geodetic System 1984 ensemble +lon_0=123 +lat_0=1.3 +type=crs: (Internal Proj Error: proj_create: Error 1027 (Invalid value for an argument): Unknown value for datum)

I see that there are now a bunch of new classes that let me do something almost equivalent via e.g.

>>> ProjectedCRS(AzumuthalEquidistantConversion(
...        longitude_natural_origin=123,
...        latitude_natural_origin=1.3),
...    geodetic_crs=CRS.from_epsg(4326)).to_wkt()
'PROJCRS["undefined",BASEGEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 (Transit)"],MEMBER["World Geodetic System 1984 (G730)"],MEMBER["World Geodetic System 1984 (G873)"],MEMBER["World Geodetic System 1984 (G1150)"],MEMBER["World Geodetic System 1984 (G1674)"],MEMBER["World Geodetic System 1984 (G1762)"],ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]],CONVERSION["unknown",METHOD["Modified Azimuthal Equidistant",ID["EPSG",9832]],PARAMETER["Latitude of natural origin",1.3,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",123,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["False easting",0,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]'

but I don't see why passing a Datum object as the datum= parameter to the CRS constructor shouldn't also work, as originally requested.

[snowman2]

I don't see why passing a Datum object as the datum= parameter to the CRS constructor shouldn't also work, as originally requested.

The datum= is for building PROJ strings and not all datums convert to PROJ string format.

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

PROJ strings have the potential to lose much of the information about a coordinate reference system (CRS).

More information: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems

The current method uses PROJ JSON, which is a better option for representing the CRS.

[zackw]

Why does it have to be only for building PROJ strings? That behavior makes sense when the value is a string, but when it's Datum object already, wouldn't it make more sense to bypass those conversions? Especially since they're lossy.

[snowman2]

@zackw do you have an example of what you are thinking that doesn't use a PROJ string?

[zackw]

The concrete thing I want to be able to do is construct arbitrary CRSes with the same datum as an existing CRS. In the particular program I am working on today, the existing CRS is the coordinate system of a map loaded by fiona.

[snowman2]

That sounds like a different problem from this issue. I would recommend creating a new issue to discuss.