pyproj output not acceptd by metpy

[yt87]

Background: I want to add projection attributes to netcdf archive files. The projection specification is non-standard ( AWIPS/GFE). I thought pyproj would allow me to do it programatically, so far I found a couple of issues.
I am not sure thoose are bugs in metpy, one could argue that pyproj should be fixed instead.

import numpy as np
import pyproj
import xarray as xr
import metpy

proj_parms = {
    'proj': 'stere',
    'lat_0': 90.0,
    'lat_ts': 60.0,
    'lon_0': -150.0,
    'R': 6371229.0,
}
crs = pyproj.CRS.from_dict(proj_parms)
cf = crs.to_cf()
cf
crs2 = pyproj.CRS.from_cf(cf)
assert crs2.to_proj4() == crs.to_proj4()
# A trivial grid 
data = np.zeros((2, 2))
y = [3000, -2000]
x = [-1000, 1000]
dims = ['y', 'x']
da = xr.DataArray(data, coords=[y, x], dims=['y', 'x']).metpy.assign_crs(cf)
#cf2 = cf.copy()
#cf2['earth_radius'] = cf['semi_major_axis']
#da = xr.DataArray(data, coords=[y, x], dims=['y', 'x']).metpy.assign_crs(cf2)
crs3 = da.metpy.cartopy_crs
crs3.proj4_params

{'crs_wkt': 'PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",ELLIPSOID["unknown",6371229,0,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",METHOD["Polar Stereographic (variant B)",ID["EPSG",9829]],PARAMETER["Latitude of standard parallel",60,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8832]],PARAMETER["Longitude of origin",-150,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8833]],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)",south,MERIDIAN[90,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",south,MERIDIAN[180,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]',
 'semi_major_axis': 6371229.0,
 'semi_minor_axis': 6371229.0,
 'inverse_flattening': 0.0,
 'reference_ellipsoid_name': 'unknown',
 'longitude_of_prime_meridian': 0.0,
 'prime_meridian_name': 'Greenwich',
 'geographic_crs_name': 'unknown',
 'horizontal_datum_name': 'unknown',
 'projected_crs_name': 'unknown',
 'grid_mapping_name': 'polar_stereographic',
 'standard_parallel': 60.0,
 'straight_vertical_longitude_from_pole': -150.0,
 'false_easting': 0.0,
 'false_northing': 0.0}

<ipython-input-2-1322490a84cc>:17: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems
  assert crs2.to_proj4() == crs.to_proj4()

---------------------------------------------------------------------------
Proj4Error                                Traceback (most recent call last)
<ipython-input-2-1322490a84cc> in <module>
     24 #cf2['earth_radius'] = cf['semi_major_axis']
     25 #da = xr.DataArray(data, coords=[y, x], dims=['y', 'x']).metpy.assign_crs(cf2)
---> 26 crs3 = da.metpy.cartopy_crs
     27 crs3.proj4_params

~/miniconda3/envs/cf/lib/python3.9/site-packages/metpy/xarray.py in cartopy_crs(self)
    236     def cartopy_crs(self):
    237         """Return the coordinate reference system (CRS) as a cartopy object."""
--> 238         return self.crs.to_cartopy()
    239 
    240     @property

~/miniconda3/envs/cf/lib/python3.9/site-packages/metpy/plots/mapping.py in to_cartopy(self)
     76             raise ValueError(f'Unhandled projection: {proj_name}') from None
     77 
---> 78         return proj_handler(self._attrs, globe)
     79 
     80     def to_pyproj(self):

~/miniconda3/envs/cf/lib/python3.9/site-packages/metpy/plots/mapping.py in make_polar_stereo(attrs_dict, globe)
    195     kwargs = CFProjection.build_projection_kwargs(attrs_dict, attr_mapping)
    196 
--> 197     return ccrs.Stereographic(globe=globe, **kwargs)

~/miniconda3/envs/cf/lib/python3.9/site-packages/cartopy/crs.py in __init__(self, central_latitude, central_longitude, false_easting, false_northing, true_scale_latitude, scale_factor, globe)
   1455                 proj4_params.append(('k_0', scale_factor))
   1456 
-> 1457         super().__init__(proj4_params, globe=globe)
   1458 
   1459         # TODO: Let the globe return the semimajor axis always.

lib/cartopy/_crs.pyx in cartopy._crs.CRS.__init__()

Proj4Error: Error from proj: b'reciprocal flattening (1/f) = 0'

• Problem description: this should explain why the current behavior is a problem and why
  the expected output is a better solution.

I would expect pyproj CF output to be accepted by metpy.assign_crs(). The relevant logic in in /metpy/plot/mapping.py, property cartopy_globe. I don't understand the comment there, WGS84 definition have specific minor and major axes, with different size (I would force ellipsoid to be a sphere when the given axes are the same).

This can be overcome by adding earth_radius to the input dictionary, as shown by the commented out lines. However that leads to another issue:crs3.proj4_params evaluate to

{'ellps': 'sphere',
 'a': 6371229.0,
 'b': 6371229.0,
 'proj': 'stere',
 'lat_0': 0.0,
 'lon_0': -150.0,
 'x_0': 0.0,
 'y_0': 0.0,
 'lat_ts': 60.0}

The value of lat_0 is lost. This is caused by missing latitude of projection centre in the CF listing. It seems that pyproj determines tis parameter based on lat_ts - that's why the assertion crs2.to_proj4() == crs.to_proj4() succeedes.

• Which platform (Linux, Windows, Mac, etc.)

Fedora 33

• Versions. Include the output of:
  • python --version
Python 3.9.2
  • python -c 'import metpy; print(metpy.__version__)'
1.0.1
  • python -c 'import pyproj; print(pyproj.__version__)'
3.0.1

[dopplershift]

This definitely looks like a bug in the metadata mapping.

[dopplershift]

Thanks for the detailed run down and example code. That error is because the conversion from PyProj to CF is giving us a value of 0 for inverse_flattening. While CF says this is legal for a sphere, PROJ gets angry when we pass it this value (through CartoPy). The solution is to update our mapping to interpret the 0 inverse_flattening as a spherical datum and to not pass that value on.

The problem with the lat_0 being dropped is a problem where the conversion from PROJ to CF does not include a value for latitude_of_projection_origin (opened as a bug upstream pyproj4/pyproj#856)--according to the CF Document this value is required and we are relying on it when creating the projection.

[dopplershift]

I think there's still something to be done about the stereographic projection here, since it's not at all clear that Polar Stereographic (style-B) requires a latitude of the origin.