Updates to support wkt with new proj version

[snowman2]

Just wondering what the plans are for updating in connection with https://gdalbarn.com.

[snowman2]

If there is interest, I am working on an updated version to proj.4 version 6.0.0 here

[snowman2]

The functionality in the CRS class is based on other fantastic projects:

• https://github.com/opendatacube/datacube-core/blob/83bae20d2a2469a6417097168fd4ede37fd2abe5/datacube/utils/geometry/_base.py
• https://github.com/mapbox/rasterio/blob/c13f0943b95c0eaa36ff3f620bd91107aa67b381/rasterio/_crs.pyx

Current improvements:

>>> from pyproj import CRS
>>> geog = CRS({'init': 'EPSG:4326'})
>>> geog.to_string()
'+proj=pipeline +step +proj=longlat +ellps=WGS84 +step +proj=unitconvert +xy_in=rad +xy_out=rad'
>>> geog.to_wkt()
'GEOGCRS["WGS 84",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]],CS[ellipsoidal,2],AXIS["longitude",east,ORDER[1],ANGLEUNIT["radian",1,ID["EPSG",9101]]],AXIS["latitude",north,ORDER[2],ANGLEUNIT["radian",1,ID["EPSG",9101]]],USAGE[SCOPE["unknown"],AREA["World"],BBOX[-90,-180,90,180]]]'
>>> geog.to_epsg()
4326

Ellipsoid & area of use information

>>> from pyproj import CRS
>>> crs1 = CRS.from_epsg(4326)
>>> crs1.to_string(4)
'+proj=longlat +datum=WGS84 +no_defs'
>>> crs1.ellipsoid.inverse_flattening
298.257223563
>>> crs1.ellipsoid.semi_major_metre
6378137.0
>>> crs1.ellipsoid.semi_minor_metre
6356752.314245179
>>> crs1.area_of_use.bounds
(-180.0, -90.0, 180.0, 90.0)
>>> crs1.area_of_use.name
'World'

WKT String support:

>>> from pyproj import CRS
>>> projection_string = (
...     'PROJCS["USA_Contiguous_Albers_Equal_Area_Conic_USGS_version",'
...     'GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",'
...     'SPHEROID["GRS_1980",6378137.0,298.257222101]],'
...     'PRIMEM["Greenwich",0.0],'
...     'UNIT["Degree",0.0174532925199433]],'
...     'PROJECTION["Albers"],'
...     'PARAMETER["false_easting",0.0],'
...     'PARAMETER["false_northing",0.0],'
...     'PARAMETER["central_meridian",-96.0],'
...     'PARAMETER["standard_parallel_1",29.5],'
...     'PARAMETER["standard_parallel_2",45.5],'
...     'PARAMETER["latitude_of_origin",23.0],'
...     'UNIT["Meter",1.0],'
...     'VERTCS["NAVD_1988",'
...     'VDATUM["North_American_Vertical_Datum_1988"],'
...     'PARAMETER["Vertical_Shift",0.0],'
...     'PARAMETER["Direction",1.0],UNIT["Centimeter",0.01]]]')
>>>
>>> proj = CRS(projection_string)
>>> proj.to_string()
'+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +ellps=GRS80'
>>> proj.to_string(4)
'+proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs'
>>> proj.to_wkt()
'PROJCRS["USA_Contiguous_Albers_Equal_Area_Conic_USGS_version",BASEGEOGCRS["NAD83",DATUM["North American Datum 1983",ELLIPSOID["GRS 1980",6378137,298.257222101,LENGTHUNIT["metre",1]],ID["EPSG",6269]],PRIMEM["Greenwich",0,ANGLEUNIT["Degree",0.0174532925199433]]],CONVERSION["unnamed",METHOD["Albers Equal Area",ID["EPSG",9822]],PARAMETER["Latitude of false origin",23,ANGLEUNIT["Degree",0.0174532925199433],ID["EPSG",8821]],PARAMETER["Longitude of false origin",-96,ANGLEUNIT["Degree",0.0174532925199433],ID["EPSG",8822]],PARAMETER["Latitude of 1st standard parallel",29.5,ANGLEUNIT["Degree",0.0174532925199433],ID["EPSG",8823]],PARAMETER["Latitude of 2nd standard parallel",45.5,ANGLEUNIT["Degree",0.0174532925199433],ID["EPSG",8824]],PARAMETER["Easting at false origin",0,LENGTHUNIT["metre",1],ID["EPSG",8826]],PARAMETER["Northing at false origin",0,LENGTHUNIT["metre",1],ID["EPSG",8827]]],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]]]]'
>>> proj.to_wkt("WKT1_GDAL")
'PROJCS["USA_Contiguous_Albers_Equal_Area_Conic_USGS_version",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0],UNIT["Degree",0.0174532925199433]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["latitude_of_center",23],PARAMETER["longitude_of_center",-96],PARAMETER["standard_parallel_1",29.5],PARAMETER["standard_parallel_2",45.5],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]'

>>> from pyproj import CRS
>>> eng_wkt = 'ENGCRS["A construction site CRS",\nEDATUM["P1",ANCHOR["Peg in south corner"]],\nCS[Cartesian,2],\nAXIS["site east",southWest,ORDER[1]],\nAXIS["site north",southEast,ORDER[2]],\nLENGTHUNIT["metre",1.0],\nTIMEEXTENT["date/time t1","date/time t2"]]'
>>> eng = CRS(eng_wkt)
proj_obj_get_ellipsoid: CRS has no geodetic CRS
proj_obj_get_ellipsoid: Object is not a CRS or GeodeticReferenceFrame
>>> eng.to_wkt()
'ENGCRS["A construction site CRS",EDATUM["P1",ANCHOR["Peg in south corner"]],CS[Cartesian,2],AXIS["site east",southWest,ORDER[1],LENGTHUNIT["metre",1]],AXIS["site north",southEast,ORDER[2],LENGTHUNIT["metre",1]],USAGE[SCOPE["unknown"]]]'
>>> eng.to_wkt("WKT1_GDAL")
'LOCAL_CS["A construction site CRS",LOCAL_DATUM["P1",32767]]'

The Proj object now also supports initialization from various methods:

>>> from pyproj import Proj, CRS
>>> crs = CRS({'init': 'EPSG:4326'})
>>> proj1 = Proj(crs)
>>> proj2 = Proj(crs.to_string())
>>> proj3 = Proj(crs.to_wkt())
>>> proj4 = Proj(4326)

[jswhit]

@snowman2, this is very nice! Would you like to merge the changes in https://github.com/snowman2/pyproj/tree/gdalbarn here? If so, a PR would be welcome.

[snowman2]

More updates:

>>> from pyproj import CRS
>>> crs = CRS('+init=epsg:32667')
>>> crs
<CRS: +init=epsg:32667>
Ellipsoid:
- semi_major_metre: 6378137.00
- semi_minor_metre: 6356752.31
- inverse_flattening: 298.26
Area of Use:
- name: USA - GoM OCS - east of 84°W
- bounds: (-84.09, 23.82, -81.17, 29.94)
Prime Meridian:
- longitude: 0.0000
- unit_name: degree
- unit_conversion_factor: 0.01745329
Axis Info:
- Easting[E] (east) EPSG:9003 (US survey foot)
- Northing[N] (north) EPSG:9003 (US survey foot)

>>> crs.get_geod()
Geod('+a=6378137 +f=0.0033528106647475126')