Unidata
pyproj output not acceptd by metpy
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 --versionPython 3.9.2python -c 'import metpy; print(metpy.__version__)'1.0.1python -c 'import pyproj; print(pyproj.__version__)'3.0.1
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.
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.