docs: def infer_datetime_units(dates):
      """Given an array of datetimes, returns a CF compatible time-unit string of
      the form "{time_unit} since {date[0]}", where `time_unit` is 'days',
      'hours', 'minutes' or 'seconds' (the first one that can evenly divide all
      unique time deltas in `dates`)
      """
      dates = np.asarray(dates).ravel()
      if np.asarray(dates).dtype == 'datetime64[ns]':
          dates = pd.to_datetime(dates, box=False)
          dates = dates[pd.notnull(dates)]
          reference_date = dates[0] if len(dates) > 0 else '1970-01-01'
          reference_date = pd.Timestamp(reference_date)
      else:
          reference_date = dates[0] if len(dates) > 0 else '1970-01-01'
          reference_date = format_cftime_datetime(reference_date)
      unique_timedeltas = np.unique(np.diff(dates))
      if unique_timedeltas.dtype == np.dtype('O'):
          # Convert to np.timedelta64 objects using pandas to work around a
          # NumPy casting bug: https://github.com/numpy/numpy/issues/11096
          unique_timedeltas = pd.to_timedelta(unique_timedeltas, box=False)
      units = _infer_time_units_from_diff(unique_timedeltas)
      return '%s since %s' % (units, reference_date)