Unexpected behaviour for a CF conforming datasets

I am trying to plot a CF conforming dataset. A simplified definition:

import os
from pathlib import Path

import cartopy.crs as ccrs
import cf_xarray
import hvplot.xarray
import numpy as np
import pyproj

import xarray as xr

proj_cf = {
    "semi_major_axis": 6378137.0,
    "semi_minor_axis": 6356752.314245179,
    "inverse_flattening": 298.257223563,
    "reference_ellipsoid_name": "WGS 84",
    "longitude_of_prime_meridian": 0.0,
    "prime_meridian_name": "Greenwich",
    "geographic_crs_name": "unnamed",
    "horizontal_datum_name": "Unknown based on WGS84 ellipsoid",
    "grid_mapping_name": "rotated_latitude_longitude",
    "grid_north_pole_latitude": 23.5,
    "grid_north_pole_longitude": 15.0,
    "north_pole_grid_longitude": 0.0,

crs = pyproj.CRS.from_cf(proj_cf)
proj_cf["crs_wkt"] = crs.to_wkt()

xlim = (-2, 2)
ylim = (-2, 2)
eps = 1e-5
rlons = np.arange(xlim[0] + 1 / 2, xlim[1] + eps)
rlats = np.arange(ylim[0] + 1 / 2, ylim[1] + eps)

ds = xr.Dataset(
        "dummy": (
            ("rlats", "rlons"),
            np.ones((rlats.size, rlons.size)),
            {"grid_mapping": "rotated_latitude_longitude"},
        "rlon": (
                "standard_name": "grid_longitude",
                "units": "degrees",
                "long_name": "longitude in rotated pole grid",
                "axis": "X",
        "rlat": (
                "standard_name": "grid_latitude",
                "units": "degrees",
                "long_name": "latitude in rotated pole grid",
                "axis": "Y",
        "rotated_latitude_longitude": ((), 0, proj_cf),  # CF convention

proj_attrs = ds.coords["rotated_latitude_longitude"].attrs
crs_cartopy = ccrs.RotatedPole(

The first attempt to plot it:


Getting a traceback

DataError                                 Traceback (most recent call last)
Cell In[2], line 2
      1 'trying ds'
----> 2 ds["dummy"].hvplot.quadmesh(projection=crs_cartopy)

File ~/mambaforge/envs/icec/lib/python3.10/site-packages/hvplot/plotting/core.py:2067, in hvPlot.quadmesh(self, x, y, z, colorbar, **kwds)
   2018 def quadmesh(self, x=None, y=None, z=None, colorbar=True, **kwds):
   2019     """
   2020     QuadMesh plot
   2065     - HoloViews: https://holoviews.org/reference/elements/bokeh/QuadMesh.html
   2066     """
-> 2067     return self(x, y, z=z, kind="quadmesh", colorbar=colorbar, **kwds)

File ~/mambaforge/envs/icec/lib/python3.10/site-packages/hvplot/plotting/core.py:92, in hvPlotBase.__call__(self, x, y, kind, **kwds)
     89         plot = self._get_converter(x, y, kind, **kwds)(kind, x, y)
     90         return pn.panel(plot, **panel_dict)
---> 92 return self._get_converter(x, y, kind, **kwds)(kind, x, y)

File ~/mambaforge/envs/icec/lib/python3.10/site-packages/hvplot/converter.py:1241, in HoloViewsConverter.__call__(self, kind, x, y)
   1239     else:
   1240         name = data.name or self.label or self.value_label
-> 1241         dataset = Dataset(data, self.indexes, name)
   1242 else:
   1243     try:

File ~/mambaforge/envs/icec/lib/python3.10/site-packages/holoviews/core/data/__init__.py:331, in Dataset.__init__(self, data, kdims, vdims, **kwargs)
    328 kdims, vdims = kwargs.get('kdims'), kwargs.get('vdims')
    330 validate_vdims = kwargs.pop('_validate_vdims', True)
--> 331 initialized = Interface.initialize(type(self), data, kdims, vdims,
    332                                    datatype=kwargs.get('datatype'))
    333 (data, self.interface, dims, extra_kws) = initialized
    334 super().__init__(data, **dict(kwargs, **dict(dims, **extra_kws)))

File ~/mambaforge/envs/icec/lib/python3.10/site-packages/holoviews/core/data/interface.py:253, in Interface.initialize(cls, eltype, data, kdims, vdims, datatype)
    251     continue
    252 try:
--> 253     (data, dims, extra_kws) = interface.init(eltype, data, kdims, vdims)
    254     break
    255 except DataError:

File ~/mambaforge/envs/icec/lib/python3.10/site-packages/holoviews/core/data/xarray.py:207, in XArrayInterface.init(cls, eltype, data, kdims, vdims)
    205     raise TypeError('Data must be be an xarray Dataset type.')
    206 elif not_found:
--> 207     raise DataError("xarray Dataset must define coordinates "
    208                     "for all defined kdims, %s coordinates not found."
    209                     % not_found, cls)
    211 for vdim in vdims:
    212     if packed:

DataError: xarray Dataset must define coordinates for all defined kdims, [Dimension('rlons'), Dimension('rlats')] coordinates not found.

XArrayInterface expects gridded data, for more information on supported datatypes see http://holoviews.org/user_guide/Gridded_Datasets.html

So, it looks that holoviews can not get projection from the dataset directly.
Second attempt: adding crs=crs_cartopy to the arguments.

ds["dummy"].hvplot.quadmesh(crs=crs_cartopy, projection=crs_cartopy)

Surprisingly, the traceback is the same.

Third attempt: removing projection definition in the dataset

ds2 = ds.drop_vars('rotated_latitude_longitude')
ds2['dummy'].attrs = {}

It works, however the plot is not right. As expected, since source projection is not defined.

The final attempt produces correct output:

ds2["dummy"].hvplot.quadmesh(crs=crs_cartopy, projection=crs_cartopy)

That is, I have to specify both crs and projection arguments as cartopy projections, but also have to drop reference to the projection definition within the dataset. This looks to me as a bug, specifying crs as an argument should override search for projection in the dataset. Should I open an issue on GitHub?