Skip to content

Coordinate systems with non-Earth radii do not work properly in iris 3.2.rc0 #4582

@dennissergeev

Description

@dennissergeev

🐛 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.

  1. (mostly iris-side) Using rotate_winds() on cubes with a non-Earth radius GeogCS - 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)
  1. (mostly cartopy-side) Just calling as_cartopy_projection() on a non-Earth radius CS results in an error in Jupyter Lab when the Projection._repr_html_(self) method is called, which in turn calls the transform_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)

Metadata

Metadata

Assignees

No one assigned

    Labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions