Pyproj return "None" when "towgs" exists

xycarto commented 1 year ago

Pyproj returns None when towgs exists

Code Sample, a copy-pastable example if possible

from pyproj import CRS
proj4 = "+proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
print(CRS.to_epsg(CRS.from_proj4(proj4)))

Returns None

Problem description

Ideally this should return 2193.

This appears to be common proj4 string for 'epsg:2193`.

gdalsrsinfo -o proj4 epsg:2193

I suspect this is an issue with the addition of the +towgs84=0,0,0,0,0,0,0 in the proj4 and not an issue specific to epsg:2193

Normally, I would not rely on the proj4 string for operations, however, in my use case, I am forced to use it

Expected Output

proj4 = "+proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +ellps=GRS80 +units=m +no_defs"
print(CRS.to_epsg(CRS.from_proj4(proj4)))
2193

Environment Information

pip3 show pyproj
name: pyproj
Version: 3.5.0
Summary: Python interface to PROJ (cartographic projections and coordinate transformations library)
Home-page: https://github.com/pyproj4/pyproj
Author: Jeff Whitaker
Author-email: jeffrey.s.whitaker@noaa.gov
License: MIT
Location: /usr/local/lib/python3.10/dist-packages
Requires: certifi
Required-by: geopandas

pip3 show geopandas
Name: geopandas
Version: 0.12.2
Summary: Geographic pandas extensions
Home-page: http://geopandas.org
Author: GeoPandas contributors
Author-email: kjordahl@alum.mit.edu
License: BSD
Location: /usr/local/lib/python3.10/dist-packages
Requires: fiona, packaging, pandas, pyproj, shapely

gdalinfo --version
GDAL 3.4.1, released 2021/12/27

Installation method

pip3 install geopandas --upgrade

Also tried:

pip3 install pyproj --upgrade --force-reinstall

snowman2 commented 1 year ago

The EPSG code is for a Projected CRS:

>>> from pyproj import CRS
>>> crs = CRS("EPSG:2193")
>>> crs
<Projected CRS: EPSG:2193>
Name: NZGD2000 / New Zealand Transverse Mercator 2000
Axis Info [cartesian]:
- N[north]: Northing (metre)
- E[east]: Easting (metre)
Area of Use:
- name: New Zealand - North Island, South Island, Stewart Island - onshore.
- bounds: (166.37, -47.33, 178.63, -34.1)
Coordinate Operation:
- name: New Zealand Transverse Mercator 2000
- method: Transverse Mercator
Datum: New Zealand Geodetic Datum 2000
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich

>>> crs.to_proj4()
'+proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +ellps=GRS80 +units=m +no_defs +type=crs'

When you add the towsg84 the type of CRS is a Bound CRS:

>>> from pyproj import CRS
>>> proj4 = "+proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
>>> crs = CRS(proj4)
>>> crs
<Bound CRS: +proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=160 ...>
Name: unknown
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- undefined
Coordinate Operation:
- name: Transformation from unknown to WGS84
- method: Position Vector transformation (geog2D domain)
Datum: Unknown based on GRS80 ellipsoid using towgs84=0,0,0,0,0,0,0
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich
Source CRS: unknown

To get the EPSG code of the projected CRS, you use the source_crs attribute (see to_epsg docstring):

>>> crs = CRS(proj4)
>>> crs.source_crs.list_authority(min_confidence=0)
[AuthorityMatchInfo(auth_name='EPSG', code='2193', confidence=50)]
>>> crs.source_crs.to_epsg(min_confidence=50)
2193

xycarto commented 1 year ago

Thank you for the very clear explanation. Very helpful! Implemented into my process

pyproj4 / pyproj