Problems plotting long model control runs with gregorian calendar

[jbusecke]

I noticed a problem the other day, when I tried to plot some very long control run data from CMIP6 on ocean.pangeo.io.

The control run of this model (CSIRO's ACCESS-ESM1-5), starts in year 101, but runs for 900 years. If I try to plot the full run as a timeseries at any point in the grid:

import xarray as xr
import matplotlib.pyplot as plt
import intake
%matplotlib inline

col = intake.open_esm_datastore("https://raw.githubusercontent.com/NCAR/intake-esm-datastore/master/catalogs/pangeo-cmip6.json")
cat = col.search(variable_id='o2', source_id='ACCESS-ESM1-5', experiment_id='piControl', table_id='Omon')
data_dict = cat.to_dataset_dict(zarr_kwargs={'consolidated': True})

# This needs a good amount of dask workers!
ds = data_dict['CMIP.CSIRO.ACCESS-ESM1-5.piControl.Omon.gn']
ds.o2.isel(i=180, j=100, lev=0).plot()

I get the following error:

Error message

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
/srv/conda/envs/notebook/lib/python3.7/site-packages/ipykernel/pylab/backend_inline.py in show(close, block)
     41             display(
     42                 figure_manager.canvas.figure,
---> 43                 metadata=_fetch_figure_metadata(figure_manager.canvas.figure)
     44             )
     45     finally:

/srv/conda/envs/notebook/lib/python3.7/site-packages/ipykernel/pylab/backend_inline.py in _fetch_figure_metadata(fig)
    179         # the background is transparent
    180         ticksLight = _is_light([label.get_color()
--> 181                                 for axes in fig.axes
    182                                 for axis in (axes.xaxis, axes.yaxis)
    183                                 for label in axis.get_ticklabels()])

/srv/conda/envs/notebook/lib/python3.7/site-packages/ipykernel/pylab/backend_inline.py in <listcomp>(.0)
    181                                 for axes in fig.axes
    182                                 for axis in (axes.xaxis, axes.yaxis)
--> 183                                 for label in axis.get_ticklabels()])
    184         if ticksLight.size and (ticksLight == ticksLight[0]).all():
    185             # there are one or more tick labels, all with the same lightness

/srv/conda/envs/notebook/lib/python3.7/site-packages/matplotlib/axis.py in get_ticklabels(self, minor, which)
   1294         if minor:
   1295             return self.get_minorticklabels()
-> 1296         return self.get_majorticklabels()
   1297 
   1298     def get_majorticklines(self):

/srv/conda/envs/notebook/lib/python3.7/site-packages/matplotlib/axis.py in get_majorticklabels(self)
   1250     def get_majorticklabels(self):
   1251         'Return a list of Text instances for the major ticklabels.'
-> 1252         ticks = self.get_major_ticks()
   1253         labels1 = [tick.label1 for tick in ticks if tick.label1.get_visible()]
   1254         labels2 = [tick.label2 for tick in ticks if tick.label2.get_visible()]

/srv/conda/envs/notebook/lib/python3.7/site-packages/matplotlib/axis.py in get_major_ticks(self, numticks)
   1405         'Get the tick instances; grow as necessary.'
   1406         if numticks is None:
-> 1407             numticks = len(self.get_majorticklocs())
   1408 
   1409         while len(self.majorTicks) < numticks:

/srv/conda/envs/notebook/lib/python3.7/site-packages/matplotlib/axis.py in get_majorticklocs(self)
   1322     def get_majorticklocs(self):
   1323         """Get the array of major tick locations in data coordinates."""
-> 1324         return self.major.locator()
   1325 
   1326     def get_minorticklocs(self):

/srv/conda/envs/notebook/lib/python3.7/site-packages/nc_time_axis/__init__.py in __call__(self)
    136     def __call__(self):
    137         vmin, vmax = self.axis.get_view_interval()
--> 138         return self.tick_values(vmin, vmax)
    139 
    140     def tick_values(self, vmin, vmax):

/srv/conda/envs/notebook/lib/python3.7/site-packages/nc_time_axis/__init__.py in tick_values(self, vmin, vmax)
    192             raise ValueError(msg)
    193 
--> 194         return utime.date2num(ticks)
    195 
    196 

cftime/_cftime.pyx in cftime._cftime.utime.date2num()

cftime/_cftime.pyx in cftime._cftime.JulianDayFromDate()

cftime/_cftime.pyx in cftime._cftime._IntJulianDayFromDate()

ValueError: year zero does not exist in the proleptic_gregorian calendar

When I restrict the timeseries to a shorter time it works fine:

ds = data_dict['CMIP.CSIRO.ACCESS-ESM1-5.piControl.Omon.gn']
ds.o2.isel(i=180, j=100, lev=0, time=slice(0,500)).plot()

I assume that the internal logic tries to set the xlimit on the left side to some negative year if the full range is large. Is there a way to suppress that behavior? Or could the plotting routine default to a leftmost date of year 1 for any cftime data with gregorian calendar?

[max-sixty]

Thanks for the report.

This isn't immediately helpful to you but: any ideas why xarray isn't in the stack trace? Is this a plotting oddity?

[mathause]

It looks like you can set xlim in da.plot.line - could you test this?

[spencerkclark]

I think this is probably best fixed in nc-time-axis. It looks like it is creating generic datetimes with year 0 somewhere in here, which cftime does not allow for the proleptic Gregorian calendar. The error gets thrown when trying to apply the calendar type in utime.date2num(ticks).

I'm fine though with you posting this here for visibility; I'm not sure how active the nc-time-axis repo is.

[jbusecke]

Apologies for not noticing this earlier. I (wrongly) assumed that xarray would handle the axis limits. Should I raise an issue over there? I am currently quite busy so it might take a while for me to be able to work on a PR.

It looks like you can set xlim in da.plot.line - could you test this?

I tried this:

# This needs a good amount of dask workers!
ds = data_dict['CMIP.CSIRO.ACCESS-ESM1-5.piControl.Omon.gn']
ts.plot(xlim=['0005', '2000'])

and am getting the same error.

[spencerkclark]

No worries! Yeah, it would be great if you could post an issue there; I think you can probably come up with a more minimal example (i.e. take xarray out of it and use some synthetic data with cftime dates near year 0 as the x-coordinate).

The fact that cftime.DatetimeProlepticGregorian supports years <= -1 and >= 1, but not year 0, makes things a bit awkward, but I think is based on CF conventions.

[jbusecke]

Closing in favor of SciTools/nc-time-axis#44