Confusing behavior when plotting size==0 array, from dataset with ndim>1

[Sevans711]

Version

v2026.6.1

How did you install UXarray?

Source

What happened?

Attempting to plot a DataArray with size=0 led to confusing behavior, but only when it is connected(?) to an array with ndim>1. Noticed this while working on #1285, where the array comes from subset.bounding_box with a region that does not overlap anywhere with the grid.

Example:

ds = ...  # v1.nc from meshfiles/ugrid/geoflow-small; this includes 'time' and 'meshLayers' dims.
arr0 = ds['v1'].isel(meshLayers=0, time=0)
# pick subset where there is no data:
arr1 = arr0.subset.bounding_box(lon_bounds=(-1000, -999), lat_bounds=(-1,1))
print(arr1.size)   # prints 0
print(arr1.ndim)   # prints 1
print(arr1.dims)  # prints (n_node,)
print(arr1)   # this works just fine, it shows there are 0 values.
arr1.plot()

This crashes with ValueError: Per-column arrays must each be 1-dimensional, which is the error that shows up when trying to plot a 2D+ array. That's confusing; this array is only 1D. Even more confusing is the following. After running that, I ran:

print(arr1.size)  # this is still 0
print(arr1.ndim)   # this is now 3. How did that happen?
print(arr1)   # this now crashes, with ValueError: zip() argument 2 is longer than argument 1

Additional notes:

• Using the same example above but only .isel(time=0) leads to printout of arr1.ndim==2, arr1.dims=('meshLayers', 'n_node') initially, but then still the same crash, and the same arr1.ndim==3 after the crash.
• Using an example that initially has only 1D instead, e.g. the files from meshfiles/ugrid/quad-hexagon, actually creates a plot without complaint (the plot shows no data, as expected) and the printouts above all work as expected, even after the call to .plot().

Maybe related to #1461.

What did you expect to happen?

I expected:

• when array.size==0, either successfully plot nothing or make a clearer error message, regardless of array.ndim.
• if crashing from plot(), array should not be altered in-place; printouts about public aspects of array (like print(array) and print(array.ndim)) before plot() crash and after crash should be the same.
• plot() should behave the same way for data that "used to be 2D+ but is now 1D" as it behaves for data that "was always 1D". I did not expect plot() to care about where the data came from, only the current state of the array.

Can you provide a MCVE to repoduce the bug?

See code above. To load the example files, I used:

import os
import uxarray as ux
def meshfile_path(*args):
    """returns abspath to uxarray/test/meshfiles, joined with *args if any provided.
    E.g., meshfiles_path("ugrid", "quad-hexagon", "grid.nc")
        --> uxarray/test/meshfiles/ugrid/quad-hexagon/grid.nc
    """
    result = os.path.join(ux.__file__, '..', '..', 'test', 'meshfiles', *args)
    return os.path.abspath(result)

ds = ux.open_dataset(meshfile_path("ugrid", "geoflow-small", "grid.nc"),
                meshfile_path("ugrid", "geoflow-small", "v1.nc"))


For the quad-hexagon example, the call looks like this instead:


data = ux.open_dataset(meshfile_path("ugrid", "quad-hexagon", "grid.nc"),
                meshfile_path("ugrid", "quad-hexagon", "data.nc"))

[Sevans711]

I noticed this seems to be related to how .subset interacts with array.load() in this case. The .load() does occur inside array.plot.points() by default, because it gets array.values.

If I call arr0.load() right after loading arr0, before doing the subsetting, The first example actually behaves just fine. Something about load() is getting confused by the subsetting in case subsetting returns array with size 0. Maybe this is primarily a subsetting issue, not a plotting issue?