Repository: pydata/xarray

Issue:
Drawing only one contour in FacetGrid
#### MCVE Code Sample
<!-- In order for the maintainers to efficiently understand and prioritize issues, we ask you post a "Minimal, Complete and Verifiable Example" (MCVE): http://matthewrocklin.com/blog/work/2018/02/28/minimal-bug-reports -->

```python
import numpy as np
import matplotlib.pyplot as plt
import xarray as xr

x, y = np.meshgrid(np.arange(12), np.arange(12))
z = xr.DataArray(np.sqrt(x**2 + y**2))
z2 = xr.DataArray(np.sqrt(x**2 + y**2)+1)
ds = xr.concat([z,z2], dim = 'time')
ds['time'] = [0,1]

f, ax = plt.subplots()
ds[0].plot.contour(ax = ax, levels=[4], colors=['k'])
f, ax = plt.subplots()
ds[1].plot.contour(ax = ax, levels=[4], colors=['k'])
ds.plot.contour(col = 'time', levels=[4], colors=['k'])
```

#### Output
![obrazek](https://user-images.githubusercontent.com/6815953/69541063-8f28b500-0f88-11ea-805c-8478f49a0fd2.png)

#### Expected Output
![obrazek](https://user-images.githubusercontent.com/6815953/69546120-6528c000-0f93-11ea-8873-71cab41eaade.png)



#### Problem Description
While drawing only one contour was fixed some time ago ([#866](https://github.com/pydata/xarray/issues/866)), single contour plots in FacetGrid do not work for me.

#### Output of ``xr.show_versions()``
<details>
INSTALLED VERSIONS
------------------
commit: None
python: 3.6.8 (default, Oct  7 2019, 12:59:55) 
[GCC 8.3.0]
python-bits: 64
OS: Linux
OS-release: 4.15.0-65-generic
machine: x86_64
processor: x86_64
byteorder: little
LC_ALL: None
LANG: en_US.UTF-8
LOCALE: en_US.UTF-8
libhdf5: 1.10.2
libnetcdf: 4.6.3

xarray: 0.14.1
pandas: 0.24.2
numpy: 1.17.3
scipy: 1.2.1
netCDF4: 1.5.1.2
pydap: installed
h5netcdf: 0.7.4
h5py: 2.9.0
Nio: None
zarr: 2.3.2
cftime: 1.0.3.4
nc_time_axis: 1.2.0
PseudoNetCDF: None
rasterio: 1.0.5
cfgrib: 0.9.6.2
iris: None
bottleneck: 1.2.1
dask: 2.6.0
distributed: 2.6.0
matplotlib: 3.0.3
cartopy: 0.16.0
seaborn: 0.9.0
numbagg: None
setuptools: 41.0.1
pip: 19.3.1
conda: None
pytest: 4.4.1
IPython: 7.1.1
sphinx: 2.0.1
</details>


Discussion:
dcherian: @mathause any thoughts on how we could test this?
mathause: The MWE in  #3569 throws an error after #3916 (`ValueError("Cannot supply vmin and a norm with a different vmin.")`) - maybe it is enough to test that it works again?

I reverted to before #3569 and tried to figure out the difference in the `QuadContourSet` in terms of norm and some other properties of `h` and `facet._mappables[0]` is and couldn't figure it out.

```python
h = ds[0].plot.contour(levels=[4], colors=['k'])
facet = ds.plot.contour(col='time', levels=[4], colors=['k'])

```


dcherian: > The MWE in #3569 throws an error after #3916 (ValueError("Cannot supply vmin and a norm with a different vmin.")) - maybe it is enough to test that it works again?

Confirmed. I have added this as a test. Thanks.

> I reverted to before #3569 and tried to figure out the difference in the QuadContourSet in terms of norm and some other properties of h and facet._mappables[0] is and couldn't figure it out.

Yes I could not figure what the difference was either...

Patch:
diff --git a/doc/examples/monthly-means.ipynb b/doc/examples/monthly-means.ipynb
index fad40e01..bc88f4a9 100644
--- a/doc/examples/monthly-means.ipynb
+++ b/doc/examples/monthly-means.ipynb
@@ -29,89 +29,9 @@
     "import numpy as np\n",
     "import pandas as pd\n",
     "import xarray as xr\n",
-    "from netCDF4 import num2date\n",
     "import matplotlib.pyplot as plt "
    ]
   },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "#### Some calendar information so we can support any netCDF calendar.  "
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "ExecuteTime": {
-     "end_time": "2018-11-28T20:51:35.991620Z",
-     "start_time": "2018-11-28T20:51:35.960336Z"
-    }
-   },
-   "outputs": [],
-   "source": [
-    "dpm = {'noleap': [0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31],\n",
-    "       '365_day': [0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31],\n",
-    "       'standard': [0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31],\n",
-    "       'gregorian': [0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31],\n",
-    "       'proleptic_gregorian': [0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31],\n",
-    "       'all_leap': [0, 31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31],\n",
-    "       '366_day': [0, 31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31],\n",
-    "       '360_day': [0, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30]} "
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "#### A few calendar functions to determine the number of days in each month\n",
-    "If you were just using the standard calendar, it would be easy to use the `calendar.month_range` function."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "ExecuteTime": {
-     "end_time": "2018-11-28T20:51:36.015151Z",
-     "start_time": "2018-11-28T20:51:35.994079Z"
-    }
-   },
-   "outputs": [],
-   "source": [
-    "def leap_year(year, calendar='standard'):\n",
-    "    \"\"\"Determine if year is a leap year\"\"\"\n",
-    "    leap = False\n",
-    "    if ((calendar in ['standard', 'gregorian',\n",
-    "        'proleptic_gregorian', 'julian']) and\n",
-    "        (year % 4 == 0)):\n",
-    "        leap = True\n",
-    "        if ((calendar == 'proleptic_gregorian') and\n",
-    "            (year % 100 == 0) and\n",
-    "            (year % 400 != 0)):\n",
-    "            leap = False\n",
-    "        elif ((calendar in ['standard', 'gregorian']) and\n",
-    "                 (year % 100 == 0) and (year % 400 != 0) and\n",
-    "                 (year < 1583)):\n",
-    "            leap = False\n",
-    "    return leap\n",
-    "\n",
-    "def get_dpm(time, calendar='standard'):\n",
-    "    \"\"\"\n",
-    "    return a array of days per month corresponding to the months provided in `months`\n",
-    "    \"\"\"\n",
-    "    month_length = np.zeros(len(time), dtype=np.int)\n",
-    "    \n",
-    "    cal_days = dpm[calendar]\n",
-    "    \n",
-    "    for i, (month, year) in enumerate(zip(time.month, time.year)):\n",
-    "        month_length[i] = cal_days[month]\n",
-    "        if leap_year(year, calendar=calendar) and month == 2:\n",
-    "            month_length[i] += 1\n",
-    "    return month_length"
-   ]
-  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -131,7 +51,7 @@
    "outputs": [],
    "source": [
     "ds = xr.tutorial.open_dataset('rasm').load()\n",
-    "print(ds)"
+    "ds"
    ]
   },
   {
@@ -143,7 +63,17 @@
     "- calculate the month lengths for each monthly data record\n",
     "- calculate weights using `groupby('time.season')`\n",
     "\n",
-    "Finally, we just need to multiply our weights by the `Dataset` and sum allong the time dimension.  "
+    "Finally, we just need to multiply our weights by the `Dataset` and sum allong the time dimension.  Creating a `DataArray` for the month length is as easy as using the `days_in_month` accessor on the time coordinate.  The calendar type, in this case `'noleap'`, is automatically considered in this operation."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "month_length = ds.time.dt.days_in_month\n",
+    "month_length"
    ]
   },
   {
@@ -157,13 +87,8 @@
    },
    "outputs": [],
    "source": [
-    "# Make a DataArray with the number of days in each month, size = len(time)\n",
-    "month_length = xr.DataArray(get_dpm(ds.time.to_index(), calendar='noleap'),\n",
-    "                            coords=[ds.time], name='month_length')\n",
-    "\n",
     "# Calculate the weights by grouping by 'time.season'.\n",
-    "# Conversion to float type ('astype(float)') only necessary for Python 2.x\n",
-    "weights = month_length.groupby('time.season') / month_length.astype(float).groupby('time.season').sum()\n",
+    "weights = month_length.groupby('time.season') / month_length.groupby('time.season').sum()\n",
     "\n",
     "# Test that the sum of the weights for each season is 1.0\n",
     "np.testing.assert_allclose(weights.groupby('time.season').sum().values, np.ones(4))\n",
@@ -183,7 +108,7 @@
    },
    "outputs": [],
    "source": [
-    "print(ds_weighted)"
+    "ds_weighted"
    ]
   },
   {
@@ -262,13 +187,9 @@
    "source": [
     "# Wrap it into a simple function\n",
     "def season_mean(ds, calendar='standard'):\n",
-    "    # Make a DataArray of season/year groups\n",
-    "    year_season = xr.DataArray(ds.time.to_index().to_period(freq='Q-NOV').to_timestamp(how='E'),\n",
-    "                               coords=[ds.time], name='year_season')\n",
-    "\n",
     "    # Make a DataArray with the number of days in each month, size = len(time)\n",
-    "    month_length = xr.DataArray(get_dpm(ds.time.to_index(), calendar=calendar),\n",
-    "                                coords=[ds.time], name='month_length')\n",
+    "    month_length = ds.time.dt.days_in_month\n",
+    "\n",
     "    # Calculate the weights by grouping by 'time.season'\n",
     "    weights = month_length.groupby('time.season') / month_length.groupby('time.season').sum()\n",
     "\n",
@@ -278,13 +199,6 @@
     "    # Calculate the weighted average\n",
     "    return (ds * weights).groupby('time.season').sum(dim='time')"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {
@@ -304,7 +218,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.8"
+   "version": "3.7.3"
   },
   "toc": {
    "base_numbering": 1,
@@ -321,5 +235,5 @@
   }
  },
  "nbformat": 4,
- "nbformat_minor": 2
+ "nbformat_minor": 4
 }
diff --git a/doc/groupby.rst b/doc/groupby.rst
index 927e192e..223185bd 100644
--- a/doc/groupby.rst
+++ b/doc/groupby.rst
@@ -157,7 +157,7 @@ This last line is roughly equivalent to the following::
 
     results = []
     for label, group in ds.groupby('letters'):
-        results.append(group - alt.sel(x=label))
+        results.append(group - alt.sel(letters=label))
     xr.concat(results, dim='x')
 
 Squeezing
diff --git a/doc/io.rst b/doc/io.rst
index 6064aa35..0c666099 100644
--- a/doc/io.rst
+++ b/doc/io.rst
@@ -281,7 +281,7 @@ library::
         # glob expands paths with * to a list of files, like the unix shell
         paths = sorted(glob(files))
         datasets = [xr.open_dataset(p) for p in paths]
-        combined = xr.concat(dataset, dim)
+        combined = xr.concat(datasets, dim)
         return combined
 
     combined = read_netcdfs('/all/my/files/*.nc', dim='time')
diff --git a/doc/roadmap.rst b/doc/roadmap.rst
index e336b35b..401dac77 100644
--- a/doc/roadmap.rst
+++ b/doc/roadmap.rst
@@ -219,6 +219,11 @@ Current core developers
 -  Fabien Maussion
 -  Keisuke Fujii
 -  Maximilian Roos
+-  Deepak Cherian
+-  Spencer Clark
+-  Tom Nicholas
+-  Guido Imperiale
+-  Justus Magin
 
 NumFOCUS
 ~~~~~~~~
diff --git a/doc/weather-climate.rst b/doc/weather-climate.rst
index 9e7c0f1d..768cf655 100644
--- a/doc/weather-climate.rst
+++ b/doc/weather-climate.rst
@@ -95,7 +95,7 @@ For data indexed by a :py:class:`~xarray.CFTimeIndex` xarray currently supports:
 
 - Access of basic datetime components via the ``dt`` accessor (in this case
   just "year", "month", "day", "hour", "minute", "second", "microsecond",
-  "season", "dayofyear", and "dayofweek"):
+  "season", "dayofyear", "dayofweek", and "days_in_month"):
 
 .. ipython:: python
 
@@ -104,6 +104,7 @@ For data indexed by a :py:class:`~xarray.CFTimeIndex` xarray currently supports:
    da.time.dt.season
    da.time.dt.dayofyear
    da.time.dt.dayofweek
+   da.time.dt.days_in_month
 
 - Rounding of datetimes to fixed frequencies via the ``dt`` accessor:
 
diff --git a/doc/whats-new.rst b/doc/whats-new.rst
index dac63600..e4c3a4d5 100644
--- a/doc/whats-new.rst
+++ b/doc/whats-new.rst
@@ -48,30 +48,52 @@ New Features
   By `Todd Jennings <https://github.com/toddrjen>`_
 - Allow plotting of boolean arrays. (:pull:`3766`)
   By `Marek Jacob <https://github.com/MeraX>`_
+- A ``days_in_month`` accessor for :py:class:`xarray.CFTimeIndex`, analogous to
+  the ``days_in_month`` accessor for a :py:class:`pandas.DatetimeIndex`, which
+  returns the days in the month each datetime in the index.  Now days in month
+  weights for both standard and non-standard calendars can be obtained using
+  the :py:class:`~core.accessor_dt.DatetimeAccessor` (:pull:`3935`).  This
+  feature requires cftime version 1.1.0 or greater.  By
+  `Spencer Clark <https://github.com/spencerkclark>`_.
 
 Bug fixes
 ~~~~~~~~~
+- Fix wrong order in converting a ``pd.Series`` with a MultiIndex to ``DataArray``. (:issue:`3951`)
+  By `Keisuke Fujii <https://github.com/fujiisoup>`_.
 - Fix renaming of coords when one or more stacked coords is not in
   sorted order during stack+groupby+apply operations. (:issue:`3287`,
   :pull:`3906`) By `Spencer Hill <https://github.com/spencerahill>`_
 - Fix a regression where deleting a coordinate from a copied :py:class:`DataArray`
   can affect the original :py:class:`Dataarray`.  (:issue:`3899`, :pull:`3871`)
   By `Todd Jennings <https://github.com/toddrjen>`_
+- Fix :py:class:`~xarray.plot.FacetGrid` plots with a single contour. (:issue:`3569`, :pull:`3915`).
+  By `Deepak Cherian <https://github.com/dcherian>`_
 - Use divergent colormap if ``levels`` spans 0. (:issue:`3524`)
   By `Deepak Cherian <https://github.com/dcherian>`_
-- Fix ``FacetGrid`` when ``vmin == vmax``. (:issue:`3734`)
+- Fix :py:class:`~xarray.plot.FacetGrid` when ``vmin == vmax``. (:issue:`3734`)
   By `Deepak Cherian <https://github.com/dcherian>`_
 - Fix bug where plotting line plots with 2D coordinates depended on dimension
   order. (:issue:`3933`)
   By `Tom Nicholas <https://github.com/TomNicholas>`_.
+- Fix ``RasterioDeprecationWarning`` when using a ``vrt`` in ``open_rasterio``. (:issue:`3964`)
+  By `Taher Chegini <https://github.com/cheginit>`_.
 
 Documentation
 ~~~~~~~~~~~~~
+- update the docstring of :py:meth:`DataArray.assign_coords` : clarify how to
+  add a new coordinate to an existing dimension and illustrative example
+  (:issue:`3952`, :pull:`3958`) By
+  `Etienne Combrisson <https://github.com/EtienneCmb>`_.
 - update the docstring of :py:meth:`Dataset.diff` and
   :py:meth:`DataArray.diff` so it does document the ``dim``
   parameter as required. (:issue:`1040`, :pull:`3909`)
   By `Justus Magin <https://github.com/keewis>`_.
-
+- Updated :doc:`Calculating Seasonal Averages from Timeseries of Monthly Means
+  <examples/monthly-means>` example notebook to take advantage of the new
+  ``days_in_month`` accessor for :py:class:`xarray.CFTimeIndex`
+  (:pull:`3935`). By `Spencer Clark <https://github.com/spencerkclark>`_.
+- Updated the list of current core developers. (:issue:`3892`)
+  By `Tom Nicholas <https://github.com/TomNicholas>`_.
 
 Internal Changes
 ~~~~~~~~~~~~~~~~
diff --git a/xarray/backends/rasterio_.py b/xarray/backends/rasterio_.py
index d041e430..77beffd0 100644
--- a/xarray/backends/rasterio_.py
+++ b/xarray/backends/rasterio_.py
@@ -224,7 +224,7 @@ def open_rasterio(filename, parse_coordinates=None, chunks=None, cache=None, loc
             crs=vrt.crs.to_string(),
             resampling=vrt.resampling,
             src_nodata=vrt.src_nodata,
-            dst_nodata=vrt.dst_nodata,
+            nodata=vrt.nodata,
             tolerance=vrt.tolerance,
             transform=vrt.transform,
             width=vrt.width,
diff --git a/xarray/coding/cftimeindex.py b/xarray/coding/cftimeindex.py
index 2e42702c..6fc28d21 100644
--- a/xarray/coding/cftimeindex.py
+++ b/xarray/coding/cftimeindex.py
@@ -243,6 +243,9 @@ class CFTimeIndex(pd.Index):
         "dayofyr", "The ordinal day of year of the datetime", "1.0.2.1"
     )
     dayofweek = _field_accessor("dayofwk", "The day of week of the datetime", "1.0.2.1")
+    days_in_month = _field_accessor(
+        "daysinmonth", "The number of days in the month of the datetime", "1.1.0.0"
+    )
     date_type = property(get_date_type)
 
     def __new__(cls, data, name=None):
diff --git a/xarray/core/common.py b/xarray/core/common.py
index a0036420..8f6d57e9 100644
--- a/xarray/core/common.py
+++ b/xarray/core/common.py
@@ -399,10 +399,14 @@ class DataWithCoords(SupportsArithmetic, AttrAccessMixin):
         Parameters
         ----------
         coords : dict, optional
-            A dict with keys which are variables names. If the values are
-            callable, they are computed on this object and assigned to new
-            coordinate variables. If the values are not callable,
-            (e.g. a ``DataArray``, scalar, or array), they are simply assigned.
+            A dict where the keys are the names of the coordinates
+            with the new values to assign. If the values are callable, they are
+            computed on this object and assigned to new coordinate variables.
+            If the values are not callable, (e.g. a ``DataArray``, scalar, or
+            array), they are simply assigned. A new coordinate can also be
+            defined and attached to an existing dimension using a tuple with
+            the first element the dimension name and the second element the
+            values for this new coordinate.
 
         **coords_kwargs : keyword, value pairs, optional
             The keyword arguments form of ``coords``.
@@ -440,10 +444,24 @@ class DataWithCoords(SupportsArithmetic, AttrAccessMixin):
         Coordinates:
           * lon      (lon) int64 -2 -1 0 1
 
+        New coordinate can also be attached to an existing dimension:
+
+        >>> lon_2 = np.array([300, 289, 0, 1])
+        >>> da.assign_coords(lon_2=('lon', lon_2))
+        <xarray.DataArray (lon: 4)>
+        array([0.28298 , 0.667347, 0.657938, 0.177683])
+        Coordinates:
+          * lon      (lon) int64 358 359 0 1
+            lon_2    (lon) int64 300 289 0 1
+
+        Note that the same result can also be obtained with a dict e.g.
+
+        >>> _ = da.assign_coords({"lon_2": ('lon', lon_2)})
+
         Notes
         -----
-        Since ``coords_kwargs`` is a dictionary, the order of your arguments may
-        not be preserved, and so the order of the new variables is not well
+        Since ``coords_kwargs`` is a dictionary, the order of your arguments
+        may not be preserved, and so the order of the new variables is not well
         defined. Assigning multiple variables within the same ``assign_coords``
         is possible, but you cannot reference other variables created within
         the same ``assign_coords`` call.
diff --git a/xarray/core/dataset.py b/xarray/core/dataset.py
index 97b3caf2..dd1e31cc 100644
--- a/xarray/core/dataset.py
+++ b/xarray/core/dataset.py
@@ -4604,8 +4604,7 @@ class Dataset(Mapping, ImplementsDatasetReduce, DataWithCoords):
         if not dataframe.columns.is_unique:
             raise ValueError("cannot convert DataFrame with non-unique columns")
 
-        idx = remove_unused_levels_categories(dataframe.index)
-        dataframe = dataframe.set_index(idx)
+        idx, dataframe = remove_unused_levels_categories(dataframe.index, dataframe)
         obj = cls()
 
         if isinstance(idx, pd.MultiIndex):
diff --git a/xarray/core/indexes.py b/xarray/core/indexes.py
index dea1767d..a4a5fa2c 100644
--- a/xarray/core/indexes.py
+++ b/xarray/core/indexes.py
@@ -9,7 +9,7 @@ from .utils import is_scalar
 from .variable import Variable
 
 
-def remove_unused_levels_categories(index):
+def remove_unused_levels_categories(index, dataframe=None):
     """
     Remove unused levels from MultiIndex and unused categories from CategoricalIndex
     """
@@ -28,7 +28,11 @@ def remove_unused_levels_categories(index):
             index = pd.MultiIndex.from_arrays(levels, names=index.names)
     elif isinstance(index, pd.CategoricalIndex):
         index = index.remove_unused_categories()
-    return index
+
+    if dataframe is None:
+        return index
+    dataframe = dataframe.set_index(index)
+    return dataframe.index, dataframe
 
 
 class Indexes(collections.abc.Mapping):
diff --git a/xarray/plot/facetgrid.py b/xarray/plot/facetgrid.py
index 4f3268c1..819eded6 100644
--- a/xarray/plot/facetgrid.py
+++ b/xarray/plot/facetgrid.py
@@ -273,7 +273,9 @@ class FacetGrid:
             # None is the sentinel value
             if d is not None:
                 subset = self.data.loc[d]
-                mappable = func(subset, x=x, y=y, ax=ax, **func_kwargs)
+                mappable = func(
+                    subset, x=x, y=y, ax=ax, **func_kwargs, _is_facetgrid=True
+                )
                 self._mappables.append(mappable)
 
         self._finalize_grid(x, y)
diff --git a/xarray/plot/plot.py b/xarray/plot/plot.py
index 0fe30283..4657bee9 100644
--- a/xarray/plot/plot.py
+++ b/xarray/plot/plot.py
@@ -693,7 +693,10 @@ def _plot2d(plotfunc):
         _ensure_plottable(xplt, yplt, zval)
 
         cmap_params, cbar_kwargs = _process_cmap_cbar_kwargs(
-            plotfunc, zval.data, **locals()
+            plotfunc,
+            zval.data,
+            **locals(),
+            _is_facetgrid=kwargs.pop("_is_facetgrid", False),
         )
 
         if "contour" in plotfunc.__name__:
diff --git a/xarray/plot/utils.py b/xarray/plot/utils.py
index a756a1c1..c3512828 100644
--- a/xarray/plot/utils.py
+++ b/xarray/plot/utils.py
@@ -153,6 +153,7 @@ def _determine_cmap_params(
     levels=None,
     filled=True,
     norm=None,
+    _is_facetgrid=False,
 ):
     """
     Use some heuristics to set good defaults for colorbar and range.
@@ -581,14 +582,15 @@ def _is_numeric(arr):
 
 
 def _add_colorbar(primitive, ax, cbar_ax, cbar_kwargs, cmap_params):
-    plt = import_matplotlib_pyplot()
+
     cbar_kwargs.setdefault("extend", cmap_params["extend"])
     if cbar_ax is None:
         cbar_kwargs.setdefault("ax", ax)
     else:
         cbar_kwargs.setdefault("cax", cbar_ax)
 
-    cbar = plt.colorbar(primitive, **cbar_kwargs)
+    fig = ax.get_figure()
+    cbar = fig.colorbar(primitive, **cbar_kwargs)
 
     return cbar
 
@@ -735,6 +737,7 @@ def _process_cmap_cbar_kwargs(
     colors=None,
     cbar_kwargs: Union[Iterable[Tuple[str, Any]], Mapping[str, Any]] = None,
     levels=None,
+    _is_facetgrid=False,
     **kwargs,
 ):
     """
@@ -781,6 +784,12 @@ def _process_cmap_cbar_kwargs(
 
     cmap_args = getfullargspec(_determine_cmap_params).args
     cmap_kwargs.update((a, kwargs[a]) for a in cmap_args if a in kwargs)
-    cmap_params = _determine_cmap_params(**cmap_kwargs)
+    if not _is_facetgrid:
+        cmap_params = _determine_cmap_params(**cmap_kwargs)
+    else:
+        cmap_params = {
+            k: cmap_kwargs[k]
+            for k in ["vmin", "vmax", "cmap", "extend", "levels", "norm"]
+        }
 
     return cmap_params, cbar_kwargs


Test patch:
diff --git a/xarray/tests/__init__.py b/xarray/tests/__init__.py
index df86b571..40c5cfa2 100644
--- a/xarray/tests/__init__.py
+++ b/xarray/tests/__init__.py
@@ -62,6 +62,7 @@ has_h5netcdf, requires_h5netcdf = _importorskip("h5netcdf")
 has_pynio, requires_pynio = _importorskip("Nio")
 has_pseudonetcdf, requires_pseudonetcdf = _importorskip("PseudoNetCDF")
 has_cftime, requires_cftime = _importorskip("cftime")
+has_cftime_1_1_0, requires_cftime_1_1_0 = _importorskip("cftime", minversion="1.1.0.0")
 has_dask, requires_dask = _importorskip("dask")
 has_bottleneck, requires_bottleneck = _importorskip("bottleneck")
 has_nc_time_axis, requires_nc_time_axis = _importorskip("nc_time_axis")
diff --git a/xarray/tests/test_cftimeindex.py b/xarray/tests/test_cftimeindex.py
index d31bf947..b30e32c9 100644
--- a/xarray/tests/test_cftimeindex.py
+++ b/xarray/tests/test_cftimeindex.py
@@ -15,7 +15,7 @@ from xarray.coding.cftimeindex import (
 )
 from xarray.tests import assert_array_equal, assert_identical
 
-from . import raises_regex, requires_cftime
+from . import raises_regex, requires_cftime, requires_cftime_1_1_0
 from .test_coding_times import (
     _ALL_CALENDARS,
     _NON_STANDARD_CALENDARS,
@@ -229,6 +229,13 @@ def test_cftimeindex_dayofweek_accessor(index):
     assert_array_equal(result, expected)
 
 
+@requires_cftime_1_1_0
+def test_cftimeindex_days_in_month_accessor(index):
+    result = index.days_in_month
+    expected = [date.daysinmonth for date in index]
+    assert_array_equal(result, expected)
+
+
 @requires_cftime
 @pytest.mark.parametrize(
     ("string", "date_args", "reso"),
diff --git a/xarray/tests/test_dataarray.py b/xarray/tests/test_dataarray.py
index cf31182e..c3e5aafa 100644
--- a/xarray/tests/test_dataarray.py
+++ b/xarray/tests/test_dataarray.py
@@ -3515,6 +3515,14 @@ class TestDataArray:
             expected_da, DataArray.from_series(actual).drop_vars(["x", "y"])
         )
 
+    def test_from_series_multiindex(self):
+        # GH:3951
+        df = pd.DataFrame({"B": [1, 2, 3], "A": [4, 5, 6]})
+        df = df.rename_axis("num").rename_axis("alpha", axis=1)
+        actual = df.stack("alpha").to_xarray()
+        assert (actual.sel(alpha="B") == [1, 2, 3]).all()
+        assert (actual.sel(alpha="A") == [4, 5, 6]).all()
+
     @requires_sparse
     def test_from_series_sparse(self):
         import sparse
@@ -4524,7 +4532,7 @@ class TestReduce1D(TestReduce):
 
     def test_idxmin(self, x, minindex, maxindex, nanindex):
         ar0 = xr.DataArray(
-            x, dims=["x"], coords={"x": np.arange(x.size) * 4}, attrs=self.attrs,
+            x, dims=["x"], coords={"x": np.arange(x.size) * 4}, attrs=self.attrs
         )
 
         # dim doesn't exist
@@ -4620,7 +4628,7 @@ class TestReduce1D(TestReduce):
 
     def test_idxmax(self, x, minindex, maxindex, nanindex):
         ar0 = xr.DataArray(
-            x, dims=["x"], coords={"x": np.arange(x.size) * 4}, attrs=self.attrs,
+            x, dims=["x"], coords={"x": np.arange(x.size) * 4}, attrs=self.attrs
         )
 
         # dim doesn't exist
diff --git a/xarray/tests/test_plot.py b/xarray/tests/test_plot.py
index d0a8a02a..bf1f9ed6 100644
--- a/xarray/tests/test_plot.py
+++ b/xarray/tests/test_plot.py
@@ -2319,3 +2319,15 @@ def test_plot_transposes_properly(plotfunc):
     # pcolormesh returns 1D array but imshow returns a 2D array so it is necessary
     # to ravel() on the LHS
     assert np.all(hdl.get_array().ravel() == da.to_masked_array().ravel())
+
+
+@requires_matplotlib
+def test_facetgrid_single_contour():
+    # regression test for GH3569
+    x, y = np.meshgrid(np.arange(12), np.arange(12))
+    z = xr.DataArray(np.sqrt(x ** 2 + y ** 2))
+    z2 = xr.DataArray(np.sqrt(x ** 2 + y ** 2) + 1)
+    ds = xr.concat([z, z2], dim="time")
+    ds["time"] = [0, 1]
+
+    ds.plot.contour(col="time", levels=[4], colors=["k"])