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.

[Mikejmnez]

Just to add to this discussion, and to make it clear, there is no bug when making this type of plots. This has to do with "transforming" the square-cube geometry onto the arctic-crown geometry that emulates lat-lon grid. Making horizontal section in index space (use_coords=False)always works and, in most cases, making horizontal section with use_coords=True works for latitudes south of 60ish where the native grid has little distortion (in ECCO). This also applied when plotting LLC4320 data.

The particular problem/issue has to do with the presence of NaNs in the coordinates XC, YC(or XU, YU etc) used as x, y axes in the horizontal section. This issue DOES NOT have to do with the presence of nans on the data field being plotted (defined by varName), even though in this particular case, the presence of nans in the coords XC, YC is related to the presence of nans in the data variables (like Depth). That is, you can have nans in the variable varName and successfully make plots using use_coords=True, as long as there are no nans in the XC, YC fields.

One potential way to work around this is to "map" all nans in the coordinates (XC, YC) to their nearest non-nan neighbor or close to it, within machine precision, when making use of horizontal_section. So, there is no need to modify the dataset or anything. That way we collapse all the nans that arise from the transformation to an artificial near singularity in the horizontal section plot.