potential bug: regression in cartopy 0.19 ?

[neutrinoceros]

Description

yt's test suite exercises a couple test cases using cartopy. Since August yt has been pinning cartopy to 0.18 to mitigate some compatibility issues while #1140 was being worked on.

Now that the original problem seems to be fixed on your side, I've allowed more recent cartopy versions to install, and discovered that one of our tests is now broken. See downstream issue yt-project/yt#3705

See bellow for a reproducer that doesn't use yt.

I'm very open to considering that we're doing something wrong here and that this script shouldn't work, which would imply that the bug is actually on the yt side, but given my limited understanding of it, I wanted to check in with cartopy devs first.

Code to reproduce

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
import cartopy

data = np.random.random_sample((10, 10))
transform = cartopy.crs.Miller()
projection = cartopy.crs.Mollweide()

fig = plt.figure()
ax = fig.add_axes((0, 0, 1, 1), projection=projection)
cax = fig.add_axes((1, 1, 1, 1))

im = ax.imshow(
    data,
    extent=[-180.0, 180.0, -90.0, 90.0],
    norm=LogNorm(),
    transform=transform,
)
cb = fig.colorbar(im, cax=cax)
fig.savefig("/tmp/demo.png")

Traceback

Traceback (most recent call last):
  File "/private/tmp/t_cartopy.py", line 20, in <module>
    cb = fig.colorbar(im, cax=cax)
  File "/Users/robcleme/.pyenv/versions/cartopy_tmp/lib/python3.10/site-packages/matplotlib/figure.py", line 1158, in colorbar
    cb = cbar.Colorbar(cax, mappable, **cb_kw)
  File "/Users/robcleme/.pyenv/versions/cartopy_tmp/lib/python3.10/site-packages/matplotlib/colorbar.py", line 482, in __init__
    self._reset_locator_formatter_scale()
  File "/Users/robcleme/.pyenv/versions/cartopy_tmp/lib/python3.10/site-packages/matplotlib/colorbar.py", line 1166, in _reset_locator_formatter_scale
    self._process_values()
  File "/Users/robcleme/.pyenv/versions/cartopy_tmp/lib/python3.10/site-packages/matplotlib/colorbar.py", line 1110, in _process_values
    b = self.norm.inverse(b)
  File "/Users/robcleme/.pyenv/versions/cartopy_tmp/lib/python3.10/site-packages/matplotlib/colors.py", line 1554, in inverse
    raise ValueError("Invalid vmin or vmax")
ValueError: Invalid vmin or vmax

Full environment definition

Operating system

I'm seeing this both on Windows (on CI) and MacOS

Cartopy version

0.20.1
I think I saw it too in 0.19.0.post1 but I've lost track of it and I can't easily install this version on my machine.

conda list

pip list

Package         Version
--------------- ---------
Cartopy         0.20.1
certifi         2021.10.8
cycler          0.11.0
fonttools       4.28.3
kiwisolver      1.3.2
matplotlib      3.5.0
numpy           1.21.4
packaging       21.3
Pillow          8.4.0
pip             21.2.3
pyparsing       3.0.6
pyproj          3.3.0
pyshp           2.1.3
python-dateutil 2.8.2
scipy           1.7.3
setuptools      57.4.0
setuptools-scm  6.3.2
Shapely         1.8.0
six             1.16.0
tomli           1.2.2

[dopplershift]

Does this work if you remove the CartoPy projections and use plain matplotlib?

[neutrinoceros]

yes.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm

data = np.random.random_sample((10, 10))

fig = plt.figure()
ax = fig.add_axes((0, 0, 1, 1))
cax = fig.add_axes((1, 1, 1, 1))

im = ax.imshow(
    data,
    extent=[-180.0, 180.0, -90.0, 90.0],
    norm=LogNorm(),
)
cb = fig.colorbar(im, cax=cax)
fig.savefig("/tmp/demo.png")

produces this image

[dopplershift]

I'd say definitely a bug in cartopy then.

[greglucas]

My guess is that this is because the reprojection of the image fails and produces an all-masked result. If this is only a test, then you probably want to change Miller() to PlateCarree() for the transform. The extents of PlateCarree are (-180, 180), the extents of Miller are (-2e7, 2e7).

[neutrinoceros]

I think there might be several things going on here

• the affected test worked in cartopy 0.18, but as I'm not able to install it locally, I can't see what the reproducer script would produce. From @greglucas' response I'd guess that the test should not pass, so maybe 0.20's behaviour is actually correct ?

If anyone could try out my script with an older version of cartopy that doesn't crash and post the image it would be helpful.

My guess is that this is because the reprojection of the image fails and produces an all-masked result.

Indeed, if we remove the norm=LogNorm() argument in imshow the script doesn't fail but produces an empty image.
The final error is from the norm failing to produce finite extremas from all negative (i.e. masked ?) values.

The extents of PlateCarree are (-180, 180), the extents of Miller are (-2e7, 2e7).

ok this part is probably a bug in yt, because the extent are not user-defined as in my minimal reproducer script. What's user-specified is a str "Miller", and the rest is internal.

[rcomer]

I tried the script with Cartopy 0.18 and Matplotlib 3.3.4:

[greglucas]

Unfortunately, I purposefully decided to change the projection coordinates on some projections because they were using non-sensical globes to make the numbers closer to pseudo-lat/lon (of which, Miller appears to be one of those).
commit: 40b7fc8
Newer versions of PROJ complain about conversion between globes that aren't the same, so this was to address that issue. I didn't see any way to do this without a straight break, so I am not sure it was a bug in YT, rather than a breaking change on our end.

Thanks, Ruth! I think that image will actually be different in v0.20 too due to another change in image interpolation that turns those "waves" into straight lines now. Removing the geocentric pivot in this commit: 017f989

So, indeed @neutrinoceros, you hit multiple breaking changes here, sorry about that!

[neutrinoceros]

That's okay ! I'd be happy to declare cartopy 0.20 as the new minimal compatible version.

I may need assistance to perform the necessary changes there however.
Would it be okay to ping any of you guys in a yt PR for review when I can manage a patch ?

[neutrinoceros]

@greglucas I confirm that the PlateCarre transform is still compatible with Mollweide projection, however, the test that lead me here seems to be clearly targeting Miller transforms
https://github.com/yt-project/yt/blob/84cb4b40d91136a63be67e653b0f8997174f641e/yt/visualization/tests/test_geo_projections.py#L131-L152

Actually we're imitating a user asking for a "Miller" transform, and apparently we set the Mollweide projection internally, whatever the transform.
https://github.com/yt-project/yt/blob/84cb4b40d91136a63be67e653b0f8997174f641e/yt/geometry/coordinates/geographic_coordinates.py#L375

Is there a more suited default I can associate with the Miller transform, and in general, is there any way to determine a compatible projection at runtime, given a transform ?

[greglucas]

I'm not exactly clear what the data_projection vs data_transform is in yt... Is data_projection, what the MPL Axes are in, or is there a secondary pivot/transform somewhere along the path?

I don't know where your original test was failing, because for the reproducer here I think the key is extent=[-180.0, 180.0, -90.0, 90.0], but I don't see that defined in the test. So, I guess my quick thought is that there is probably a hard-coded (-180, 180, -90, 90) somewhere in the test suite that isn't quite right because that extent should really be data_transform dependent. So, either drop the extent=... (which will then use the current extent of the axes I think), or put in the full extent for the Miller transform which can be obtained from a Miller Axes.

ax = plt.axes(..., transform=ccrs.Miller())
ax.set_global()
ax.get_extent()

[neutrinoceros]

I'm not exactly clear what the data_projection vs data_transform is in yt... Is data_projection, what the MPL Axes are in,

correct !

or is there a secondary pivot/transform somewhere along the path?

not that I can tell

I don't know where your original test was failing, because for the reproducer here I think the key is extent=[-180.0, 180.0, -90.0, 90.0], but I don't see that defined in the test

It's effectively what happens in our test but it's pretty hard to track down where these values are coming from. I won't go too much into the details but it's not from the test suite, it's definitely in the code base itself.
Note that we have a comment in the code base that seems very related to the issue being discussed:
https://github.com/yt-project/yt/blob/84cb4b40d91136a63be67e653b0f8997174f641e/yt/visualization/base_plot_types.py#L254-L272

Unfortunately the reason why the code is disabled is only vaguely documented:

This is currently commented out because it breaks in some instances.

But I think that there's a good chance that the patch we need has to be happen here (and possibly nowhere else).
That said, none of the following two solutions seems to work

attempt 1: determine extent dynamically as suggested

self.axes.set_global()
extent = self.axes.get_extent()

attempt 2: reactivating the commented code

self.axes.set_extent(extent)

The following error is logged a few dozen times

ERROR    shapely.geos:geos.py:252 IllegalArgumentException: Argument must be Polygonal or LinearRing

Both patches still lead to the original error ValueError: Invalid vmin or vmax

[greglucas]

That comment has a few interesting tidbits. The kdtree to produce inf's is essentially what I think is going on here, except now it is nan's/masked everywhere for the small domain instead of inf's, which is what is causing issues with your LogNorm scaling. Also, the comment in that code seems related to this issue here of when the extent gets updated: #1468 so we should look into that some more on Cartopy's side as well.

The extent should be defined in the transform's coordinates, so the coordinate system your data is originally in, and not the one you're going to (axes). Sorry for the confusion there. Perhaps setting extent=None is the easiest here to get the global extent.

[neutrinoceros]

Thank you so much for your insight ! For reference, I've proposed a patch on the yt side base on your reply here yt-project/yt#3708

[neutrinoceros]

The yt patch has been merged, I don't think there's anything left to do here.