Error in `plot.horizontal_section` for certain cutouts

[ThomasHaine]

The following code errors in plot.horizontal_section if the 3rd point is included (in the Arctic), but not if it's excluded:

import oceanspy as ospy
from poseidon_viewer import get_shapes
import xoak
import matplotlib.pyplot as plt

od = ospy.open_oceandataset.from_catalog('ECCO')

# From Poseidon-viewer:
ps = [{"type":"Point","coordinates":[-47.0337310394683,24.77353204723849]},{"type":"Point","coordinates":[-169.68462560083955,21.754828557562874]},{"type":"Point","coordinates":[-155.84669435179075,77.10081606690576]}]
# ps = [{"type":"Point","coordinates":[-47.0337310394683,24.77353204723849]},{"type":"Point","coordinates":[-169.68462560083955,21.754828557562874]}]
lons, lats = ospy.utils.viewer_to_range(ps)

# Plot station locations...
cut_od = od.subsample.cutout(XRange = [min(lons)-10,max(lons)+10],YRange = [min(lats)-5,max(lats)+5])
fig = plt.figure(figsize=(12, 5))
ax = cut_od.plot.horizontal_section(varName="Depth")

The error is:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[3], line 16
     14 cut_od = od.subsample.cutout(XRange = [min(lons)-10,max(lons)+10],YRange = [min(lats)-5,max(lats)+5])
     15 fig = plt.figure(figsize=(12, 5))
---> 16 ax = cut_od.plot.horizontal_section(varName="Depth")

File ~/mambaforge/envs/Oceanography/lib/python3.10/site-packages/oceanspy/plot.py:1133, in _plotMethods.horizontal_section(self, **kwargs)
   1131 @_functools.wraps(horizontal_section)
   1132 def horizontal_section(self, **kwargs):
-> 1133     return horizontal_section(self._od, **kwargs)

File ~/mambaforge/envs/Oceanography/lib/python3.10/site-packages/oceanspy/plot.py:722, in horizontal_section(od, varName, plotType, use_coords, contourName, meanAxes, intAxes, contour_kwargs, clabel_kwargs, cutout_kwargs, **kwargs)
    720 args = {"x": X_name, "y": Y_name, **kwargs}
    721 plotfunc = eval("_xr.plot." + plotType)
--> 722 p = plotfunc(da, **args)
    724 # Contour
    725 if contourName is not None:

File ~/mambaforge/envs/Oceanography/lib/python3.10/site-packages/xarray/plot/dataarray_plot.py:1597, in _plot2d.<locals>.newplotfunc(***failed resolving arguments***)
   1593     raise ValueError("plt.imshow's `aspect` kwarg is not available in xarray")
   1595 ax = get_axis(figsize, size, aspect, ax, **subplot_kws)
-> 1597 primitive = plotfunc(
   1598     xplt,
   1599     yplt,
   1600     zval,
   1601     ax=ax,
   1602     cmap=cmap_params["cmap"],
   1603     vmin=cmap_params["vmin"],
   1604     vmax=cmap_params["vmax"],
   1605     norm=cmap_params["norm"],
   1606     **kwargs,
   1607 )
   1609 # Label the plot with metadata
   1610 if add_labels:

File ~/mambaforge/envs/Oceanography/lib/python3.10/site-packages/xarray/plot/dataarray_plot.py:2320, in pcolormesh(x, y, z, ax, xscale, yscale, infer_intervals, **kwargs)
   2317         y = _infer_interval_breaks(y, axis=0, scale=yscale)
   2319 ax.grid(False)
-> 2320 primitive = ax.pcolormesh(x, y, z, **kwargs)
   2322 # by default, pcolormesh picks "round" values for bounds
   2323 # this results in ugly looking plots with lots of surrounding whitespace
   2324 if not hasattr(ax, "projection") and x.ndim == 1 and y.ndim == 1:
   2325     # not a cartopy geoaxis

File ~/mambaforge/envs/Oceanography/lib/python3.10/site-packages/matplotlib/__init__.py:1442, in _preprocess_data.<locals>.inner(ax, data, *args, **kwargs)
   1439 @functools.wraps(func)
   1440 def inner(ax, *args, data=None, **kwargs):
   1441     if data is None:
-> 1442         return func(ax, *map(sanitize_sequence, args), **kwargs)
   1444     bound = new_sig.bind(ax, *args, **kwargs)
   1445     auto_label = (bound.arguments.get(label_namer)
   1446                   or bound.kwargs.get(label_namer))

File ~/mambaforge/envs/Oceanography/lib/python3.10/site-packages/matplotlib/axes/_axes.py:6220, in Axes.pcolormesh(self, alpha, norm, cmap, vmin, vmax, shading, antialiased, *args, **kwargs)
   6217 shading = shading.lower()
   6218 kwargs.setdefault('edgecolors', 'none')
-> 6220 X, Y, C, shading = self._pcolorargs('pcolormesh', *args,
   6221                                     shading=shading, kwargs=kwargs)
   6222 coords = np.stack([X, Y], axis=-1)
   6223 # convert to one dimensional array, except for 3D RGB(A) arrays

File ~/mambaforge/envs/Oceanography/lib/python3.10/site-packages/matplotlib/axes/_axes.py:5717, in Axes._pcolorargs(self, funcname, shading, *args, **kwargs)
   5715 if funcname == 'pcolormesh':
   5716     if np.ma.is_masked(X) or np.ma.is_masked(Y):
-> 5717         raise ValueError(
   5718             'x and y arguments to pcolormesh cannot have '
   5719             'non-finite values or be of type '
   5720             'numpy.ma.core.MaskedArray with masked values')
   5721     # safe_masked_invalid() returns an ndarray for dtypes other
   5722     # than floating point.
   5723     if isinstance(X, np.ma.core.MaskedArray):

ValueError: x and y arguments to pcolormesh cannot have non-finite values or be of type numpy.ma.core.MaskedArray with masked values

[Mikejmnez]

the error means there are NaNs in the XC, YC variables used by cartopy when making plots. This occurs frequently when the cutout includes some arctic (face=6) data in the ECCO/LLC4320, where the grid greatly diverts from lat/lon. I found that this occurs even more frequently in ECCO than in LLC4320 (due to grid resolution and subtle differences in the respective grid generation).

So this is one of two scenarios when we don't want to use the model coordinates (XC, YC) when making horizontal maps via oceanspy. The other scenario is when making plots of cutouts localized to the Pacific Ocean involving data across the discontinuity at the +- 180 degree..

Solution

Set use_coords=False when making horizontal_sections. So in your example:

fig = plt.figure(figsize=(12, 5))
ax = cut_od.plot.horizontal_section(varName="Depth", use_coords=False)

that makes the plot in index space, instead of using the map projection.

[ThomasHaine]

Thanks @Mikejmnez! That solves the issue. Maybe convert this issue to discussion?

[Mikejmnez]

Sure! I am glad this issue came up so it can be documented. It is not clear from the error one gets what the error is (which not very helpful ha) and how to work around it. But yeah, no error here, no bug, just the product of making cutouts in these grids.