Error when plotting NetCDF data from 1300

[SF-N]

What happened?

Hi, I am trying to plot some data from ICON climate simulations (XPP) which start in 1300.
Because of the range of np.datime46, which is 1678 to 2262, this results in the following error:

import dask
import requests
import aiohttp
import cftime
import numpy as np
import pandas as pd
import earthkit.plots
import earthkit.plots.quickmap
dask.config.set(array__chunk_size="4MiB");

BASE_URL = "https://object-store.os-api.cci1.ecmwf.int/esiwacebucket"

filename2 = "slo1802_atm_2d_ml_13000101T000000Z.nc"
import kerchunk.netCDF3
url2 = f"{BASE_URL}/ICON-XPP/{filename2}"
print(url2)
kc2 = kerchunk.netCDF3.NetCDF3ToZarr(
    url2,
    inline_threshold=0, error="raise", storage_options=dict(block_size=512),
    max_chunk_size=2**22,  # 4MiB
).translate()

ds2 = xr.open_dataset(
    "reference://", engine="zarr", backend_kwargs=dict(
        storage_options=dict(fo=kc2),
    ), consolidated=False, chunks=dict(),
)


earthkit.plots.quickmap.plot( ds2["t_s"].sel(time=cftime.DatetimeGregorian(1300, 4, 1, 0, 0, 0))[::10], 
                             units="degC", x="clon", y="clat" )

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[38], line 1
----> 1 earthkit.plots.quickmap.plot( ds2["t_s"].sel(time=cftime.DatetimeGregorian(1300, 4, 1, 0, 0, 0))[::10], 
      2                              units="degC", x="clon", y="clat" )

File [~/Documents/ECMWF/Code/earthkit-plots/src/earthkit/plots/quickmap.py:11](http://localhost:8888/lab/tree/tmp_compression_lab_notebooks/compression-lab-notebooks/04-example-datasets/~/Documents/ECMWF/Code/earthkit-plots/src/earthkit/plots/quickmap.py#line=10), in _quickmap.<locals>.wrapper(return_subplot, domain, *args, **kwargs)
      9 getattr(subplot, function.__name__)(*args, **kwargs)
     10 for method in schema.quickmap_workflow:
---> 11     getattr(subplot, method)()
     12 return subplot

File [~/Documents/ECMWF/Code/earthkit-plots/src/earthkit/plots/schemas.py:136](http://localhost:8888/lab/tree/tmp_compression_lab_notebooks/compression-lab-notebooks/04-example-datasets/~/Documents/ECMWF/Code/earthkit-plots/src/earthkit/plots/schemas.py#line=135), in Schema.apply.<locals>.decorator.<locals>.wrapper(*args, **kwargs)
    135 def wrapper(*args, **kwargs):
--> 136     return function(*args, **self._update_kwargs(kwargs, keys))

File [~/Documents/ECMWF/Code/earthkit-plots/src/earthkit/plots/components/subplots.py:917](http://localhost:8888/lab/tree/tmp_compression_lab_notebooks/compression-lab-notebooks/04-example-datasets/~/Documents/ECMWF/Code/earthkit-plots/src/earthkit/plots/components/subplots.py#line=916), in Subplot.title(self, label, unique, wrap, capitalize, **kwargs)
    915 if label is None:
    916     label = self._default_title_template
--> 917 label = self.format_string(label, unique)
    918 plt.sca(self.ax)
    919 if capitalize:

File [~/Documents/ECMWF/Code/earthkit-plots/src/earthkit/plots/components/subplots.py:945](http://localhost:8888/lab/tree/tmp_compression_lab_notebooks/compression-lab-notebooks/04-example-datasets/~/Documents/ECMWF/Code/earthkit-plots/src/earthkit/plots/components/subplots.py#line=944), in Subplot.format_string(self, string, unique, grouped)
    941     return string_utils.list_to_human(
    942         [LayerFormatter(layer).format(string) for layer in self.layers]
    943     )
    944 else:
--> 945     return SubplotFormatter(self, unique=unique).format(string)

File [~/Documents/ECMWF/Code/earthkit-plots/src/earthkit/plots/metadata/formatters.py:105](http://localhost:8888/lab/tree/tmp_compression_lab_notebooks/compression-lab-notebooks/04-example-datasets/~/Documents/ECMWF/Code/earthkit-plots/src/earthkit/plots/metadata/formatters.py#line=104), in BaseFormatter.format(self, format_string, *args, **kwargs)
    104 def format(self, format_string, [/](http://localhost:8888/), *args, **kwargs):
--> 105     kwargs = self.format_keys(format_string, kwargs)
    106     keys = list(kwargs)
    107     for key in keys:

File [~/Documents/ECMWF/Code/earthkit-plots/src/earthkit/plots/metadata/formatters.py:90](http://localhost:8888/lab/tree/tmp_compression_lab_notebooks/compression-lab-notebooks/04-example-datasets/~/Documents/ECMWF/Code/earthkit-plots/src/earthkit/plots/metadata/formatters.py#line=89), in BaseFormatter.format_keys(self, format_string, kwargs)
     88 for key in keys:
     89     main_key, *methods = key.split(".")
---> 90     result = self.format_key(main_key)
     91     for method in methods:
     92         if method in SPECIAL_METHODS:

File [~/Documents/ECMWF/Code/earthkit-plots/src/earthkit/plots/metadata/formatters.py:172](http://localhost:8888/lab/tree/tmp_compression_lab_notebooks/compression-lab-notebooks/04-example-datasets/~/Documents/ECMWF/Code/earthkit-plots/src/earthkit/plots/metadata/formatters.py#line=171), in SubplotFormatter.format_key(self, key)
    170     values = [getattr(self.subplot, self.SUBPLOT_ATTRIBUTES[key])]
    171 else:
--> 172     values = [
    173         LayerFormatter(layer).format_key(key) for layer in self.subplot.layers
    174     ]
    175 return values

File [~/Documents/ECMWF/Code/earthkit-plots/src/earthkit/plots/metadata/formatters.py:173](http://localhost:8888/lab/tree/tmp_compression_lab_notebooks/compression-lab-notebooks/04-example-datasets/~/Documents/ECMWF/Code/earthkit-plots/src/earthkit/plots/metadata/formatters.py#line=172), in <listcomp>(.0)
    170     values = [getattr(self.subplot, self.SUBPLOT_ATTRIBUTES[key])]
    171 else:
    172     values = [
--> 173         LayerFormatter(layer).format_key(key) for layer in self.subplot.layers
    174     ]
    175 return values

File [~/Documents/ECMWF/Code/earthkit-plots/src/earthkit/plots/metadata/formatters.py:145](http://localhost:8888/lab/tree/tmp_compression_lab_notebooks/compression-lab-notebooks/04-example-datasets/~/Documents/ECMWF/Code/earthkit-plots/src/earthkit/plots/metadata/formatters.py#line=144), in LayerFormatter.format_key(self, key)
    143             value = metadata.units.format_units(value)
    144 else:
--> 145     value = metadata.labels.extract(self.layer.source, key)
    146 return value

File [~/Documents/ECMWF/Code/earthkit-plots/src/earthkit/plots/metadata/labels.py:89](http://localhost:8888/lab/tree/tmp_compression_lab_notebooks/compression-lab-notebooks/04-example-datasets/~/Documents/ECMWF/Code/earthkit-plots/src/earthkit/plots/metadata/labels.py#line=88), in extract(data, attr, default)
     76 """
     77 Extract an attribute from a data object.
     78 
   (...)
     86     The default label to use if the attribute is not found.
     87 """
     88 if attr in TIME_KEYS:
---> 89     handler = TimeFormatter(data.datetime())
     90     label = getattr(handler, attr)
     91     if len(label) == 1:

File [~/Documents/ECMWF/Code/earthkit-plots/src/earthkit/plots/sources/xarray.py:86](http://localhost:8888/lab/tree/tmp_compression_lab_notebooks/compression-lab-notebooks/04-example-datasets/~/Documents/ECMWF/Code/earthkit-plots/src/earthkit/plots/sources/xarray.py#line=85), in XarraySource.datetime(self)
     84 def datetime(self):
     85     """Get the datetime of the data."""
---> 86     datetimes = [
     87         pd.to_datetime(dt).to_pydatetime()
     88         for dt in np.atleast_1d(self.data.time.values)
     89     ]
     90     return {
     91         "base_time": datetimes,
     92         "valid_time": datetimes,
     93     }

File [~/Documents/ECMWF/Code/earthkit-plots/src/earthkit/plots/sources/xarray.py:87](http://localhost:8888/lab/tree/tmp_compression_lab_notebooks/compression-lab-notebooks/04-example-datasets/~/Documents/ECMWF/Code/earthkit-plots/src/earthkit/plots/sources/xarray.py#line=86), in <listcomp>(.0)
     84 def datetime(self):
     85     """Get the datetime of the data."""
     86     datetimes = [
---> 87         pd.to_datetime(dt).to_pydatetime()
     88         for dt in np.atleast_1d(self.data.time.values)
     89     ]
     90     return {
     91         "base_time": datetimes,
     92         "valid_time": datetimes,
     93     }

File [~/Documents/ECMWF/Code/data_compression/venv/lib/python3.10/site-packages/pandas/core/tools/datetimes.py:1101](http://localhost:8888/lab/tree/tmp_compression_lab_notebooks/compression-lab-notebooks/04-example-datasets/~/Documents/ECMWF/Code/data_compression/venv/lib/python3.10/site-packages/pandas/core/tools/datetimes.py#line=1100), in to_datetime(arg, errors, dayfirst, yearfirst, utc, format, exact, unit, infer_datetime_format, origin, cache)
   1099         result = convert_listlike(argc, format)
   1100 else:
-> 1101     result = convert_listlike(np.array([arg]), format)[0]
   1102     if isinstance(arg, bool) and isinstance(result, np.bool_):
   1103         result = bool(result)  # TODO: avoid this kludge.

File [~/Documents/ECMWF/Code/data_compression/venv/lib/python3.10/site-packages/pandas/core/tools/datetimes.py:435](http://localhost:8888/lab/tree/tmp_compression_lab_notebooks/compression-lab-notebooks/04-example-datasets/~/Documents/ECMWF/Code/data_compression/venv/lib/python3.10/site-packages/pandas/core/tools/datetimes.py#line=434), in _convert_listlike_datetimes(arg, format, name, utc, unit, errors, dayfirst, yearfirst, exact)
    432 if format is not None and format != "mixed":
    433     return _array_strptime_with_fallback(arg, name, utc, format, exact, errors)
--> 435 result, tz_parsed = objects_to_datetime64(
    436     arg,
    437     dayfirst=dayfirst,
    438     yearfirst=yearfirst,
    439     utc=utc,
    440     errors=errors,
    441     allow_object=True,
    442 )
    444 if tz_parsed is not None:
    445     # We can take a shortcut since the datetime64 numpy array
    446     # is in UTC
    447     out_unit = np.datetime_data(result.dtype)[0]

File [~/Documents/ECMWF/Code/data_compression/venv/lib/python3.10/site-packages/pandas/core/arrays/datetimes.py:2398](http://localhost:8888/lab/tree/tmp_compression_lab_notebooks/compression-lab-notebooks/04-example-datasets/~/Documents/ECMWF/Code/data_compression/venv/lib/python3.10/site-packages/pandas/core/arrays/datetimes.py#line=2397), in objects_to_datetime64(data, dayfirst, yearfirst, utc, errors, allow_object, out_unit)
   2395 # if str-dtype, convert
   2396 data = np.array(data, copy=False, dtype=np.object_)
-> 2398 result, tz_parsed = tslib.array_to_datetime(
   2399     data,
   2400     errors=errors,
   2401     utc=utc,
   2402     dayfirst=dayfirst,
   2403     yearfirst=yearfirst,
   2404     creso=abbrev_to_npy_unit(out_unit),
   2405 )
   2407 if tz_parsed is not None:
   2408     # We can take a shortcut since the datetime64 numpy array
   2409     #  is in UTC
   2410     return result, tz_parsed

File tslib.pyx:414, in pandas._libs.tslib.array_to_datetime()

File tslib.pyx:596, in pandas._libs.tslib.array_to_datetime()

File tslib.pyx:588, in pandas._libs.tslib.array_to_datetime()

TypeError: <class 'cftime._cftime.DatetimeGregorian'> is not convertible to datetime, at position 0

Plotting the data from e.g. 2020 from the same simulation works fine.

The data is available at https://object-store.os-api.cci1.ecmwf.int/esiwacebucket/ICON-XPP/slo1802_atm_2d_ml_13000101T000000Z.nc.
A python notebook with the 2020 plot and the error for 1300 is available at: https://github.com/climet-eu/compression-lab-notebooks/blob/ICON-XPP_dataset/04-example-datasets/ICONXPP.ipynb

Could you please let me know how to plot the data from 1300 without any errors?

Furthermore, I am wondering why the coastline is not shown in the plots.

What are the steps to reproduce the bug?

Please run the notebook mentioned above.

Version

PR 31 (#31)

Platform (OS and architecture)

macOS 14.4.1

Relevant log output

No response

Accompanying data

No response

Organisation

No response

[SF-N]

The problem with the coastlines not showing up resloves when converting clon and clat to degrees.

[JamesVarndell]

Hi @SF-N, thanks for raising this issue. The ICON grid is currently not supported by earthkit-regrid (which is what earthkit-plots uses under-the-hood for interpolating onto a 2D map), so in your attached example, earthkit-plots is falling back to simple bilinear interpolation, treating the ICON grid points as an unstructured point cloud. This might be okay for your purposes, but I just wanted to clarify that that's what's happening in this case. If you would like more robust ICON support, you could raise an issue in the earthkit-regrid repository.

As for your datetime issue, I think that the problem is that earthkit-plots doesn't currently support cftime datetime objects. I'll raise a feature issue for this now.

Many thanks!