Plotting NEMO data with Multidimensional Coordinates breaks with ValueError

I want to plot NEMO data on tripolar ORCA grid using the xarray plot functionality. Following @rabernat's tutorial at http://xarray.pydata.org/en/stable/examples/multidimensional-coords.html, I used the x='nav_lon', y='nav_lat' to specify the (two-dimensional) coordinates that contain the longitudes and latitudes. However, I'm getting a

ValueError: ('y', 'x') must be a permuted list of ('time_counter', 'y', 'x'), unless `...` is included

See below a simple code that breaks. The ORCA025-test.nc file is 6MB in size and can be downloaded from https://surfdrive.surf.nl/files/index.php/s/DYWAOzc3WmN7DcS

> import xarray as xr
> print(xr.__version__)
0.15.1

> T = xr.open_dataset('ORCA025-test.nc').thetao.isel(deptht=0)
> T

> T.plot.pcolormesh(x='nav_lon', y='nav_lat')
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-3-ce97772e3e2b> in <module>
----> 1 T.plot.pcolormesh(x='nav_lon', y='nav_lat')

~/anaconda3/envs/py3_parcels/lib/python3.6/site-packages/xarray/plot/plot.py in plotmethod(_PlotMethods_obj, x, y, figsize, size, aspect, ax, row, col, col_wrap, xincrease, yincrease, add_colorbar, add_labels, vmin, vmax, cmap, colors, center, robust, extend, levels, infer_intervals, subplot_kws, cbar_ax, cbar_kwargs, xscale, yscale, xticks, yticks, xlim, ylim, norm, **kwargs)
    814         for arg in ["_PlotMethods_obj", "newplotfunc", "kwargs"]:
    815             del allargs[arg]
--> 816         return newplotfunc(**allargs)
    817 
    818     # Add to class _PlotMethods

~/anaconda3/envs/py3_parcels/lib/python3.6/site-packages/xarray/plot/plot.py in newplotfunc(darray, x, y, figsize, size, aspect, ax, row, col, col_wrap, xincrease, yincrease, add_colorbar, add_labels, vmin, vmax, cmap, center, robust, extend, levels, infer_intervals, colors, subplot_kws, cbar_ax, cbar_kwargs, xscale, yscale, xticks, yticks, xlim, ylim, norm, **kwargs)
    681 
    682         if dims != darray.dims:
--> 683             darray = darray.transpose(*dims, transpose_coords=True)
    684 
    685         # Pass the data as a masked ndarray too

~/anaconda3/envs/py3_parcels/lib/python3.6/site-packages/xarray/core/dataarray.py in transpose(self, transpose_coords, *dims)
   1900         """
   1901         if dims:
-> 1902             dims = tuple(utils.infix_dims(dims, self.dims))
   1903         variable = self.variable.transpose(*dims)
   1904         if transpose_coords:

~/anaconda3/envs/py3_parcels/lib/python3.6/site-packages/xarray/core/utils.py in infix_dims(dims_supplied, dims_all)
    687         if set(dims_supplied) ^ set(dims_all):
    688             raise ValueError(
--> 689                 f"{dims_supplied} must be a permuted list of {dims_all}, unless `...` is included"
    690             )
    691         yield from dims_supplied

ValueError: ('y', 'x') must be a permuted list of ('time_counter', 'y', 'x'), unless `...` is included

Also: if I add a cartopy transformation, as in @rabernat's tutorial, I get exactly the same error.

Reproduce the failure with xarray 0.15.1 using the linked ORCA025-test.nc dataset and the T.plot.pcolormesh(x='nav_lon', y='nav_lat') entry point. Start from the plotting code and the transpose path shown in the traceback; done means the multidimensional-coordinate plot no longer raises this ValueError, including when a cartopy transformation is supplied.