Repository: pydata/xarray

Issue:
"Timestamp subtraction must have the same timezones or no timezones" when saving a NetCDF
I have an issue when saving a Dataset to NetCDF. This is [the example NetCDF I am using](https://www.dropbox.com/s/ab9qcerqr3q8p12/example.nc?dl=0). 

```python
import xarray as xr
d = xr.open_dataset('example.nc')      
d.to_netcdf('out.nc')
```
Then I get:
```python
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-4-baf698f1bf45> in <module>
----> 1 d.to_netcdf('out.nc')

~/miniconda2/envs/cds/lib/python3.6/site-packages/xarray/core/dataset.py in to_netcdf(self, path, mode, format, group, engine, encoding, unlimited_dims, compute)
   1241                          engine=engine, encoding=encoding,
   1242                          unlimited_dims=unlimited_dims,
-> 1243                          compute=compute)
   1244 
   1245     def to_zarr(self, store=None, mode='w-', synchronizer=None, group=None,

~/miniconda2/envs/cds/lib/python3.6/site-packages/xarray/backends/api.py in to_netcdf(dataset, path_or_file, mode, format, group, engine, encoding, unlimited_dims, compute, multifile)
    747         # to be parallelized with dask
    748         dump_to_store(dataset, store, writer, encoding=encoding,
--> 749                       unlimited_dims=unlimited_dims)
    750         if autoclose:
    751             store.close()

~/miniconda2/envs/cds/lib/python3.6/site-packages/xarray/backends/api.py in dump_to_store(dataset, store, writer, encoder, encoding, unlimited_dims)
    790 
    791     store.store(variables, attrs, check_encoding, writer,
--> 792                 unlimited_dims=unlimited_dims)
    793 
    794 

~/miniconda2/envs/cds/lib/python3.6/site-packages/xarray/backends/common.py in store(self, variables, attributes, check_encoding_set, writer, unlimited_dims)
    259             writer = ArrayWriter()
    260 
--> 261         variables, attributes = self.encode(variables, attributes)
    262 
    263         self.set_attributes(attributes)

~/miniconda2/envs/cds/lib/python3.6/site-packages/xarray/backends/common.py in encode(self, variables, attributes)
    345         # All NetCDF files get CF encoded by default, without this attempting
    346         # to write times, for example, would fail.
--> 347         variables, attributes = cf_encoder(variables, attributes)
    348         variables = OrderedDict([(k, self.encode_variable(v))
    349                                  for k, v in variables.items()])

~/miniconda2/envs/cds/lib/python3.6/site-packages/xarray/conventions.py in cf_encoder(variables, attributes)
    603     """
    604     new_vars = OrderedDict((k, encode_cf_variable(v, name=k))
--> 605                            for k, v in iteritems(variables))
    606     return new_vars, attributes

~/miniconda2/envs/cds/lib/python3.6/site-packages/xarray/conventions.py in <genexpr>(.0)
    603     """
    604     new_vars = OrderedDict((k, encode_cf_variable(v, name=k))
--> 605                            for k, v in iteritems(variables))
    606     return new_vars, attributes

~/miniconda2/envs/cds/lib/python3.6/site-packages/xarray/conventions.py in encode_cf_variable(var, needs_copy, name)
    233                   variables.CFMaskCoder(),
    234                   variables.UnsignedIntegerCoder()]:
--> 235         var = coder.encode(var, name=name)
    236 
    237     # TODO(shoyer): convert all of these to use coders, too:

~/miniconda2/envs/cds/lib/python3.6/site-packages/xarray/coding/times.py in encode(self, variable, name)
    393                 data,
    394                 encoding.pop('units', None),
--> 395                 encoding.pop('calendar', None))
    396             safe_setitem(attrs, 'units', units, name=name)
    397             safe_setitem(attrs, 'calendar', calendar, name=name)

~/miniconda2/envs/cds/lib/python3.6/site-packages/xarray/coding/times.py in encode_cf_datetime(dates, units, calendar)
    363         # an OverflowError is raised if the ref_date is too far away from
    364         # dates to be encoded (GH 2272).
--> 365         num = (pd.DatetimeIndex(dates.ravel()) - ref_date) / time_delta
    366         num = num.values.reshape(dates.shape)
    367 

~/miniconda2/envs/cds/lib/python3.6/site-packages/pandas/core/indexes/datetimelike.py in __sub__(self, other)
    898                 result = self._add_offset(-other)
    899             elif isinstance(other, (datetime, np.datetime64)):
--> 900                 result = self._sub_datelike(other)
    901             elif is_integer(other):
    902                 # This check must come after the check for np.timedelta64

~/miniconda2/envs/cds/lib/python3.6/site-packages/pandas/core/indexes/datetimes.py in _sub_datelike(self, other)
    876             # require tz compat
    877             elif not self._has_same_tz(other):
--> 878                 raise TypeError("Timestamp subtraction must have the same "
    879                                 "timezones or no timezones")
    880             else:

TypeError: Timestamp subtraction must have the same timezones or no timezones
```

I have tried with Python 3.7 and 3.6. I have also installed the latest version of xarray hoping that this issue was linked with [#2630](https://github.com/pydata/xarray/pull/2630). Apparently, with other similar NetCDFs I don't get the error but however this is not supposed to happen, given that the same exact code was working a couple of months ago. 

<details>
INSTALLED VERSIONS
------------------
commit: None
python: 3.6.7 | packaged by conda-forge | (default, Nov 20 2018, 18:20:05) 
[GCC 4.2.1 Compatible Apple LLVM 9.0.0 (clang-900.0.37)]
python-bits: 64
OS: Darwin
OS-release: 18.2.0
machine: x86_64
processor: i386
byteorder: little
LC_ALL: en_US.UTF-8
LANG: en_US.UTF-8
LOCALE: en_US.UTF-8
libhdf5: 1.10.3
libnetcdf: 4.6.1

xarray: 0.11.1+9.g06244df
pandas: 0.23.4
numpy: 1.15.4
scipy: 1.1.0
netCDF4: 1.4.2
pydap: None
h5netcdf: None
h5py: None
Nio: None
zarr: None
cftime: 1.0.3.4
PseudonetCDF: None
rasterio: None
cfgrib: 0.9.5.1
iris: None
bottleneck: None
cyordereddict: None
dask: None
distributed: None
matplotlib: 3.0.2
cartopy: 0.17.0
seaborn: None
setuptools: 40.6.3
pip: 18.1
conda: None
pytest: None
IPython: 7.2.0
sphinx: None
</details>


Discussion:

Patch:
diff --git a/doc/whats-new.rst b/doc/whats-new.rst
index 27d0dc61..bfe6e57e 100644
--- a/doc/whats-new.rst
+++ b/doc/whats-new.rst
@@ -39,7 +39,10 @@ Bug fixes
   as an argument to ``scipy.interpolate.interp1d``, allowing for interpolation
   from higher frequencies to lower frequencies.  Datapoints outside the bounds
   of the original time coordinate are now filled with NaN (:issue:`2197`). By
-  `Spencer Clark <https://github.com/spencerkclark>`_. 
+  `Spencer Clark <https://github.com/spencerkclark>`_.
+- Saving files with times encoded with reference dates with timezones
+  (e.g. '2000-01-01T00:00:00-05:00') no longer raises an error
+  (:issue:`2649`).  By `Spencer Clark <https://github.com/spencerkclark>`_.
 
 .. _whats-new.0.11.2:
 
diff --git a/xarray/backends/api.py b/xarray/backends/api.py
index 3cd62e82..5f88783b 100644
--- a/xarray/backends/api.py
+++ b/xarray/backends/api.py
@@ -10,7 +10,8 @@ import numpy as np
 
 from .. import Dataset, backends, conventions
 from ..core import indexing
-from ..core.combine import _auto_combine, _infer_concat_order_from_positions
+from ..core.combine import (
+    _CONCAT_DIM_DEFAULT, _auto_combine, _infer_concat_order_from_positions)
 from ..core.pycompat import basestring, path_type
 from ..core.utils import close_on_error, is_grib_path, is_remote_uri
 from .common import ArrayWriter
@@ -483,9 +484,6 @@ class _MultiFileCloser(object):
             f.close()
 
 
-_CONCAT_DIM_DEFAULT = '__infer_concat_dim__'
-
-
 def open_mfdataset(paths, chunks=None, concat_dim=_CONCAT_DIM_DEFAULT,
                    compat='no_conflicts', preprocess=None, engine=None,
                    lock=None, data_vars='all', coords='different',
@@ -606,7 +604,7 @@ def open_mfdataset(paths, chunks=None, concat_dim=_CONCAT_DIM_DEFAULT,
     # Coerce 1D input into ND to maintain backwards-compatible API until API
     # for N-D combine decided
     # (see https://github.com/pydata/xarray/pull/2553/#issuecomment-445892746)
-    if concat_dim is None or concat_dim == _CONCAT_DIM_DEFAULT:
+    if concat_dim is None or concat_dim is _CONCAT_DIM_DEFAULT:
         concat_dims = concat_dim
     elif not isinstance(concat_dim, list):
         concat_dims = [concat_dim]
diff --git a/xarray/coding/times.py b/xarray/coding/times.py
index 0f2045cf..c337a42e 100644
--- a/xarray/coding/times.py
+++ b/xarray/coding/times.py
@@ -359,6 +359,11 @@ def encode_cf_datetime(dates, units=None, calendar=None):
         time_delta = np.timedelta64(1, delta_units).astype('timedelta64[ns]')
         ref_date = pd.Timestamp(ref_date)
 
+        # If the ref_date Timestamp is timezone-aware, convert to UTC and
+        # make it timezone-naive (GH 2649).
+        if ref_date.tz is not None:
+            ref_date = ref_date.tz_convert(None)
+
         # Wrap the dates in a DatetimeIndex to do the subtraction to ensure
         # an OverflowError is raised if the ref_date is too far away from
         # dates to be encoded (GH 2272).
diff --git a/xarray/core/combine.py b/xarray/core/combine.py
index abb2e2c1..e552d8d9 100644
--- a/xarray/core/combine.py
+++ b/xarray/core/combine.py
@@ -368,7 +368,7 @@ def _auto_concat(datasets, dim=None, data_vars='all', coords='different'):
         return concat(datasets, dim=dim, data_vars=data_vars, coords=coords)
 
 
-_CONCAT_DIM_DEFAULT = '__infer_concat_dim__'
+_CONCAT_DIM_DEFAULT = utils.ReprObject('<inferred>')
 
 
 def _infer_concat_order_from_positions(datasets, concat_dims):


Test patch:
diff --git a/xarray/tests/test_backends.py b/xarray/tests/test_backends.py
index 02ac0690..57e875bb 100644
--- a/xarray/tests/test_backends.py
+++ b/xarray/tests/test_backends.py
@@ -2367,6 +2367,29 @@ class TestDask(DatasetIOBase):
             with open_mfdataset([tmp], concat_dim=dim) as actual:
                 assert_identical(expected, actual)
 
+    def test_open_multi_dataset(self):
+        # Test for issue GH #1988 and #2647. This makes sure that the
+        # concat_dim is utilized when specified in open_mfdataset().
+        # The additional wrinkle is to ensure that a length greater
+        # than one is tested as well due to numpy's implicit casting
+        # of 1-length arrays to booleans in tests, which allowed
+        # #2647 to still pass the test_open_single_dataset(),
+        # which is itself still needed as-is because the original
+        # bug caused one-length arrays to not be used correctly
+        # in concatenation.
+        rnddata = np.random.randn(10)
+        original = Dataset({'foo': ('x', rnddata)})
+        dim = DataArray([100, 150], name='baz', dims='baz')
+        expected = Dataset({'foo': (('baz', 'x'),
+                                    np.tile(rnddata[np.newaxis, :], (2, 1)))},
+                           {'baz': [100, 150]})
+        with create_tmp_file() as tmp1, \
+                create_tmp_file() as tmp2:
+            original.to_netcdf(tmp1)
+            original.to_netcdf(tmp2)
+            with open_mfdataset([tmp1, tmp2], concat_dim=dim) as actual:
+                assert_identical(expected, actual)
+
     def test_dask_roundtrip(self):
         with create_tmp_file() as tmp:
             data = create_test_data()
diff --git a/xarray/tests/test_coding_times.py b/xarray/tests/test_coding_times.py
index d9a40c23..198f40ae 100644
--- a/xarray/tests/test_coding_times.py
+++ b/xarray/tests/test_coding_times.py
@@ -750,3 +750,19 @@ def test_encode_cf_datetime_pandas_min():
     np.testing.assert_array_equal(num, expected_num)
     assert units == expected_units
     assert calendar == expected_calendar
+
+
+def test_encode_cf_datetime_units_with_tz():
+    # Regression test for GH 2649
+    units = 'days since 2000-01-01T00:00:00-05:00'
+    calendar = 'proleptic_gregorian'
+    dates = pd.date_range('2000', periods=3, tz='US/Eastern').values
+    num, units, calendar = encode_cf_datetime(dates,
+                                              units=units,
+                                              calendar=calendar)
+    expected_num = np.array([0, 1, 2])
+    expected_units = 'days since 2000-01-01T00:00:00-05:00'
+    expected_calendar = 'proleptic_gregorian'
+    np.testing.assert_array_equal(num, expected_num)
+    assert units == expected_units
+    assert calendar == expected_calendar