-
Notifications
You must be signed in to change notification settings - Fork 297
Closed
Labels
Description
🐛 Bug Report
Unfortunately there are still bugs related to #4408 in iris 3.2.0rc0, even with @wjbenfold's fix (#4497, thanks for that by the way). I tested the iris.analysis.cartography.rotate_winds()
function in the release candidate today and discovered that conversion between iris & cartopy projections / CRS continues to assume an Earth radius, failing for coordinate systems with semi-major axes different to iris.fileformats.pp.EARTH_RADIUS
.
How To Reproduce
Steps to reproduce the behaviour:
There are actually 2 seemingly related bugs, both of which trigger an error in a call of the transform_points()
method of a CRS.
- (mostly
iris
-side) Usingrotate_winds()
on cubes with a non-Earth radiusGeogCS
- see the error message below. It seems to be caused by this line in the cartography module. To reproduce:
import cartopy.crs as ccrs
from iris.coord_systems import GeogCS
import numpy as np
x, y = np.meshgrid([1,2,3], [4,5,6])
non_earth_cs = GeogCS(123456.7)
crs_latlon = ccrs.Geodetic(globe=ccrs.Globe(ellipse="sphere")) # Same as L930 in cartography.py
crs_latlon.transform_points(non_earth_cs.as_cartopy_projection(), x, y)
- (mostly
cartopy
-side) Just callingas_cartopy_projection()
on a non-Earth radius CS results in an error in Jupyter Lab when theProjection._repr_html_(self)
method is called, which in turn calls thetransform_points()
method in the same way as above.
from iris.coord_systems import GeogCS
non_earth_cs = GeogCS(123456.7)
non_earth_cs.as_cartopy_projection()
Environment
- OS & Version: e.g., Ubuntu 20.04 LTS
- Python Version: 3.10.2
- Iris Version: 3.2.0.rc0
- Cartopy Version: 0.20.2
- Pyproj Version: 3.3.0
Additional context
Error message
---------------------------------------------------------------------------
ProjError Traceback (most recent call last)
Input In [169], in <module>
----> 1 tl_u, tl_v = rotate_winds(u, v, tl_cs)
File /lustre/dirac3/home/dc-serg1/mambaforge/envs/exo/lib/python3.10/site-packages/iris/analysis/cartography.py:1179, in rotate_winds(u_cube, v_cube, target_cs)
1176 vt_cube.rename("transformed_{}".format(v_cube.name()))
1178 # Get distance scalings for source crs.
-> 1179 ds_dx1, ds_dy1 = _crs_distance_differentials(src_crs, x, y)
1181 # Get distance scalings for target crs.
1182 x2, y2 = _transform_xy(src_crs, x, y, target_crs)
File /lustre/dirac3/home/dc-serg1/mambaforge/envs/exo/lib/python3.10/site-packages/iris/analysis/cartography.py:932, in _crs_distance_differentials(crs, x, y)
930 crs_latlon = ccrs.Geodetic(globe=ccrs.Globe(ellipse="sphere"))
931 # Transform points to true-latlon (just to get the true latitudes).
--> 932 _, true_lat = _transform_xy(crs, x, y, crs_latlon)
933 # Get coordinate differentials w.r.t. true-latlon.
934 dlon_dx, dlat_dx, dlon_dy, dlat_dy = _inter_crs_differentials(
935 crs, x, y, crs_latlon
936 )
File /lustre/dirac3/home/dc-serg1/mambaforge/envs/exo/lib/python3.10/site-packages/iris/analysis/cartography.py:854, in _transform_xy(crs_from, x, y, crs_to)
838 def _transform_xy(crs_from, x, y, crs_to):
839 """
840 Shorthand function to transform 2d points between coordinate
841 reference systems.
(...)
852
853 """
--> 854 pts = crs_to.transform_points(crs_from, x, y)
855 return pts[..., 0], pts[..., 1]
File /lustre/dirac3/home/dc-serg1/mambaforge/envs/exo/lib/python3.10/site-packages/cartopy/crs.py:408, in CRS.transform_points(self, src_crs, x, y, z, trap)
405 x[to_180] = (((x[to_180] + 180) % 360) - 180)
406 try:
407 result[:, 0], result[:, 1], result[:, 2] = \
--> 408 _safe_pj_transform(src_crs, self, x, y, z, trap=trap)
409 except ProjError as err:
410 msg = str(err).lower()
File /lustre/dirac3/home/dc-serg1/mambaforge/envs/exo/lib/python3.10/site-packages/cartopy/crs.py:50, in _safe_pj_transform(src_crs, tgt_crs, x, y, z, trap)
49 def _safe_pj_transform(src_crs, tgt_crs, x, y, z=None, trap=True):
---> 50 transformer = _get_transformer_from_crs(src_crs, tgt_crs)
51 transformed_coords = transformer.transform(x, y, z, errcheck=trap)
52 if z is None:
File /lustre/dirac3/home/dc-serg1/mambaforge/envs/exo/lib/python3.10/site-packages/cartopy/crs.py:46, in _get_transformer_from_crs(src_crs, tgt_crs)
44 @lru_cache()
45 def _get_transformer_from_crs(src_crs, tgt_crs):
---> 46 return Transformer.from_crs(src_crs, tgt_crs, always_xy=True)
File /lustre/dirac3/home/dc-serg1/mambaforge/envs/exo/lib/python3.10/site-packages/pyproj/transformer.py:565, in Transformer.from_crs(crs_from, crs_to, skip_equivalent, always_xy, area_of_interest, authority, accuracy, allow_ballpark)
558 if skip_equivalent:
559 warnings.warn(
560 "skip_equivalent is deprecated.",
561 DeprecationWarning,
562 stacklevel=2,
563 )
--> 565 return Transformer(
566 TransformerFromCRS(
567 cstrencode(CRS.from_user_input(crs_from).srs),
568 cstrencode(CRS.from_user_input(crs_to).srs),
569 always_xy=always_xy,
570 area_of_interest=area_of_interest,
571 authority=authority,
572 accuracy=accuracy if accuracy is None else str(accuracy),
573 allow_ballpark=allow_ballpark,
574 )
575 )
File /lustre/dirac3/home/dc-serg1/mambaforge/envs/exo/lib/python3.10/site-packages/pyproj/transformer.py:310, in Transformer.__init__(self, transformer_maker)
304 raise ProjError(
305 "Transformer must be initialized using: "
306 "'from_crs', 'from_pipeline', or 'from_proj'."
307 )
309 self._local = TransformerLocal()
--> 310 self._local.transformer = transformer_maker()
311 self._transformer_maker = transformer_maker
File /lustre/dirac3/home/dc-serg1/mambaforge/envs/exo/lib/python3.10/site-packages/pyproj/transformer.py:97, in TransformerFromCRS.__call__(self)
91 def __call__(self) -> _Transformer:
92 """
93 Returns
94 -------
95 _Transformer
96 """
---> 97 return _Transformer.from_crs(
98 self.crs_from,
99 self.crs_to,
100 always_xy=self.always_xy,
101 area_of_interest=self.area_of_interest,
102 authority=self.authority,
103 accuracy=self.accuracy,
104 allow_ballpark=self.allow_ballpark,
105 )
File pyproj/_transformer.pyx:1001, in pyproj._transformer._Transformer.from_crs()
ProjError: Error creating Transformer from CRS.: (Internal Proj Error: proj_create_operations: Source and target ellipsoid do not belong to the same celestial body)