BUG: incompatibility with scipy 1.11.0

[neutrinoceros]

Description

Scipy 1.11.0 (released yesterday) appear to break use cases exercised in yt's CI

More specifically the breaking change seems to be scipy/scipy#18502

Code to reproduce

import numpy as np
import cartopy.crs as ccrs
from matplotlib.figure import Figure

fig = Figure()
ax = fig.add_axes((0, 0, 1, 1), projection=ccrs.Mollweide())
ax.imshow([[0, 1], [0, 1]], transform=ccrs.PlateCarree())

Traceback

Traceback (most recent call last):
  File "/private/tmp/t.py", line 7, in <module>
    ax.imshow([[0, 1], [0, 1]], transform=ccrs.PlateCarree())
  File "/private/tmp/.venv/lib/python3.11/site-packages/cartopy/mpl/geoaxes.py", line 318, in wrapper
    return func(self, *args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/private/tmp/.venv/lib/python3.11/site-packages/cartopy/mpl/geoaxes.py", line 1331, in imshow
    img, extent = warp_array(img,
                  ^^^^^^^^^^^^^^^
  File "/private/tmp/.venv/lib/python3.11/site-packages/cartopy/img_transform.py", line 192, in warp_array
    array = regrid(array, source_native_xy[0], source_native_xy[1],
            ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/private/tmp/.venv/lib/python3.11/site-packages/cartopy/img_transform.py", line 278, in regrid
    _, indices = kdtree.query(target_xyz, k=1)
                 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "_ckdtree.pyx", line 795, in scipy.spatial._ckdtree.cKDTree.query
ValueError: 'x' must be finite, check for nan or inf values

Full environment definition

Operating system

Discovered in CI on Ubuntu, and reproduced on macOS

Cartopy version

0.21.1

pip list

Package         Version
--------------- --------
Cartopy         0.21.1
certifi         2023.5.7
contourpy       1.1.0
cycler          0.11.0
fonttools       4.40.0
kiwisolver      1.4.4
matplotlib      3.7.1
numpy           1.25.0
packaging       23.1
Pillow          9.5.0
pip             23.1.2
pyparsing       3.1.0
pyproj          3.6.0
pyshp           2.3.1
python-dateutil 2.8.2
scipy           1.11.0
setuptools      65.5.0
shapely         2.0.1
six             1.16.0

[chrishavlin]

I was able to reproduce this error with the cartopy example for plotting images (link) for a number of projection choices: Mollweide, Orthographic, LambertAzimuthalEqualArea, LambertCylindrical all result in the same error as above (simply changing the ax = plt.axes(projection=ccrs.PlateCarree()) line to use different projections).

[chrishavlin]

So I think the issue is simply that previously cartopy didn't worry if the proj4 projection returned invalid target_xyz values and relied on masking after any kdtree calls. Here's a diff that attempts to fix the error in this PR by filtering the target_xyz values before evaluating with the kdtree:

diff --git a/lib/cartopy/img_transform.py b/lib/cartopy/img_transform.py
index 14af226e..a1d8b726 100644
--- a/lib/cartopy/img_transform.py
+++ b/lib/cartopy/img_transform.py
@@ -268,6 +268,7 @@ def regrid(array, source_x_coords, source_y_coords, source_proj, target_proj,
                                               target_x_points.flatten(),
                                               target_y_points.flatten())
 
     if _is_pykdtree:
         kdtree = pykdtree.kdtree.KDTree(xyz)
         # Use sqr_dists=True because we don't care about distances,
@@ -277,7 +278,11 @@ def regrid(array, source_x_coords, source_y_coords, source_proj, target_proj,
         # Versions of scipy >= v0.16 added the balanced_tree argument,
         # which caused the KDTree to hang with this input.
         kdtree = scipy.spatial.cKDTree(xyz, balanced_tree=False)
-        _, indices = kdtree.query(target_xyz, k=1)
+
+        # precompute finite mask to avoid scipy kdtree error
+        target_finite_mask = np.array([np.isfinite(target_xyz[:, i]) for i in range(3)]).all(axis=0)
+        indices = np.zeros(target_xyz.shape[0], dtype=int)
+        _, indices[target_finite_mask] = kdtree.query(target_xyz[target_finite_mask, :], k=1)
     mask = indices >= len(xyz)
     indices[mask] = 0

In my limited testing that seems to fix the problem, but I'm not familiar enough with the cartopy codebase to feel comfortable submitting this as PR. Still, thought it might be helpful for you all.

[greglucas]

@chrishavlin there is only one way to get comfortable and that is to contribute :) Your PR would be more than welcome!

Seems like a pretty reasonable fix, and may even be appropriate to put above that block and use it regardless of whether we are using scipy or pykdtree. We might need to update the mask below too? I think you can probably do an np.any(..., axis=-1) or something similar rather than iterating over indices too.

[chrishavlin]

well I was worried about having the bandwidth in the short term to work out the related mask refactoring that you mention. I can certainly submit a PR based on the above diff to get it started, I just might be a little slow to respond to reviews and thought this might be a critical enough issue that someone more familiar might want to get a fix in faster :) But I'll go ahead and submit it.

[CommonClimate]

+1 for a fix on this; it's affecting the Pyleoclim package, which has both scipy and cartopy as dependencies.

[GermainRV]

I was all weekend with the same error. I also think it is due to some incompatibilities of some packages. Basically, I'm transforming a geostationary projection image to PlateCarreee (for example), it was working fine until a few days ago (today is July 2nd), now I get this error (in the last line). ValueError: 'x' must be finite, check for nan or inf values . I have even used the example from cartopy documentation to verify if it was my code problem but I get the same error.
Here is the link to the example from cartopy documentation so you can try and see if it runs: Cartopy example: Reprojecting images from a Geostationary projection

Just in case, I also indicate the versions of some of my packages

``` (satimg) C:\Users\germa>conda list | findstr "numpy matplotlib scipy cartopy" cartopy 0.21.1 py311h178a126_1 conda-forge matplotlib 3.7.1 py311h1ea47a8_0 conda-forge matplotlib-base 3.7.1 py311h6e989c2_0 conda-forge matplotlib-inline 0.1.6 pyhd8ed1ab_0 conda-forge numpy 1.25.0 py311h0b4df5a_0 conda-forge scipy 1.11.1 py311h37ff6ca_0 conda-forge ```

Solution

If you need to run your code and this new scipy version is the problem, install:
conda install scipy=1.10.0

[chrishavlin]

To add to @GermainRV 's comment -- if there are reasons you need scipy>=1.11.0, you can also install pykdtree and cartopy will not use scipy for the kdtree (which is what is causing the error).