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"},
)
},
coords={
"rlon": (
("rlon"),
rlons,
{
"standard_name": "grid_longitude",
"units": "degrees",
"long_name": "longitude in rotated pole grid",
"axis": "X",
},
),
"rlat": (
("rlat"),
rlats,
{
"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(
pole_longitude=proj_attrs["grid_north_pole_longitude"],
pole_latitude=proj_attrs["grid_north_pole_latitude"],
central_rotated_longitude=proj_attrs["north_pole_grid_longitude"],
globe=ccrs.Globe(ellipse="WGS84"),
)
The first attempt to plot it:
ds["dummy"].hvplot.quadmesh(projection=crs_cartopy)
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
2021
(...)
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 = {}
ds2["dummy"].hvplot.quadmesh(projection=crs_cartopy)
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?