Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Gridliner inline label projection hard-coded for WGS84 #1902

Closed
Bob131 opened this issue Oct 6, 2021 · 3 comments · Fixed by #1903
Closed

Gridliner inline label projection hard-coded for WGS84 #1902

Bob131 opened this issue Oct 6, 2021 · 3 comments · Fixed by #1903

Comments

@Bob131
Copy link

Bob131 commented Oct 6, 2021

Description

Calling cartopy.crs.Projection.gridlines() with a CRS using a non-WGS84 Globe and with x_inline=True raises an exception during drawing, apparently due to gridliner.py:800: when either axis is specified as inline, the transform used comes from a newly-instantiated PlateCarree projection with a default Globe (i.e., a WGS84 globe).

The raised exception message is ProjError: Error creating Transformer from CRS.: (Internal Proj Error: proj_create_operations: Source and target ellipsoid do not belong to the same celestial body).

Code to reproduce

import matplotlib.pyplot as plt
import cartopy.crs as ccrs

globe = ccrs.Globe(semimajor_axis=1, semiminor_axis=1, ellipse=None)
projection = data_ccrs = ccrs.PlateCarree(globe=globe)

plt.figure()
ax = plt.subplot(111, projection=projection)
ax.set_global()
ax.gridlines(crs=data_ccrs, draw_labels=['x'], x_inline=True)
plt.show()

Traceback

Traceback (most recent call last):
  File "/usr/lib64/python3.10/site-packages/matplotlib/backends/backend_gtk3cairo.py", line 31, in on_draw_event
    self.figure.draw(self._renderer)
  File "/usr/lib64/python3.10/site-packages/matplotlib/artist.py", line 72, in draw_wrapper
    result = draw(artist, renderer, *args, **kwargs)
  File "/usr/lib64/python3.10/site-packages/matplotlib/artist.py", line 49, in draw_wrapper
    return draw(artist, renderer)
  File "/usr/lib64/python3.10/site-packages/matplotlib/figure.py", line 2799, in draw
    mimage._draw_list_compositing_images(
  File "/usr/lib64/python3.10/site-packages/matplotlib/image.py", line 132, in _draw_list_compositing_images
    a.draw(renderer)
  File "/usr/lib64/python3.10/site-packages/matplotlib/artist.py", line 49, in draw_wrapper
    return draw(artist, renderer)
  File "/usr/lib64/python3.10/site-packages/cartopy/mpl/geoaxes.py", line 543, in draw
    self._draw_preprocess(renderer)
  File "/usr/lib64/python3.10/site-packages/cartopy/mpl/geoaxes.py", line 517, in _draw_preprocess
    gl._draw_gridliner(renderer=renderer)
  File "/usr/lib64/python3.10/site-packages/cartopy/mpl/gridliner.py", line 948, in _draw_gridliner
    this_path = update_artist(artist, renderer)
  File "/usr/lib64/python3.10/site-packages/cartopy/mpl/gridliner.py", line 775, in update_artist
    artist.update_bbox_position_size(renderer)
  File "/usr/lib64/python3.10/site-packages/matplotlib/text.py", line 506, in update_bbox_position_size
    posx, posy = self.get_transform().transform((posx, posy))
  File "/usr/lib64/python3.10/site-packages/matplotlib/transforms.py", line 1463, in transform
    res = self.transform_affine(self.transform_non_affine(values))
  File "/usr/lib64/python3.10/site-packages/matplotlib/transforms.py", line 2386, in transform_non_affine
    return self._a.transform_non_affine(points)
  File "/usr/lib64/python3.10/site-packages/cartopy/mpl/geoaxes.py", line 134, in transform_non_affine
    return prj.transform_points(self.source_projection,
  File "/usr/lib64/python3.10/site-packages/cartopy/crs.py", line 401, in transform_points
    _safe_pj_transform(src_crs, self, x, y, z, trap=trap)
  File "/usr/lib64/python3.10/site-packages/cartopy/crs.py", line 43, in _safe_pj_transform
    transformer = Transformer.from_crs(src_crs, tgt_crs, always_xy=True)
  File "/usr/lib64/python3.10/site-packages/pyproj/transformer.py", line 534, in from_crs
    return Transformer(
  File "/usr/lib64/python3.10/site-packages/pyproj/transformer.py", line 310, in __init__
    self._local.transformer = transformer_maker()
  File "/usr/lib64/python3.10/site-packages/pyproj/transformer.py", line 97, in __call__
    return _Transformer.from_crs(
  File "pyproj/_transformer.pyx", line 643, in pyproj._transformer._Transformer.from_crs
pyproj.exceptions.ProjError: Error creating Transformer from CRS.: (Internal Proj Error: proj_create_operations: Source and target ellipsoid do not belong to the same celestial body)
Full environment definition

Operating system

Fedora 35

Cartopy version

v0.20.0, but appears to be an issue on master.

@greglucas
Copy link
Contributor

That makes it seem like all inline grid labels would have to be in lat/lon coordinates currently. I'm not entirely clear why the crs transform the user inputs with the gridliner isn't used for that. A PR to address this would be welcome!

@QuLogic
Copy link
Member

QuLogic commented Oct 7, 2021

We obviously don't test any inconsistent globes very much. Changing to the CRS doesn't break any test at all, but seems the right thing to do.

@Bob131
Copy link
Author

Bob131 commented Oct 10, 2021

Thanks for getting this fixed!

Sorry for the radio silence; I noticed I had made a mistake in my description of the problem and wanted to get to the bottom of it. It turns out that it's not so much that proj objects to transformations with different geodetic systems or ellipsoids (the same reproducer appears to work fine specifying an Earth ellipsoid known to proj) but actually due to a heuristic it uses to classify ellipsoids.

The error is thrown after checking equality of celestialBody of two CRSs: coordinateoperationfactory.cpp:3006.

The following traceback snippet shows how our newly-created ellipsoid gets its celestialBody:

#10 0x00007f43f9ea1e42 in osgeo::proj::io::PROJStringParser::Private::buildDatum(osgeo::proj::io::Step&, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&) (this=0x7f43f81ecef0, step=<optimized out>, title="") at /usr/src/debug/proj-8.1.1-1.fc35.x86_64/src/iso19111/io.cpp:9102
#11 0x00007f43f9ea598b in osgeo::proj::io::PROJStringParser::Private::buildGeographicCRS(int, int, int, bool) (this=0x7f43f81ecef0, iStep=0, iUnitConvert=-1, iAxisSwap=-1, ignorePROJAxis=<optimized out>) at /usr/src/debug/proj-8.1.1-1.fc35.x86_64/src/iso19111/io.cpp:9414
#12 0x00007f43f9eb0275 in osgeo::proj::io::PROJStringParser::createFromPROJString(std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&) (this=0x7fff18018700, projString=<optimized out>) at /usr/src/debug/proj-8.1.1-1.fc35.x86_64/src/iso19111/io.cpp:10511
#13 0x00007f43f9e9831c in osgeo::proj::io::createFromUserInput(std::string const&, osgeo::proj::io::DatabaseContextPtr const&, bool, PJ_CONTEXT*) (text=<optimized out>, dbContext=std::shared_ptr<osgeo::proj::io::DatabaseContext> (empty) = {...}, usePROJ4InitRules=<optimized out>, ctx=0x7f43f2b688d0) at /usr/src/debug/proj-8.1.1-1.fc35.x86_64/src/iso19111/io.cpp:6385
#14 0x00007f43f9e9b073 in osgeo::proj::io::createFromUserInput(std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, pj_ctx*) (text="+proj=eqc +a=1 +b=1 +lon_0=0.0 +to_meter=0.0174532925199433 +vto_meter=1 +no_defs +type=crs", ctx=0x7f43f2b688d0) at /usr/src/debug/proj-8.1.1-1.fc35.x86_64/src/iso19111/io.cpp:6972

The relevant line of buildDatum is io.cpp:9302, calling Ellipsoid::guessBodyName passing the length of the semi-major axis.

I realise this information probably isn't very useful for the (now-solved) problem at hand, but it seemed like a good idea to document it here for future reference.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging a pull request may close this issue.

3 participants