Repository: pydata/xarray

Issue:
apply_ufunc and Datasets with variables without the core dimension
### Is your feature request related to a problem?

Consider this example 

```python
ds = xr.Dataset({"a": ("x", [1, 2, 3]), "b": ("y", [1, 2, 3])})
xr.apply_ufunc(np.mean, ds, input_core_dims=[["x"]])
```

This raises 
```
ValueError: operand to apply_ufunc has required core dimensions ['x'], but some of these dimensions are absent on an input variable: ['x']
```

because core dimension `x` is missing on variable `b`. This behaviour makes it annoying to use `apply_ufunc` on Datasets.

### Describe the solution you'd like

Add a new kwarg to `apply_ufunc` called `missing_core_dim` that controls how to handle variables without *all* input core dimensions. This kwarg could take one of two values:
1. `"raise"` - raise an error, current behaviour
2. `"copy"` - skip applying the function and copy the variable from input to output.
3. `"drop"`- skip applying the function and drop the variable.

### Describe alternatives you've considered

_No response_

### Additional context

_No response_

Discussion:
max-sixty: Any thoughts on this, team?
max-sixty: What's the best name for this arg? 
- `missing_core_dim`? Possibly too synonymous with other args which actually list dims, like `exclude_dims`?
- `on_missing_core_dim`? Arguably more fluent with `on_missing_core_dim="raise"`, but less standard...
max-sixty: OK, I'm reasonably happy with this now — it's ready to go, pending docs & changelog. Any thoughts prior?
max-sixty: > We may have to think about how we handle multiple input_vars and what happens if multiple `ds` are passed (e.g `apply_ufunc(func, ds1, ds2)`). But that are probably rather edge-cases.

Sorry, just saw this. Yes, I actually had a TODO prior re adding tests for this, which I just added.

In the case of `drop` & `raise`, the order doesn't matter. In the case of `copy`, we copy from the first argument. Do you think that's reasonable?

Patch:
diff --git a/doc/whats-new.rst b/doc/whats-new.rst
index f6094103..afd493d2 100644
--- a/doc/whats-new.rst
+++ b/doc/whats-new.rst
@@ -30,6 +30,10 @@ New Features
   By `Martin Raspaud <https://github.com/mraspaud>`_.
 - Improved static typing of reduction methods (:pull:`6746`).
   By `Richard Kleijn <https://github.com/rhkleijn>`_.
+- Added `on_missing_core_dims` to :py:meth:`apply_ufunc` to allow for copying or
+  dropping a :py:class:`Dataset`'s variables with missing core dimensions.
+  (:pull:`8138`)
+  By `Maximilian Roos <https://github.com/max-sixty>`_.
 
 Breaking changes
 ~~~~~~~~~~~~~~~~
@@ -40,6 +44,9 @@ Breaking changes
   extracts and add the indexes from another :py:class:`Coordinates` object
   passed via ``coords`` (:pull:`8107`).
   By `Benoît Bovy <https://github.com/benbovy>`_.
+- Static typing of ``xlim`` and ``ylim`` arguments in plotting functions now must
+  be ``tuple[float, float]`` to align with matplotlib requirements. (:issue:`7802`, :pull:`8030`).
+  By `Michael Niklas <https://github.com/headtr1ck>`_.
 
 Deprecations
 ~~~~~~~~~~~~
@@ -70,6 +77,13 @@ Bug fixes
   special case ``NaT`` handling in :py:meth:`~core.accessor_dt.DatetimeAccessor.isocalendar()`
   (:issue:`7928`, :pull:`8084`).
   By `Kai Mühlbauer <https://github.com/kmuehlbauer>`_.
+- Calling plot with kwargs ``col``, ``row`` or ``hue`` no longer squeezes dimensions passed via these arguments
+  (:issue:`7552`, :pull:`8174`).
+  By `Wiktor Kraśnicki <https://github.com/wkrasnicki>`_.
+- Fixed a bug where casting from ``float`` to ``int64`` (undefined for ``NaN``) led to varying
+  issues (:issue:`7817`, :issue:`7942`, :issue:`7790`, :issue:`6191`, :issue:`7096`,
+  :issue:`1064`, :pull:`7827`).
+  By `Kai Mühlbauer <https://github.com/kmuehlbauer>`_.
 
 Documentation
 ~~~~~~~~~~~~~
@@ -80,6 +94,8 @@ Internal Changes
 
 - Many error messages related to invalid dimensions or coordinates now always show the list of valid dims/coords (:pull:`8079`).
   By `András Gunyhó <https://github.com/mgunyho>`_.
+- Refactor of encoding and decoding times/timedeltas to preserve nanosecond resolution in arrays that contain missing values (:pull:`7827`).
+  By `Kai Mühlbauer <https://github.com/kmuehlbauer>`_.
 
 .. _whats-new.2023.08.0:
 
diff --git a/requirements.txt b/requirements.txt
deleted file mode 100644
index 95271051..00000000
--- a/requirements.txt
+++ /dev/null
@@ -1,7 +0,0 @@
-# This file is redundant with setup.cfg;
-# it exists to let GitHub build the repository dependency graph
-# https://help.github.com/en/github/visualizing-repository-data-with-graphs/listing-the-packages-that-a-repository-depends-on
-
-numpy >= 1.21
-packaging >= 21.3
-pandas >= 1.4
diff --git a/xarray/backends/netcdf3.py b/xarray/backends/netcdf3.py
index ef389eef..db00ef19 100644
--- a/xarray/backends/netcdf3.py
+++ b/xarray/backends/netcdf3.py
@@ -88,13 +88,30 @@ def encode_nc3_attrs(attrs):
     return {k: encode_nc3_attr_value(v) for k, v in attrs.items()}
 
 
+def _maybe_prepare_times(var):
+    # checks for integer-based time-like and
+    # replaces np.iinfo(np.int64).min with _FillValue or np.nan
+    # this keeps backwards compatibility
+
+    data = var.data
+    if data.dtype.kind in "iu":
+        units = var.attrs.get("units", None)
+        if units is not None:
+            if coding.variables._is_time_like(units):
+                mask = data == np.iinfo(np.int64).min
+                if mask.any():
+                    data = np.where(mask, var.attrs.get("_FillValue", np.nan), data)
+    return data
+
+
 def encode_nc3_variable(var):
     for coder in [
         coding.strings.EncodedStringCoder(allows_unicode=False),
         coding.strings.CharacterArrayCoder(),
     ]:
         var = coder.encode(var)
-    data = coerce_nc3_dtype(var.data)
+    data = _maybe_prepare_times(var)
+    data = coerce_nc3_dtype(data)
     attrs = encode_nc3_attrs(var.attrs)
     return Variable(var.dims, data, attrs, var.encoding)
 
diff --git a/xarray/coding/times.py b/xarray/coding/times.py
index b531dc97..79efbecf 100644
--- a/xarray/coding/times.py
+++ b/xarray/coding/times.py
@@ -25,6 +25,7 @@ from xarray.core.common import contains_cftime_datetimes, is_np_datetime_like
 from xarray.core.formatting import first_n_items, format_timestamp, last_item
 from xarray.core.pdcompat import nanosecond_precision_timestamp
 from xarray.core.pycompat import is_duck_dask_array
+from xarray.core.utils import emit_user_level_warning
 from xarray.core.variable import Variable
 
 try:
@@ -122,6 +123,18 @@ def _netcdf_to_numpy_timeunit(units: str) -> str:
     }[units]
 
 
+def _numpy_to_netcdf_timeunit(units: str) -> str:
+    return {
+        "ns": "nanoseconds",
+        "us": "microseconds",
+        "ms": "milliseconds",
+        "s": "seconds",
+        "m": "minutes",
+        "h": "hours",
+        "D": "days",
+    }[units]
+
+
 def _ensure_padded_year(ref_date: str) -> str:
     # Reference dates without a padded year (e.g. since 1-1-1 or since 2-3-4)
     # are ambiguous (is it YMD or DMY?). This can lead to some very odd
@@ -171,6 +184,20 @@ def _unpack_netcdf_time_units(units: str) -> tuple[str, str]:
     return delta_units, ref_date
 
 
+def _unpack_time_units_and_ref_date(units: str) -> tuple[str, pd.Timestamp]:
+    # same us _unpack_netcdf_time_units but finalizes ref_date for
+    # processing in encode_cf_datetime
+    time_units, _ref_date = _unpack_netcdf_time_units(units)
+    # TODO: the strict enforcement of nanosecond precision Timestamps can be
+    # relaxed when addressing GitHub issue #7493.
+    ref_date = nanosecond_precision_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)
+    return time_units, ref_date
+
+
 def _decode_cf_datetime_dtype(
     data, units: str, calendar: str, use_cftime: bool | None
 ) -> np.dtype:
@@ -222,8 +249,8 @@ def _decode_datetime_with_pandas(
             "pandas."
         )
 
-    delta, ref_date = _unpack_netcdf_time_units(units)
-    delta = _netcdf_to_numpy_timeunit(delta)
+    time_units, ref_date = _unpack_netcdf_time_units(units)
+    time_units = _netcdf_to_numpy_timeunit(time_units)
     try:
         # TODO: the strict enforcement of nanosecond precision Timestamps can be
         # relaxed when addressing GitHub issue #7493.
@@ -237,8 +264,8 @@ def _decode_datetime_with_pandas(
         warnings.filterwarnings("ignore", "invalid value encountered", RuntimeWarning)
         if flat_num_dates.size > 0:
             # avoid size 0 datetimes GH1329
-            pd.to_timedelta(flat_num_dates.min(), delta) + ref_date
-            pd.to_timedelta(flat_num_dates.max(), delta) + ref_date
+            pd.to_timedelta(flat_num_dates.min(), time_units) + ref_date
+            pd.to_timedelta(flat_num_dates.max(), time_units) + ref_date
 
     # To avoid integer overflow when converting to nanosecond units for integer
     # dtypes smaller than np.int64 cast all integer and unsigned integer dtype
@@ -251,9 +278,12 @@ def _decode_datetime_with_pandas(
 
     # Cast input ordinals to integers of nanoseconds because pd.to_timedelta
     # works much faster when dealing with integers (GH 1399).
-    flat_num_dates_ns_int = (flat_num_dates * _NS_PER_TIME_DELTA[delta]).astype(
-        np.int64
-    )
+    # properly handle NaN/NaT to prevent casting NaN to int
+    nan = np.isnan(flat_num_dates) | (flat_num_dates == np.iinfo(np.int64).min)
+    flat_num_dates = flat_num_dates * _NS_PER_TIME_DELTA[time_units]
+    flat_num_dates_ns_int = np.zeros_like(flat_num_dates, dtype=np.int64)
+    flat_num_dates_ns_int[nan] = np.iinfo(np.int64).min
+    flat_num_dates_ns_int[~nan] = flat_num_dates[~nan].astype(np.int64)
 
     # Use pd.to_timedelta to safely cast integer values to timedeltas,
     # and add those to a Timestamp to safely produce a DatetimeIndex.  This
@@ -364,6 +394,10 @@ def _infer_time_units_from_diff(unique_timedeltas) -> str:
     return "seconds"
 
 
+def _time_units_to_timedelta64(units: str) -> np.timedelta64:
+    return np.timedelta64(1, _netcdf_to_numpy_timeunit(units)).astype("timedelta64[ns]")
+
+
 def infer_calendar_name(dates) -> CFCalendar:
     """Given an array of datetimes, infer the CF calendar name"""
     if is_np_datetime_like(dates.dtype):
@@ -572,9 +606,12 @@ def _should_cftime_be_used(
 
 
 def _cleanup_netcdf_time_units(units: str) -> str:
-    delta, ref_date = _unpack_netcdf_time_units(units)
+    time_units, ref_date = _unpack_netcdf_time_units(units)
+    time_units = time_units.lower()
+    if not time_units.endswith("s"):
+        time_units = f"{time_units}s"
     try:
-        units = f"{delta} since {format_timestamp(ref_date)}"
+        units = f"{time_units} since {format_timestamp(ref_date)}"
     except (OutOfBoundsDatetime, ValueError):
         # don't worry about reifying the units if they're out of bounds or
         # formatted badly
@@ -633,32 +670,32 @@ def encode_cf_datetime(
     """
     dates = np.asarray(dates)
 
+    data_units = infer_datetime_units(dates)
+
     if units is None:
-        units = infer_datetime_units(dates)
+        units = data_units
     else:
         units = _cleanup_netcdf_time_units(units)
 
     if calendar is None:
         calendar = infer_calendar_name(dates)
 
-    delta, _ref_date = _unpack_netcdf_time_units(units)
     try:
         if not _is_standard_calendar(calendar) or dates.dtype.kind == "O":
             # parse with cftime instead
             raise OutOfBoundsDatetime
         assert dates.dtype == "datetime64[ns]"
 
-        delta_units = _netcdf_to_numpy_timeunit(delta)
-        time_delta = np.timedelta64(1, delta_units).astype("timedelta64[ns]")
-
-        # TODO: the strict enforcement of nanosecond precision Timestamps can be
-        # relaxed when addressing GitHub issue #7493.
-        ref_date = nanosecond_precision_timestamp(_ref_date)
+        time_units, ref_date = _unpack_time_units_and_ref_date(units)
+        time_delta = _time_units_to_timedelta64(time_units)
 
-        # 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)
+        # retrieve needed units to faithfully encode to int64
+        needed_units, data_ref_date = _unpack_time_units_and_ref_date(data_units)
+        if data_units != units:
+            # this accounts for differences in the reference times
+            ref_delta = abs(data_ref_date - ref_date).to_timedelta64()
+            if ref_delta > np.timedelta64(0, "ns"):
+                needed_units = _infer_time_units_from_diff(ref_delta)
 
         # 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
@@ -666,29 +703,60 @@ def encode_cf_datetime(
         dates_as_index = pd.DatetimeIndex(dates.ravel())
         time_deltas = dates_as_index - ref_date
 
-        # Use floor division if time_delta evenly divides all differences
-        # to preserve integer dtype if possible (GH 4045).
-        if np.all(time_deltas % time_delta == np.timedelta64(0, "ns")):
-            num = time_deltas // time_delta
+        # needed time delta to encode faithfully to int64
+        needed_time_delta = _time_units_to_timedelta64(needed_units)
+        if time_delta <= needed_time_delta:
+            # calculate int64 floor division
+            # to preserve integer dtype if possible (GH 4045, GH7817).
+            num = time_deltas // time_delta.astype(np.int64)
+            num = num.astype(np.int64, copy=False)
         else:
+            emit_user_level_warning(
+                f"Times can't be serialized faithfully with requested units {units!r}. "
+                f"Resolution of {needed_units!r} needed. "
+                f"Serializing timeseries to floating point."
+            )
             num = time_deltas / time_delta
         num = num.values.reshape(dates.shape)
 
     except (OutOfBoundsDatetime, OverflowError, ValueError):
         num = _encode_datetime_with_cftime(dates, units, calendar)
+        # do it now only for cftime-based flow
+        # we already covered for this in pandas-based flow
+        num = cast_to_int_if_safe(num)
 
-    num = cast_to_int_if_safe(num)
     return (num, units, calendar)
 
 
 def encode_cf_timedelta(timedeltas, units: str | None = None) -> tuple[np.ndarray, str]:
-    if units is None:
-        units = infer_timedelta_units(timedeltas)
+    data_units = infer_timedelta_units(timedeltas)
 
-    np_unit = _netcdf_to_numpy_timeunit(units)
-    num = 1.0 * timedeltas / np.timedelta64(1, np_unit)
-    num = np.where(pd.isnull(timedeltas), np.nan, num)
-    num = cast_to_int_if_safe(num)
+    if units is None:
+        units = data_units
+
+    time_delta = _time_units_to_timedelta64(units)
+    time_deltas = pd.TimedeltaIndex(timedeltas.ravel())
+
+    # retrieve needed units to faithfully encode to int64
+    needed_units = data_units
+    if data_units != units:
+        needed_units = _infer_time_units_from_diff(np.unique(time_deltas.dropna()))
+
+    # needed time delta to encode faithfully to int64
+    needed_time_delta = _time_units_to_timedelta64(needed_units)
+    if time_delta <= needed_time_delta:
+        # calculate int64 floor division
+        # to preserve integer dtype if possible
+        num = time_deltas // time_delta.astype(np.int64)
+        num = num.astype(np.int64, copy=False)
+    else:
+        emit_user_level_warning(
+            f"Timedeltas can't be serialized faithfully with requested units {units!r}. "
+            f"Resolution of {needed_units!r} needed. "
+            f"Serializing timedeltas to floating point."
+        )
+        num = time_deltas / time_delta
+    num = num.values.reshape(timedeltas.shape)
     return (num, units)
 
 
@@ -702,9 +770,10 @@ class CFDatetimeCoder(VariableCoder):
         ) or contains_cftime_datetimes(variable):
             dims, data, attrs, encoding = unpack_for_encoding(variable)
 
-            (data, units, calendar) = encode_cf_datetime(
-                data, encoding.pop("units", None), encoding.pop("calendar", None)
-            )
+            units = encoding.pop("units", None)
+            calendar = encoding.pop("calendar", None)
+            (data, units, calendar) = encode_cf_datetime(data, units, calendar)
+
             safe_setitem(attrs, "units", units, name=name)
             safe_setitem(attrs, "calendar", calendar, name=name)
 
diff --git a/xarray/coding/variables.py b/xarray/coding/variables.py
index 58c4739b..d694c531 100644
--- a/xarray/coding/variables.py
+++ b/xarray/coding/variables.py
@@ -215,6 +215,21 @@ def _apply_mask(
     return np.where(condition, decoded_fill_value, data)
 
 
+def _is_time_like(units):
+    # test for time-like
+    time_strings = [
+        "since",
+        "days",
+        "hours",
+        "minutes",
+        "seconds",
+        "milliseconds",
+        "microseconds",
+        "nanoseconds",
+    ]
+    return any(tstr in str(units) for tstr in time_strings)
+
+
 class CFMaskCoder(VariableCoder):
     """Mask or unmask fill values according to CF conventions."""
 
@@ -236,19 +251,32 @@ class CFMaskCoder(VariableCoder):
                 f"Variable {name!r} has conflicting _FillValue ({fv}) and missing_value ({mv}). Cannot encode data."
             )
 
+        # special case DateTime to properly handle NaT
+        is_time_like = _is_time_like(attrs.get("units"))
+
         if fv_exists:
             # Ensure _FillValue is cast to same dtype as data's
             encoding["_FillValue"] = dtype.type(fv)
             fill_value = pop_to(encoding, attrs, "_FillValue", name=name)
             if not pd.isnull(fill_value):
-                data = duck_array_ops.fillna(data, fill_value)
+                if is_time_like and data.dtype.kind in "iu":
+                    data = duck_array_ops.where(
+                        data != np.iinfo(np.int64).min, data, fill_value
+                    )
+                else:
+                    data = duck_array_ops.fillna(data, fill_value)
 
         if mv_exists:
             # Ensure missing_value is cast to same dtype as data's
             encoding["missing_value"] = dtype.type(mv)
             fill_value = pop_to(encoding, attrs, "missing_value", name=name)
             if not pd.isnull(fill_value) and not fv_exists:
-                data = duck_array_ops.fillna(data, fill_value)
+                if is_time_like and data.dtype.kind in "iu":
+                    data = duck_array_ops.where(
+                        data != np.iinfo(np.int64).min, data, fill_value
+                    )
+                else:
+                    data = duck_array_ops.fillna(data, fill_value)
 
         return Variable(dims, data, attrs, encoding, fastpath=True)
 
@@ -275,7 +303,11 @@ class CFMaskCoder(VariableCoder):
                     stacklevel=3,
                 )
 
-            dtype, decoded_fill_value = dtypes.maybe_promote(data.dtype)
+            # special case DateTime to properly handle NaT
+            if _is_time_like(attrs.get("units")) and data.dtype.kind in "iu":
+                dtype, decoded_fill_value = np.int64, np.iinfo(np.int64).min
+            else:
+                dtype, decoded_fill_value = dtypes.maybe_promote(data.dtype)
 
             if encoded_fill_values:
                 transform = partial(
diff --git a/xarray/core/computation.py b/xarray/core/computation.py
index 2068cdca..6dfff2a2 100644
--- a/xarray/core/computation.py
+++ b/xarray/core/computation.py
@@ -32,6 +32,8 @@ if TYPE_CHECKING:
     from xarray.core.dataset import Dataset
     from xarray.core.types import CombineAttrsOptions, JoinOptions
 
+    MissingCoreDimOptions = Literal["raise", "copy", "drop"]
+
 _NO_FILL_VALUE = utils.ReprObject("<no-fill-value>")
 _DEFAULT_NAME = utils.ReprObject("<default-name>")
 _JOINS_WITHOUT_FILL_VALUES = frozenset({"inner", "exact"})
@@ -42,6 +44,7 @@ def _first_of_type(args, kind):
     for arg in args:
         if isinstance(arg, kind):
             return arg
+
     raise ValueError("This should be unreachable.")
 
 
@@ -347,7 +350,7 @@ def assert_and_return_exact_match(all_keys):
         if keys != first_keys:
             raise ValueError(
                 "exact match required for all data variable names, "
-                f"but {keys!r} != {first_keys!r}"
+                f"but {list(keys)} != {list(first_keys)}: {set(keys) ^ set(first_keys)} are not in both."
             )
     return first_keys
 
@@ -376,7 +379,7 @@ def collect_dict_values(
     ]
 
 
-def _as_variables_or_variable(arg):
+def _as_variables_or_variable(arg) -> Variable | tuple[Variable]:
     try:
         return arg.variables
     except AttributeError:
@@ -396,8 +399,39 @@ def _unpack_dict_tuples(
     return out
 
 
+def _check_core_dims(signature, variable_args, name):
+    """
+    Chcek if an arg has all the core dims required by the signature.
+
+    Slightly awkward design, of returning the error message. But we want to
+    give a detailed error message, which requires inspecting the variable in
+    the inner loop.
+    """
+    missing = []
+    for i, (core_dims, variable_arg) in enumerate(
+        zip(signature.input_core_dims, variable_args)
+    ):
+        # Check whether all the dims are on the variable. Note that we need the
+        # `hasattr` to check for a dims property, to protect against the case where
+        # a numpy array is passed in.
+        if hasattr(variable_arg, "dims") and set(core_dims) - set(variable_arg.dims):
+            missing += [[i, variable_arg, core_dims]]
+    if missing:
+        message = ""
+        for i, variable_arg, core_dims in missing:
+            message += f"Missing core dims {set(core_dims) - set(variable_arg.dims)} from arg number {i + 1} on a variable named `{name}`:\n{variable_arg}\n\n"
+        message += "Either add the core dimension, or if passing a dataset alternatively pass `on_missing_core_dim` as `copy` or `drop`. "
+        return message
+    return True
+
+
 def apply_dict_of_variables_vfunc(
-    func, *args, signature: _UFuncSignature, join="inner", fill_value=None
+    func,
+    *args,
+    signature: _UFuncSignature,
+    join="inner",
+    fill_value=None,
+    on_missing_core_dim: MissingCoreDimOptions = "raise",
 ):
     """Apply a variable level function over dicts of DataArray, DataArray,
     Variable and ndarray objects.
@@ -408,7 +442,20 @@ def apply_dict_of_variables_vfunc(
 
     result_vars = {}
     for name, variable_args in zip(names, grouped_by_name):
-        result_vars[name] = func(*variable_args)
+        core_dim_present = _check_core_dims(signature, variable_args, name)
+        if core_dim_present is True:
+            result_vars[name] = func(*variable_args)
+        else:
+            if on_missing_core_dim == "raise":
+                raise ValueError(core_dim_present)
+            elif on_missing_core_dim == "copy":
+                result_vars[name] = variable_args[0]
+            elif on_missing_core_dim == "drop":
+                pass
+            else:
+                raise ValueError(
+                    f"Invalid value for `on_missing_core_dim`: {on_missing_core_dim!r}"
+                )
 
     if signature.num_outputs > 1:
         return _unpack_dict_tuples(result_vars, signature.num_outputs)
@@ -441,6 +488,7 @@ def apply_dataset_vfunc(
     fill_value=_NO_FILL_VALUE,
     exclude_dims=frozenset(),
     keep_attrs="override",
+    on_missing_core_dim: MissingCoreDimOptions = "raise",
 ) -> Dataset | tuple[Dataset, ...]:
     """Apply a variable level function over Dataset, dict of DataArray,
     DataArray, Variable and/or ndarray objects.
@@ -467,7 +515,12 @@ def apply_dataset_vfunc(
     args = tuple(getattr(arg, "data_vars", arg) for arg in args)
 
     result_vars = apply_dict_of_variables_vfunc(
-        func, *args, signature=signature, join=dataset_join, fill_value=fill_value
+        func,
+        *args,
+        signature=signature,
+        join=dataset_join,
+        fill_value=fill_value,
+        on_missing_core_dim=on_missing_core_dim,
     )
 
     out: Dataset | tuple[Dataset, ...]
@@ -595,17 +648,9 @@ def broadcast_compat_data(
         return data
 
     set_old_dims = set(old_dims)
-    missing_core_dims = [d for d in core_dims if d not in set_old_dims]
-    if missing_core_dims:
-        raise ValueError(
-            "operand to apply_ufunc has required core dimensions {}, but "
-            "some of these dimensions are absent on an input variable: {}".format(
-                list(core_dims), missing_core_dims
-            )
-        )
-
     set_new_dims = set(new_dims)
     unexpected_dims = [d for d in old_dims if d not in set_new_dims]
+
     if unexpected_dims:
         raise ValueError(
             "operand to apply_ufunc encountered unexpected "
@@ -851,6 +896,7 @@ def apply_ufunc(
     output_sizes: Mapping[Any, int] | None = None,
     meta: Any = None,
     dask_gufunc_kwargs: dict[str, Any] | None = None,
+    on_missing_core_dim: MissingCoreDimOptions = "raise",
 ) -> Any:
     """Apply a vectorized function for unlabeled arrays on xarray objects.
 
@@ -964,6 +1010,8 @@ def apply_ufunc(
         :py:func:`dask.array.apply_gufunc`. ``meta`` should be given in the
         ``dask_gufunc_kwargs`` parameter . It will be removed as direct parameter
         a future version.
+    on_missing_core_dim : {"raise", "copy", "drop"}, default: "raise"
+        How to handle missing core dimensions on input variables.
 
     Returns
     -------
@@ -1192,6 +1240,7 @@ def apply_ufunc(
             dataset_join=dataset_join,
             fill_value=dataset_fill_value,
             keep_attrs=keep_attrs,
+            on_missing_core_dim=on_missing_core_dim,
         )
     # feed DataArray apply_variable_ufunc through apply_dataarray_vfunc
     elif any(isinstance(a, DataArray) for a in args):
diff --git a/xarray/core/options.py b/xarray/core/options.py
index eb0c56c7..a197cb4d 100644
--- a/xarray/core/options.py
+++ b/xarray/core/options.py
@@ -6,10 +6,8 @@ from typing import TYPE_CHECKING, Literal, TypedDict
 from xarray.core.utils import FrozenDict
 
 if TYPE_CHECKING:
-    try:
-        from matplotlib.colors import Colormap
-    except ImportError:
-        Colormap = str
+    from matplotlib.colors import Colormap
+
     Options = Literal[
         "arithmetic_join",
         "cmap_divergent",
@@ -164,11 +162,11 @@ class set_options:
     cmap_divergent : str or matplotlib.colors.Colormap, default: "RdBu_r"
         Colormap to use for divergent data plots. If string, must be
         matplotlib built-in colormap. Can also be a Colormap object
-        (e.g. mpl.cm.magma)
+        (e.g. mpl.colormaps["magma"])
     cmap_sequential : str or matplotlib.colors.Colormap, default: "viridis"
         Colormap to use for nondivergent data plots. If string, must be
         matplotlib built-in colormap. Can also be a Colormap object
-        (e.g. mpl.cm.magma)
+        (e.g. mpl.colormaps["magma"])
     display_expand_attrs : {"default", True, False}
         Whether to expand the attributes section for display of
         ``DataArray`` or ``Dataset`` objects. Can be
diff --git a/xarray/plot/accessor.py b/xarray/plot/accessor.py
index ff707602..203bae26 100644
--- a/xarray/plot/accessor.py
+++ b/xarray/plot/accessor.py
@@ -16,6 +16,7 @@ if TYPE_CHECKING:
     from matplotlib.container import BarContainer
     from matplotlib.contour import QuadContourSet
     from matplotlib.image import AxesImage
+    from matplotlib.patches import Polygon
     from matplotlib.quiver import Quiver
     from mpl_toolkits.mplot3d.art3d import Line3D, Poly3DCollection
     from numpy.typing import ArrayLike
@@ -47,11 +48,13 @@ class DataArrayPlotAccessor:
         return dataarray_plot.plot(self._da, **kwargs)
 
     @functools.wraps(dataarray_plot.hist)
-    def hist(self, *args, **kwargs) -> tuple[np.ndarray, np.ndarray, BarContainer]:
+    def hist(
+        self, *args, **kwargs
+    ) -> tuple[np.ndarray, np.ndarray, BarContainer | Polygon]:
         return dataarray_plot.hist(self._da, *args, **kwargs)
 
     @overload
-    def line(  # type: ignore[misc]  # None is hashable :(
+    def line(  # type: ignore[misc,unused-ignore]  # None is hashable :(
         self,
         *args: Any,
         row: None = None,  # no wrap -> primitive
@@ -69,8 +72,8 @@ class DataArrayPlotAccessor:
         yscale: ScaleOptions = None,
         xticks: ArrayLike | None = None,
         yticks: ArrayLike | None = None,
-        xlim: ArrayLike | None = None,
-        ylim: ArrayLike | None = None,
+        xlim: tuple[float, float] | None = None,
+        ylim: tuple[float, float] | None = None,
         add_legend: bool = True,
         _labels: bool = True,
         **kwargs: Any,
@@ -96,8 +99,8 @@ class DataArrayPlotAccessor:
         yscale: ScaleOptions = None,
         xticks: ArrayLike | None = None,
         yticks: ArrayLike | None = None,
-        xlim: ArrayLike | None = None,
-        ylim: ArrayLike | None = None,
+        xlim: tuple[float, float] | None = None,
+        ylim: tuple[float, float] | None = None,
         add_legend: bool = True,
         _labels: bool = True,
         **kwargs: Any,
@@ -123,20 +126,20 @@ class DataArrayPlotAccessor:
         yscale: ScaleOptions = None,
         xticks: ArrayLike | None = None,
         yticks: ArrayLike | None = None,
-        xlim: ArrayLike | None = None,
-        ylim: ArrayLike | None = None,
+        xlim: tuple[float, float] | None = None,
+        ylim: tuple[float, float] | None = None,
         add_legend: bool = True,
         _labels: bool = True,
         **kwargs: Any,
     ) -> FacetGrid[DataArray]:
         ...
 
-    @functools.wraps(dataarray_plot.line)
+    @functools.wraps(dataarray_plot.line, assigned=("__doc__",))
     def line(self, *args, **kwargs) -> list[Line3D] | FacetGrid[DataArray]:
         return dataarray_plot.line(self._da, *args, **kwargs)
 
     @overload
-    def step(  # type: ignore[misc]  # None is hashable :(
+    def step(  # type: ignore[misc,unused-ignore]  # None is hashable :(
         self,
         *args: Any,
         where: Literal["pre", "post", "mid"] = "pre",
@@ -174,12 +177,12 @@ class DataArrayPlotAccessor:
     ) -> FacetGrid[DataArray]:
         ...
 
-    @functools.wraps(dataarray_plot.step)
+    @functools.wraps(dataarray_plot.step, assigned=("__doc__",))
     def step(self, *args, **kwargs) -> list[Line3D] | FacetGrid[DataArray]:
         return dataarray_plot.step(self._da, *args, **kwargs)
 
     @overload
-    def scatter(
+    def scatter(  # type: ignore[misc,unused-ignore]  # None is hashable :(
         self,
         *args: Any,
         x: Hashable | None = None,
@@ -207,8 +210,8 @@ class DataArrayPlotAccessor:
         yscale: ScaleOptions = None,
         xticks: ArrayLike | None = None,
         yticks: ArrayLike | None = None,
-        xlim: ArrayLike | None = None,
-        ylim: ArrayLike | None = None,
+        xlim: tuple[float, float] | None = None,
+        ylim: tuple[float, float] | None = None,
         cmap=None,
         vmin: float | None = None,
         vmax: float | None = None,
@@ -248,8 +251,8 @@ class DataArrayPlotAccessor:
         yscale: ScaleOptions = None,
         xticks: ArrayLike | None = None,
         yticks: ArrayLike | None = None,
-        xlim: ArrayLike | None = None,
-        ylim: ArrayLike | None = None,
+        xlim: tuple[float, float] | None = None,
+        ylim: tuple[float, float] | None = None,
         cmap=None,
         vmin: float | None = None,
         vmax: float | None = None,
@@ -289,8 +292,8 @@ class DataArrayPlotAccessor:
         yscale: ScaleOptions = None,
         xticks: ArrayLike | None = None,
         yticks: ArrayLike | None = None,
-        xlim: ArrayLike | None = None,
-        ylim: ArrayLike | None = None,
+        xlim: tuple[float, float] | None = None,
+        ylim: tuple[float, float] | None = None,
         cmap=None,
         vmin: float | None = None,
         vmax: float | None = None,
@@ -301,12 +304,12 @@ class DataArrayPlotAccessor:
     ) -> FacetGrid[DataArray]:
         ...
 
-    @functools.wraps(dataarray_plot.scatter)
-    def scatter(self, *args, **kwargs):
+    @functools.wraps(dataarray_plot.scatter, assigned=("__doc__",))
+    def scatter(self, *args, **kwargs) -> PathCollection | FacetGrid[DataArray]:
         return dataarray_plot.scatter(self._da, *args, **kwargs)
 
     @overload
-    def imshow(
+    def imshow(  # type: ignore[misc,unused-ignore]  # None is hashable :(
         self,
         *args: Any,
         x: Hashable | None = None,
@@ -338,8 +341,8 @@ class DataArrayPlotAccessor:
         yscale: ScaleOptions = None,
         xticks: ArrayLike | None = None,
         yticks: ArrayLike | None = None,
-        xlim: ArrayLike | None = None,
-        ylim: ArrayLike | None = None,
+        xlim: tuple[float, float] | None = None,
+        ylim: tuple[float, float] | None = None,
         norm: Normalize | None = None,
         **kwargs: Any,
     ) -> AxesImage:
@@ -378,8 +381,8 @@ class DataArrayPlotAccessor:
         yscale: ScaleOptions = None,
         xticks: ArrayLike | None = None,
         yticks: ArrayLike | None = None,
-        xlim: ArrayLike | None = None,
-        ylim: ArrayLike | None = None,
+        xlim: tuple[float, float] | None = None,
+        ylim: tuple[float, float] | None = None,
         norm: Normalize | None = None,
         **kwargs: Any,
     ) -> FacetGrid[DataArray]:
@@ -418,19 +421,19 @@ class DataArrayPlotAccessor:
         yscale: ScaleOptions = None,
         xticks: ArrayLike | None = None,
         yticks: ArrayLike | None = None,
-        xlim: ArrayLike | None = None,
-        ylim: ArrayLike | None = None,
+        xlim: tuple[float, float] | None = None,
+        ylim: tuple[float, float] | None = None,
         norm: Normalize | None = None,
         **kwargs: Any,
     ) -> FacetGrid[DataArray]:
         ...
 
-    @functools.wraps(dataarray_plot.imshow)
-    def imshow(self, *args, **kwargs) -> AxesImage:
+    @functools.wraps(dataarray_plot.imshow, assigned=("__doc__",))
+    def imshow(self, *args, **kwargs) -> AxesImage | FacetGrid[DataArray]:
         return dataarray_plot.imshow(self._da, *args, **kwargs)
 
     @overload
-    def contour(
+    def contour(  # type: ignore[misc,unused-ignore]  # None is hashable :(
         self,
         *args: Any,
         x: Hashable | None = None,
@@ -462,8 +465,8 @@ class DataArrayPlotAccessor:
         yscale: ScaleOptions = None,
         xticks: ArrayLike | None = None,
         yticks: ArrayLike | None = None,
-        xlim: ArrayLike | None = None,
-        ylim: ArrayLike | None = None,
+        xlim: tuple[float, float] | None = None,
+        ylim: tuple[float, float] | None = None,
         norm: Normalize | None = None,
         **kwargs: Any,
     ) -> QuadContourSet:
@@ -502,8 +505,8 @@ class DataArrayPlotAccessor:
         yscale: ScaleOptions = None,
         xticks: ArrayLike | None = None,
         yticks: ArrayLike | None = None,
-        xlim: ArrayLike | None = None,
-        ylim: ArrayLike | None = None,
+        xlim: tuple[float, float] | None = None,
+        ylim: tuple[float, float] | None = None,
         norm: Normalize | None = None,
         **kwargs: Any,
     ) -> FacetGrid[DataArray]:
@@ -542,19 +545,19 @@ class DataArrayPlotAccessor:
         yscale: ScaleOptions = None,
         xticks: ArrayLike | None = None,
         yticks: ArrayLike | None = None,
-        xlim: ArrayLike | None = None,
-        ylim: ArrayLike | None = None,
+        xlim: tuple[float, float] | None = None,
+        ylim: tuple[float, float] | None = None,
         norm: Normalize | None = None,
         **kwargs: Any,
     ) -> FacetGrid[DataArray]:
         ...
 
-    @functools.wraps(dataarray_plot.contour)
-    def contour(self, *args, **kwargs) -> QuadContourSet:
+    @functools.wraps(dataarray_plot.contour, assigned=("__doc__",))
+    def contour(self, *args, **kwargs) -> QuadContourSet | FacetGrid[DataArray]:
         return dataarray_plot.contour(self._da, *args, **kwargs)
 
     @overload
-    def contourf(
+    def contourf(  # type: ignore[misc,unused-ignore]  # None is hashable :(
         self,
         *args: Any,
         x: Hashable | None = None,
@@ -586,8 +589,8 @@ class DataArrayPlotAccessor:
         yscale: ScaleOptions = None,
         xticks: ArrayLike | None = None,
         yticks: ArrayLike | None = None,
-        xlim: ArrayLike | None = None,
-        ylim: ArrayLike | None = None,
+        xlim: tuple[float, float] | None = None,
+        ylim: tuple[float, float] | None = None,
         norm: Normalize | None = None,
         **kwargs: Any,
     ) -> QuadContourSet:
@@ -626,8 +629,8 @@ class DataArrayPlotAccessor:
         yscale: ScaleOptions = None,
         xticks: ArrayLike | None = None,
         yticks: ArrayLike | None = None,
-        xlim: ArrayLike | None = None,
-        ylim: ArrayLike | None = None,
+        xlim: tuple[float, float] | None = None,
+        ylim: tuple[float, float] | None = None,
         norm: Normalize | None = None,
         **kwargs: Any,
     ) -> FacetGrid[DataArray]:
@@ -666,19 +669,19 @@ class DataArrayPlotAccessor:
         yscale: ScaleOptions = None,
         xticks: ArrayLike | None = None,
         yticks: ArrayLike | None = None,
-        xlim: ArrayLike | None = None,
-        ylim: ArrayLike | None = None,
+        xlim: tuple[float, float] | None = None,
+        ylim: tuple[float, float] | None = None,
         norm: Normalize | None = None,
         **kwargs: Any,
     ) -> FacetGrid:
         ...
 
-    @functools.wraps(dataarray_plot.contourf)
-    def contourf(self, *args, **kwargs) -> QuadContourSet:
+    @functools.wraps(dataarray_plot.contourf, assigned=("__doc__",))
+    def contourf(self, *args, **kwargs) -> QuadContourSet | FacetGrid[DataArray]:
         return dataarray_plot.contourf(self._da, *args, **kwargs)
 
     @overload
-    def pcolormesh(
+    def pcolormesh(  # type: ignore[misc,unused-ignore]  # None is hashable :(
         self,
         *args: Any,
         x: Hashable | None = None,
@@ -710,8 +713,8 @@ class DataArrayPlotAccessor:
         yscale: ScaleOptions = None,
         xticks: ArrayLike | None = None,
         yticks: ArrayLike | None = None,
-        xlim: ArrayLike | None = None,
-        ylim: ArrayLike | None = None,
+        xlim: tuple[float, float] | None = None,
+        ylim: tuple[float, float] | None = None,
         norm: Normalize | None = None,
         **kwargs: Any,
     ) -> QuadMesh:
@@ -750,11 +753,11 @@ class DataArrayPlotAccessor:
         yscale: ScaleOptions = None,
         xticks: ArrayLike | None = None,
         yticks: ArrayLike | None = None,
-        xlim: ArrayLike | None = None,
-        ylim: ArrayLike | None = None,
+        xlim: tuple[float, float] | None = None,
+        ylim: tuple[float, float] | None = None,
         norm: Normalize | None = None,
         **kwargs: Any,
-    ) -> FacetGrid:
+    ) -> FacetGrid[DataArray]:
         ...
 
     @overload
@@ -790,15 +793,15 @@ class DataArrayPlotAccessor:
         yscale: ScaleOptions = None,
         xticks: ArrayLike | None = None,
         yticks: ArrayLike | None = None,
-        xlim: ArrayLike | None = None,
-        ylim: ArrayLike | None = None,
+        xlim: tuple[float, float] | None = None,
+        ylim: tuple[float, float] | None = None,
         norm: Normalize | None = None,
         **kwargs: Any,
-    ) -> FacetGrid:
+    ) -> FacetGrid[DataArray]:
         ...
 
-    @functools.wraps(dataarray_plot.pcolormesh)
-    def pcolormesh(self, *args, **kwargs) -> QuadMesh:
+    @functools.wraps(dataarray_plot.pcolormesh, assigned=("__doc__",))
+    def pcolormesh(self, *args, **kwargs) -> QuadMesh | FacetGrid[DataArray]:
         return dataarray_plot.pcolormesh(self._da, *args, **kwargs)
 
     @overload
@@ -834,8 +837,8 @@ class DataArrayPlotAccessor:
         yscale: ScaleOptions = None,
         xticks: ArrayLike | None = None,
         yticks: ArrayLike | None = None,
-        xlim: ArrayLike | None = None,
-        ylim: ArrayLike | None = None,
+        xlim: tuple[float, float] | None = None,
+        ylim: tuple[float, float] | None = None,
         norm: Normalize | None = None,
         **kwargs: Any,
     ) -> Poly3DCollection:
@@ -874,8 +877,8 @@ class DataArrayPlotAccessor:
         yscale: ScaleOptions = None,
         xticks: ArrayLike | None = None,
         yticks: ArrayLike | None = None,
-        xlim: ArrayLike | None = None,
-        ylim: ArrayLike | None = None,
+        xlim: tuple[float, float] | None = None,
+        ylim: tuple[float, float] | None = None,
         norm: Normalize | None = None,
         **kwargs: Any,
     ) -> FacetGrid:
@@ -914,14 +917,14 @@ class DataArrayPlotAccessor:
         yscale: ScaleOptions = None,
         xticks: ArrayLike | None = None,
         yticks: ArrayLike | None = None,
-        xlim: ArrayLike | None = None,
-        ylim: ArrayLike | None = None,
+        xlim: tuple[float, float] | None = None,
+        ylim: tuple[float, float] | None = None,
         norm: Normalize | None = None,
         **kwargs: Any,
     ) -> FacetGrid:
         ...
 
-    @functools.wraps(dataarray_plot.surface)
+    @functools.wraps(dataarray_plot.surface, assigned=("__doc__",))
     def surface(self, *args, **kwargs) -> Poly3DCollection:
         return dataarray_plot.surface(self._da, *args, **kwargs)
 
@@ -945,7 +948,7 @@ class DatasetPlotAccessor:
         )
 
     @overload
-    def scatter(
+    def scatter(  # type: ignore[misc,unused-ignore]  # None is hashable :(
         self,
         *args: Any,
         x: Hashable | None = None,
@@ -973,8 +976,8 @@ class DatasetPlotAccessor:
         yscale: ScaleOptions = None,
         xticks: ArrayLike | None = None,
         yticks: ArrayLike | None = None,
-        xlim: ArrayLike | None = None,
-        ylim: ArrayLike | None = None,
+        xlim: tuple[float, float] | None = None,
+        ylim: tuple[float, float] | None = None,
         cmap=None,
         vmin: float | None = None,
         vmax: float | None = None,
@@ -1014,8 +1017,8 @@ class DatasetPlotAccessor:
         yscale: ScaleOptions = None,
         xticks: ArrayLike | None = None,
         yticks: ArrayLike | None = None,
-        xlim: ArrayLike | None = None,
-        ylim: ArrayLike | None = None,
+        xlim: tuple[float, float] | None = None,
+        ylim: tuple[float, float] | None = None,
         cmap=None,
         vmin: float | None = None,
         vmax: float | None = None,
@@ -1023,7 +1026,7 @@ class DatasetPlotAccessor:
         extend=None,
         levels=None,
         **kwargs: Any,
-    ) -> FacetGrid[DataArray]:
+    ) -> FacetGrid[Dataset]:
         ...
 
     @overload
@@ -1055,8 +1058,8 @@ class DatasetPlotAccessor:
         yscale: ScaleOptions = None,
         xticks: ArrayLike | None = None,
         yticks: ArrayLike | None = None,
-        xlim: ArrayLike | None = None,
-        ylim: ArrayLike | None = None,
+        xlim: tuple[float, float] | None = None,
+        ylim: tuple[float, float] | None = None,
         cmap=None,
         vmin: float | None = None,
         vmax: float | None = None,
@@ -1064,15 +1067,15 @@ class DatasetPlotAccessor:
         extend=None,
         levels=None,
         **kwargs: Any,
-    ) -> FacetGrid[DataArray]:
+    ) -> FacetGrid[Dataset]:
         ...
 
-    @functools.wraps(dataset_plot.scatter)
-    def scatter(self, *args, **kwargs) -> PathCollection | FacetGrid[DataArray]:
+    @functools.wraps(dataset_plot.scatter, assigned=("__doc__",))
+    def scatter(self, *args, **kwargs) -> PathCollection | FacetGrid[Dataset]:
         return dataset_plot.scatter(self._ds, *args, **kwargs)
 
     @overload
-    def quiver(
+    def quiver(  # type: ignore[misc,unused-ignore]  # None is hashable :(
         self,
         *args: Any,
         x: Hashable | None = None,
@@ -1142,7 +1145,7 @@ class DatasetPlotAccessor:
         extend=None,
         cmap=None,
         **kwargs: Any,
-    ) -> FacetGrid:
+    ) -> FacetGrid[Dataset]:
         ...
 
     @overload
@@ -1179,15 +1182,15 @@ class DatasetPlotAccessor:
         extend=None,
         cmap=None,
         **kwargs: Any,
-    ) -> FacetGrid:
+    ) -> FacetGrid[Dataset]:
         ...
 
-    @functools.wraps(dataset_plot.quiver)
-    def quiver(self, *args, **kwargs) -> Quiver | FacetGrid:
+    @functools.wraps(dataset_plot.quiver, assigned=("__doc__",))
+    def quiver(self, *args, **kwargs) -> Quiver | FacetGrid[Dataset]:
         return dataset_plot.quiver(self._ds, *args, **kwargs)
 
     @overload
-    def streamplot(
+    def streamplot(  # type: ignore[misc,unused-ignore]  # None is hashable :(
         self,
         *args: Any,
         x: Hashable | None = None,
@@ -1257,7 +1260,7 @@ class DatasetPlotAccessor:
         extend=None,
         cmap=None,
         **kwargs: Any,
-    ) -> FacetGrid:
+    ) -> FacetGrid[Dataset]:
         ...
 
     @overload
@@ -1294,9 +1297,9 @@ class DatasetPlotAccessor:
         extend=None,
         cmap=None,
         **kwargs: Any,
-    ) -> FacetGrid:
+    ) -> FacetGrid[Dataset]:
         ...
 
-    @functools.wraps(dataset_plot.streamplot)
-    def streamplot(self, *args, **kwargs) -> LineCollection | FacetGrid:
+    @functools.wraps(dataset_plot.streamplot, assigned=("__doc__",))
+    def streamplot(self, *args, **kwargs) -> LineCollection | FacetGrid[Dataset]:
         return dataset_plot.streamplot(self._ds, *args, **kwargs)
diff --git a/xarray/plot/dataarray_plot.py b/xarray/plot/dataarray_plot.py
index 1e932479..54f3c657 100644
--- a/xarray/plot/dataarray_plot.py
+++ b/xarray/plot/dataarray_plot.py
@@ -3,7 +3,7 @@ from __future__ import annotations
 import functools
 import warnings
 from collections.abc import Hashable, Iterable, MutableMapping
-from typing import TYPE_CHECKING, Any, Callable, Literal, cast, overload
+from typing import TYPE_CHECKING, Any, Callable, Literal, Union, cast, overload
 
 import numpy as np
 import pandas as pd
@@ -29,7 +29,6 @@ from xarray.plot.utils import (
     _resolve_intervals_2dplot,
     _update_axes,
     get_axis,
-    import_matplotlib_pyplot,
     label_from_attrs,
 )
 
@@ -40,6 +39,7 @@ if TYPE_CHECKING:
     from matplotlib.container import BarContainer
     from matplotlib.contour import QuadContourSet
     from matplotlib.image import AxesImage
+    from matplotlib.patches import Polygon
     from mpl_toolkits.mplot3d.art3d import Line3D, Poly3DCollection
     from numpy.typing import ArrayLike
 
@@ -53,7 +53,7 @@ if TYPE_CHECKING:
     )
     from xarray.plot.facetgrid import FacetGrid
 
-_styles: MutableMapping[str, Any] = {
+_styles: dict[str, Any] = {
     # Add a white border to make it easier seeing overlapping markers:
     "scatter.edgecolors": "w",
 }
@@ -264,7 +264,9 @@ def plot(
     --------
     xarray.DataArray.squeeze
     """
-    darray = darray.squeeze().compute()
+    darray = darray.squeeze(
+        d for d, s in darray.sizes.items() if s == 1 and d not in (row, col, hue)
+    ).compute()
 
     plot_dims = set(darray.dims)
     plot_dims.discard(row)
@@ -307,7 +309,7 @@ def plot(
 
 
 @overload
-def line(  # type: ignore[misc]  # None is hashable :(
+def line(  # type: ignore[misc,unused-ignore]  # None is hashable :(
     darray: DataArray,
     *args: Any,
     row: None = None,  # no wrap -> primitive
@@ -325,8 +327,8 @@ def line(  # type: ignore[misc]  # None is hashable :(
     yscale: ScaleOptions = None,
     xticks: ArrayLike | None = None,
     yticks: ArrayLike | None = None,
-    xlim: ArrayLike | None = None,
-    ylim: ArrayLike | None = None,
+    xlim: tuple[float, float] | None = None,
+    ylim: tuple[float, float] | None = None,
     add_legend: bool = True,
     _labels: bool = True,
     **kwargs: Any,
@@ -336,7 +338,7 @@ def line(  # type: ignore[misc]  # None is hashable :(
 
 @overload
 def line(
-    darray,
+    darray: T_DataArray,
     *args: Any,
     row: Hashable,  # wrap -> FacetGrid
     col: Hashable | None = None,
@@ -353,18 +355,18 @@ def line(
     yscale: ScaleOptions = None,
     xticks: ArrayLike | None = None,
     yticks: ArrayLike | None = None,
-    xlim: ArrayLike | None = None,
-    ylim: ArrayLike | None = None,
+    xlim: tuple[float, float] | None = None,
+    ylim: tuple[float, float] | None = None,
     add_legend: bool = True,
     _labels: bool = True,
     **kwargs: Any,
-) -> FacetGrid[DataArray]:
+) -> FacetGrid[T_DataArray]:
     ...
 
 
 @overload
 def line(
-    darray,
+    darray: T_DataArray,
     *args: Any,
     row: Hashable | None = None,
     col: Hashable,  # wrap -> FacetGrid
@@ -381,19 +383,19 @@ def line(
     yscale: ScaleOptions = None,
     xticks: ArrayLike | None = None,
     yticks: ArrayLike | None = None,
-    xlim: ArrayLike | None = None,
-    ylim: ArrayLike | None = None,
+    xlim: tuple[float, float] | None = None,
+    ylim: tuple[float, float] | None = None,
     add_legend: bool = True,
     _labels: bool = True,
     **kwargs: Any,
-) -> FacetGrid[DataArray]:
+) -> FacetGrid[T_DataArray]:
     ...
 
 
 # This function signature should not change so that it can use
 # matplotlib format strings
 def line(
-    darray: DataArray,
+    darray: T_DataArray,
     *args: Any,
     row: Hashable | None = None,
     col: Hashable | None = None,
@@ -410,12 +412,12 @@ def line(
     yscale: ScaleOptions = None,
     xticks: ArrayLike | None = None,
     yticks: ArrayLike | None = None,
-    xlim: ArrayLike | None = None,
-    ylim: ArrayLike | None = None,
+    xlim: tuple[float, float] | None = None,
+    ylim: tuple[float, float] | None = None,
     add_legend: bool = True,
     _labels: bool = True,
     **kwargs: Any,
-) -> list[Line3D] | FacetGrid[DataArray]:
+) -> list[Line3D] | FacetGrid[T_DataArray]:
     """
     Line plot of DataArray values.
 
@@ -459,7 +461,7 @@ def line(
         Specifies scaling for the *x*- and *y*-axis, respectively.
     xticks, yticks : array-like, optional
         Specify tick locations for *x*- and *y*-axis.
-    xlim, ylim : array-like, optional
+    xlim, ylim : tuple[float, float], optional
         Specify *x*- and *y*-axis limits.
     add_legend : bool, default: True
         Add legend with *y* axis coordinates (2D inputs only).
@@ -538,7 +540,7 @@ def line(
 
 
 @overload
-def step(  # type: ignore[misc]  # None is hashable :(
+def step(  # type: ignore[misc,unused-ignore]  # None is hashable :(
     darray: DataArray,
     *args: Any,
     where: Literal["pre", "post", "mid"] = "pre",
@@ -654,10 +656,10 @@ def hist(
     yscale: ScaleOptions = None,
     xticks: ArrayLike | None = None,
     yticks: ArrayLike | None = None,
-    xlim: ArrayLike | None = None,
-    ylim: ArrayLike | None = None,
+    xlim: tuple[float, float] | None = None,
+    ylim: tuple[float, float] | None = None,
     **kwargs: Any,
-) -> tuple[np.ndarray, np.ndarray, BarContainer]:
+) -> tuple[np.ndarray, np.ndarray, BarContainer | Polygon]:
     """
     Histogram of DataArray.
 
@@ -691,7 +693,7 @@ def hist(
         Specifies scaling for the *x*- and *y*-axis, respectively.
     xticks, yticks : array-like, optional
         Specify tick locations for *x*- and *y*-axis.
-    xlim, ylim : array-like, optional
+    xlim, ylim : tuple[float, float], optional
         Specify *x*- and *y*-axis limits.
     **kwargs : optional
         Additional keyword arguments to :py:func:`matplotlib:matplotlib.pyplot.hist`.
@@ -708,14 +710,17 @@ def hist(
     no_nan = np.ravel(darray.to_numpy())
     no_nan = no_nan[pd.notnull(no_nan)]
 
-    primitive = ax.hist(no_nan, **kwargs)
+    n, bins, patches = cast(
+        tuple[np.ndarray, np.ndarray, Union["BarContainer", "Polygon"]],
+        ax.hist(no_nan, **kwargs),
+    )
 
     ax.set_title(darray._title_for_slice())
     ax.set_xlabel(label_from_attrs(darray))
 
     _update_axes(ax, xincrease, yincrease, xscale, yscale, xticks, yticks, xlim, ylim)
 
-    return primitive
+    return n, bins, patches
 
 
 def _plot1d(plotfunc):
@@ -779,9 +784,9 @@ def _plot1d(plotfunc):
         Specify tick locations for x-axes.
     yticks : ArrayLike or None, optional
         Specify tick locations for y-axes.
-    xlim : ArrayLike or None, optional
+    xlim : tuple[float, float] or None, optional
         Specify x-axes limits.
-    ylim : ArrayLike or None, optional
+    ylim : tuple[float, float] or None, optional
         Specify y-axes limits.
     cmap : matplotlib colormap name or colormap, optional
         The mapping from data values to color space. Either a
@@ -866,8 +871,8 @@ def _plot1d(plotfunc):
         yscale: ScaleOptions = None,
         xticks: ArrayLike | None = None,
         yticks: ArrayLike | None = None,
-        xlim: ArrayLike | None = None,
-        ylim: ArrayLike | None = None,
+        xlim: tuple[float, float] | None = None,
+        ylim: tuple[float, float] | None = None,
         cmap: str | Colormap | None = None,
         vmin: float | None = None,
         vmax: float | None = None,
@@ -879,7 +884,7 @@ def _plot1d(plotfunc):
         # All 1d plots in xarray share this function signature.
         # Method signature below should be consistent.
 
-        plt = import_matplotlib_pyplot()
+        import matplotlib.pyplot as plt
 
         if subplot_kws is None:
             subplot_kws = dict()
@@ -992,9 +997,13 @@ def _plot1d(plotfunc):
 
         with plt.rc_context(_styles):
             if z is not None:
+                import mpl_toolkits
+
                 if ax is None:
                     subplot_kws.update(projection="3d")
                 ax = get_axis(figsize, size, aspect, ax, **subplot_kws)
+                assert isinstance(ax, mpl_toolkits.mplot3d.axes3d.Axes3D)
+
                 # Using 30, 30 minimizes rotation of the plot. Making it easier to
                 # build on your intuition from 2D plots:
                 ax.view_init(azim=30, elev=30, vertical_axis="y")
@@ -1078,12 +1087,12 @@ def _plot1d(plotfunc):
 
 def _add_labels(
     add_labels: bool | Iterable[bool],
-    darrays: Iterable[DataArray],
+    darrays: Iterable[DataArray | None],
     suffixes: Iterable[str],
     rotate_labels: Iterable[bool],
     ax: Axes,
 ) -> None:
-    # Set x, y, z labels:
+    """Set x, y, z labels."""
     add_labels = [add_labels] * 3 if isinstance(add_labels, bool) else add_labels
     for axis, add_label, darray, suffix, rotate_label in zip(
         ("x", "y", "z"), add_labels, darrays, suffixes, rotate_labels
@@ -1107,7 +1116,7 @@ def _add_labels(
 
 
 @overload
-def scatter(
+def scatter(  # type: ignore[misc,unused-ignore]  # None is hashable :(
     darray: DataArray,
     *args: Any,
     x: Hashable | None = None,
@@ -1150,7 +1159,7 @@ def scatter(
 
 @overload
 def scatter(
-    darray: DataArray,
+    darray: T_DataArray,
     *args: Any,
     x: Hashable | None = None,
     y: Hashable | None = None,
@@ -1186,13 +1195,13 @@ def scatter(
     extend: ExtendOptions = None,
     levels: ArrayLike | None = None,
     **kwargs,
-) -> FacetGrid[DataArray]:
+) -> FacetGrid[T_DataArray]:
     ...
 
 
 @overload
 def scatter(
-    darray: DataArray,
+    darray: T_DataArray,
     *args: Any,
     x: Hashable | None = None,
     y: Hashable | None = None,
@@ -1228,7 +1237,7 @@ def scatter(
     extend: ExtendOptions = None,
     levels: ArrayLike | None = None,
     **kwargs,
-) -> FacetGrid[DataArray]:
+) -> FacetGrid[T_DataArray]:
     ...
 
 
@@ -1257,15 +1266,24 @@ def scatter(
     if sizeplt is not None:
         kwargs.update(s=sizeplt.to_numpy().ravel())
 
-    axis_order = ["x", "y", "z"]
+    plts_or_none = (xplt, yplt, zplt)
+    _add_labels(add_labels, plts_or_none, ("", "", ""), (True, False, False), ax)
 
-    plts_dict: dict[str, DataArray | None] = dict(x=xplt, y=yplt, z=zplt)
-    plts_or_none = [plts_dict[v] for v in axis_order]
-    plts = [p for p in plts_or_none if p is not None]
-    primitive = ax.scatter(*[p.to_numpy().ravel() for p in plts], **kwargs)
-    _add_labels(add_labels, plts, ("", "", ""), (True, False, False), ax)
+    xplt_np = None if xplt is None else xplt.to_numpy().ravel()
+    yplt_np = None if yplt is None else yplt.to_numpy().ravel()
+    zplt_np = None if zplt is None else zplt.to_numpy().ravel()
+    plts_np = tuple(p for p in (xplt_np, yplt_np, zplt_np) if p is not None)
 
-    return primitive
+    if len(plts_np) == 3:
+        import mpl_toolkits
+
+        assert isinstance(ax, mpl_toolkits.mplot3d.axes3d.Axes3D)
+        return ax.scatter(xplt_np, yplt_np, zplt_np, **kwargs)
+
+    if len(plts_np) == 2:
+        return ax.scatter(plts_np[0], plts_np[1], **kwargs)
+
+    raise ValueError("At least two variables required for a scatter plot.")
 
 
 def _plot2d(plotfunc):
@@ -1374,9 +1392,9 @@ def _plot2d(plotfunc):
         Specify tick locations for x-axes.
     yticks : ArrayLike or None, optional
         Specify tick locations for y-axes.
-    xlim : ArrayLike or None, optional
+    xlim : tuple[float, float] or None, optional
         Specify x-axes limits.
-    ylim : ArrayLike or None, optional
+    ylim : tuple[float, float] or None, optional
         Specify y-axes limits.
     norm : matplotlib.colors.Normalize, optional
         If ``norm`` has ``vmin`` or ``vmax`` specified, the corresponding
@@ -1429,8 +1447,8 @@ def _plot2d(plotfunc):
         yscale: ScaleOptions = None,
         xticks: ArrayLike | None = None,
         yticks: ArrayLike | None = None,
-        xlim: ArrayLike | None = None,
-        ylim: ArrayLike | None = None,
+        xlim: tuple[float, float] | None = None,
+        ylim: tuple[float, float] | None = None,
         norm: Normalize | None = None,
         **kwargs: Any,
     ) -> Any:
@@ -1502,8 +1520,6 @@ def _plot2d(plotfunc):
             # TypeError to be consistent with pandas
             raise TypeError("No numeric data to plot.")
 
-        plt = import_matplotlib_pyplot()
-
         if (
             plotfunc.__name__ == "surface"
             and not kwargs.get("_is_facetgrid", False)
@@ -1616,6 +1632,9 @@ def _plot2d(plotfunc):
             ax.set_ylabel(label_from_attrs(darray[ylab], ylab_extra))
             ax.set_title(darray._title_for_slice())
             if plotfunc.__name__ == "surface":
+                import mpl_toolkits
+
+                assert isinstance(ax, mpl_toolkits.mplot3d.axes3d.Axes3D)
                 ax.set_zlabel(label_from_attrs(darray))
 
         if add_colorbar:
@@ -1656,7 +1675,7 @@ def _plot2d(plotfunc):
 
 
 @overload
-def imshow(
+def imshow(  # type: ignore[misc,unused-ignore]  # None is hashable :(
     darray: DataArray,
     x: Hashable | None = None,
     y: Hashable | None = None,
@@ -1698,7 +1717,7 @@ def imshow(
 
 @overload
 def imshow(
-    darray: DataArray,
+    darray: T_DataArray,
     x: Hashable | None = None,
     y: Hashable | None = None,
     *,
@@ -1733,13 +1752,13 @@ def imshow(
     ylim: ArrayLike | None = None,
     norm: Normalize | None = None,
     **kwargs: Any,
-) -> FacetGrid[DataArray]:
+) -> FacetGrid[T_DataArray]:
     ...
 
 
 @overload
 def imshow(
-    darray: DataArray,
+    darray: T_DataArray,
     x: Hashable | None = None,
     y: Hashable | None = None,
     *,
@@ -1774,7 +1793,7 @@ def imshow(
     ylim: ArrayLike | None = None,
     norm: Normalize | None = None,
     **kwargs: Any,
-) -> FacetGrid[DataArray]:
+) -> FacetGrid[T_DataArray]:
     ...
 
 
@@ -1875,7 +1894,7 @@ def imshow(
 
 
 @overload
-def contour(
+def contour(  # type: ignore[misc,unused-ignore]  # None is hashable :(
     darray: DataArray,
     x: Hashable | None = None,
     y: Hashable | None = None,
@@ -1917,7 +1936,7 @@ def contour(
 
 @overload
 def contour(
-    darray: DataArray,
+    darray: T_DataArray,
     x: Hashable | None = None,
     y: Hashable | None = None,
     *,
@@ -1952,13 +1971,13 @@ def contour(
     ylim: ArrayLike | None = None,
     norm: Normalize | None = None,
     **kwargs: Any,
-) -> FacetGrid[DataArray]:
+) -> FacetGrid[T_DataArray]:
     ...
 
 
 @overload
 def contour(
-    darray: DataArray,
+    darray: T_DataArray,
     x: Hashable | None = None,
     y: Hashable | None = None,
     *,
@@ -1993,7 +2012,7 @@ def contour(
     ylim: ArrayLike | None = None,
     norm: Normalize | None = None,
     **kwargs: Any,
-) -> FacetGrid[DataArray]:
+) -> FacetGrid[T_DataArray]:
     ...
 
 
@@ -2011,7 +2030,7 @@ def contour(
 
 
 @overload
-def contourf(
+def contourf(  # type: ignore[misc,unused-ignore]  # None is hashable :(
     darray: DataArray,
     x: Hashable | None = None,
     y: Hashable | None = None,
@@ -2053,7 +2072,7 @@ def contourf(
 
 @overload
 def contourf(
-    darray: DataArray,
+    darray: T_DataArray,
     x: Hashable | None = None,
     y: Hashable | None = None,
     *,
@@ -2088,13 +2107,13 @@ def contourf(
     ylim: ArrayLike | None = None,
     norm: Normalize | None = None,
     **kwargs: Any,
-) -> FacetGrid[DataArray]:
+) -> FacetGrid[T_DataArray]:
     ...
 
 
 @overload
 def contourf(
-    darray: DataArray,
+    darray: T_DataArray,
     x: Hashable | None = None,
     y: Hashable | None = None,
     *,
@@ -2129,7 +2148,7 @@ def contourf(
     ylim: ArrayLike | None = None,
     norm: Normalize | None = None,
     **kwargs: Any,
-) -> FacetGrid[DataArray]:
+) -> FacetGrid[T_DataArray]:
     ...
 
 
@@ -2147,7 +2166,7 @@ def contourf(
 
 
 @overload
-def pcolormesh(
+def pcolormesh(  # type: ignore[misc,unused-ignore]  # None is hashable :(
     darray: DataArray,
     x: Hashable | None = None,
     y: Hashable | None = None,
@@ -2189,7 +2208,7 @@ def pcolormesh(
 
 @overload
 def pcolormesh(
-    darray: DataArray,
+    darray: T_DataArray,
     x: Hashable | None = None,
     y: Hashable | None = None,
     *,
@@ -2224,13 +2243,13 @@ def pcolormesh(
     ylim: ArrayLike | None = None,
     norm: Normalize | None = None,
     **kwargs: Any,
-) -> FacetGrid[DataArray]:
+) -> FacetGrid[T_DataArray]:
     ...
 
 
 @overload
 def pcolormesh(
-    darray: DataArray,
+    darray: T_DataArray,
     x: Hashable | None = None,
     y: Hashable | None = None,
     *,
@@ -2265,7 +2284,7 @@ def pcolormesh(
     ylim: ArrayLike | None = None,
     norm: Normalize | None = None,
     **kwargs: Any,
-) -> FacetGrid[DataArray]:
+) -> FacetGrid[T_DataArray]:
     ...
 
 
@@ -2376,7 +2395,7 @@ def surface(
 
 @overload
 def surface(
-    darray: DataArray,
+    darray: T_DataArray,
     x: Hashable | None = None,
     y: Hashable | None = None,
     *,
@@ -2411,13 +2430,13 @@ def surface(
     ylim: ArrayLike | None = None,
     norm: Normalize | None = None,
     **kwargs: Any,
-) -> FacetGrid[DataArray]:
+) -> FacetGrid[T_DataArray]:
     ...
 
 
 @overload
 def surface(
-    darray: DataArray,
+    darray: T_DataArray,
     x: Hashable | None = None,
     y: Hashable | None = None,
     *,
@@ -2452,7 +2471,7 @@ def surface(
     ylim: ArrayLike | None = None,
     norm: Normalize | None = None,
     **kwargs: Any,
-) -> FacetGrid[DataArray]:
+) -> FacetGrid[T_DataArray]:
     ...
 
 
@@ -2465,5 +2484,8 @@ def surface(
 
     Wraps :py:meth:`matplotlib:mpl_toolkits.mplot3d.axes3d.Axes3D.plot_surface`.
     """
+    import mpl_toolkits
+
+    assert isinstance(ax, mpl_toolkits.mplot3d.axes3d.Axes3D)
     primitive = ax.plot_surface(x, y, z, **kwargs)
     return primitive
diff --git a/xarray/plot/dataset_plot.py b/xarray/plot/dataset_plot.py
index edb32b73..a3ca201e 100644
--- a/xarray/plot/dataset_plot.py
+++ b/xarray/plot/dataset_plot.py
@@ -321,7 +321,7 @@ def _dsplot(plotfunc):
 
 
 @overload
-def quiver(
+def quiver(  # type: ignore[misc,unused-ignore]  # None is hashable :(
     ds: Dataset,
     *args: Any,
     x: Hashable | None = None,
@@ -475,7 +475,7 @@ def quiver(
 
 
 @overload
-def streamplot(
+def streamplot(  # type: ignore[misc,unused-ignore]  # None is hashable :(
     ds: Dataset,
     *args: Any,
     x: Hashable | None = None,
@@ -632,7 +632,6 @@ def streamplot(
         du = du.transpose(ydim, xdim)
         dv = dv.transpose(ydim, xdim)
 
-    args = [dx.values, dy.values, du.values, dv.values]
     hue = kwargs.pop("hue")
     cmap_params = kwargs.pop("cmap_params")
 
@@ -646,7 +645,9 @@ def streamplot(
             )
 
     kwargs.pop("hue_style")
-    hdl = ax.streamplot(*args, **kwargs, **cmap_params)
+    hdl = ax.streamplot(
+        dx.values, dy.values, du.values, dv.values, **kwargs, **cmap_params
+    )
 
     # Return .lines so colorbar creation works properly
     return hdl.lines
@@ -748,7 +749,7 @@ def _temp_dataarray(ds: Dataset, y: Hashable, locals_: dict[str, Any]) -> DataAr
 
 
 @overload
-def scatter(
+def scatter(  # type: ignore[misc,unused-ignore]  # None is hashable :(
     ds: Dataset,
     *args: Any,
     x: Hashable | None = None,
diff --git a/xarray/plot/facetgrid.py b/xarray/plot/facetgrid.py
index 2b348d8b..faf809a8 100644
--- a/xarray/plot/facetgrid.py
+++ b/xarray/plot/facetgrid.py
@@ -9,7 +9,7 @@ from typing import TYPE_CHECKING, Any, Callable, Generic, Literal, TypeVar, cast
 import numpy as np
 
 from xarray.core.formatting import format_item
-from xarray.core.types import HueStyleOptions, T_Xarray
+from xarray.core.types import HueStyleOptions, T_DataArrayOrSet
 from xarray.plot.utils import (
     _LINEWIDTH_RANGE,
     _MARKERSIZE_RANGE,
@@ -21,7 +21,6 @@ from xarray.plot.utils import (
     _Normalize,
     _parse_size,
     _process_cmap_cbar_kwargs,
-    import_matplotlib_pyplot,
     label_from_attrs,
 )
 
@@ -60,7 +59,7 @@ def _nicetitle(coord, value, maxchar, template):
 T_FacetGrid = TypeVar("T_FacetGrid", bound="FacetGrid")
 
 
-class FacetGrid(Generic[T_Xarray]):
+class FacetGrid(Generic[T_DataArrayOrSet]):
     """
     Initialize the Matplotlib figure and FacetGrid object.
 
@@ -101,7 +100,7 @@ class FacetGrid(Generic[T_Xarray]):
         sometimes the rightmost grid positions in the bottom row.
     """
 
-    data: T_Xarray
+    data: T_DataArrayOrSet
     name_dicts: np.ndarray
     fig: Figure
     axs: np.ndarray
@@ -126,7 +125,7 @@ class FacetGrid(Generic[T_Xarray]):
 
     def __init__(
         self,
-        data: T_Xarray,
+        data: T_DataArrayOrSet,
         col: Hashable | None = None,
         row: Hashable | None = None,
         col_wrap: int | None = None,
@@ -166,7 +165,7 @@ class FacetGrid(Generic[T_Xarray]):
 
         """
 
-        plt = import_matplotlib_pyplot()
+        import matplotlib.pyplot as plt
 
         # Handle corner case of nonunique coordinates
         rep_col = col is not None and not data[col].to_index().is_unique
@@ -681,7 +680,10 @@ class FacetGrid(Generic[T_Xarray]):
 
     def _adjust_fig_for_guide(self, guide) -> None:
         # Draw the plot to set the bounding boxes correctly
-        renderer = self.fig.canvas.get_renderer()
+        if hasattr(self.fig.canvas, "get_renderer"):
+            renderer = self.fig.canvas.get_renderer()
+        else:
+            raise RuntimeError("MPL backend has no renderer")
         self.fig.draw(renderer)
 
         # Calculate and set the new width of the figure so the legend fits
@@ -731,6 +733,9 @@ class FacetGrid(Generic[T_Xarray]):
         if hasattr(self._mappables[-1], "extend"):
             kwargs.pop("extend", None)
         if "label" not in kwargs:
+            from xarray import DataArray
+
+            assert isinstance(self.data, DataArray)
             kwargs.setdefault("label", label_from_attrs(self.data))
         self.cbar = self.fig.colorbar(
             self._mappables[-1], ax=list(self.axs.flat), **kwargs
@@ -985,7 +990,7 @@ class FacetGrid(Generic[T_Xarray]):
         self : FacetGrid object
 
         """
-        plt = import_matplotlib_pyplot()
+        import matplotlib.pyplot as plt
 
         for ax, namedict in zip(self.axs.flat, self.name_dicts.flat):
             if namedict is not None:
@@ -1004,7 +1009,7 @@ class FacetGrid(Generic[T_Xarray]):
 
 
 def _easy_facetgrid(
-    data: T_Xarray,
+    data: T_DataArrayOrSet,
     plotfunc: Callable,
     kind: Literal["line", "dataarray", "dataset", "plot1d"],
     x: Hashable | None = None,
@@ -1020,7 +1025,7 @@ def _easy_facetgrid(
     ax: Axes | None = None,
     figsize: Iterable[float] | None = None,
     **kwargs: Any,
-) -> FacetGrid[T_Xarray]:
+) -> FacetGrid[T_DataArrayOrSet]:
     """
     Convenience method to call xarray.plot.FacetGrid from 2d plotting methods
 
diff --git a/xarray/plot/utils.py b/xarray/plot/utils.py
index b8cc4ff7..5694acc0 100644
--- a/xarray/plot/utils.py
+++ b/xarray/plot/utils.py
@@ -47,14 +47,6 @@ _MARKERSIZE_RANGE = (18.0, 36.0, 72.0)
 _LINEWIDTH_RANGE = (1.5, 1.5, 6.0)
 
 
-def import_matplotlib_pyplot():
-    """import pyplot"""
-    # TODO: This function doesn't do anything (after #6109), remove it?
-    import matplotlib.pyplot as plt
-
-    return plt
-
-
 def _determine_extend(calc_data, vmin, vmax):
     extend_min = calc_data.min() < vmin
     extend_max = calc_data.max() > vmax
@@ -505,28 +497,29 @@ def _maybe_gca(**subplot_kws: Any) -> Axes:
     return plt.axes(**subplot_kws)
 
 
-def _get_units_from_attrs(da) -> str:
+def _get_units_from_attrs(da: DataArray) -> str:
     """Extracts and formats the unit/units from a attributes."""
     pint_array_type = DuckArrayModule("pint").type
     units = " [{}]"
     if isinstance(da.data, pint_array_type):
-        units = units.format(str(da.data.units))
-    elif da.attrs.get("units"):
-        units = units.format(da.attrs["units"])
-    elif da.attrs.get("unit"):
-        units = units.format(da.attrs["unit"])
-    else:
-        units = ""
-    return units
+        return units.format(str(da.data.units))
+    if "units" in da.attrs:
+        return units.format(da.attrs["units"])
+    if "unit" in da.attrs:
+        return units.format(da.attrs["unit"])
+    return ""
 
 
-def label_from_attrs(da, extra: str = "") -> str:
+def label_from_attrs(da: DataArray | None, extra: str = "") -> str:
     """Makes informative labels if variable metadata (attrs) follows
     CF conventions."""
+    if da is None:
+        return ""
+
     name: str = "{}"
-    if da.attrs.get("long_name"):
+    if "long_name" in da.attrs:
         name = name.format(da.attrs["long_name"])
-    elif da.attrs.get("standard_name"):
+    elif "standard_name" in da.attrs:
         name = name.format(da.attrs["standard_name"])
     elif da.name is not None:
         name = name.format(da.name)
@@ -774,8 +767,8 @@ def _update_axes(
     yscale: ScaleOptions = None,
     xticks: ArrayLike | None = None,
     yticks: ArrayLike | None = None,
-    xlim: ArrayLike | None = None,
-    ylim: ArrayLike | None = None,
+    xlim: tuple[float, float] | None = None,
+    ylim: tuple[float, float] | None = None,
 ) -> None:
     """
     Update axes with provided parameters
@@ -1166,7 +1159,7 @@ def legend_elements(
 
 def _legend_add_subtitle(handles, labels, text):
     """Add a subtitle to legend handles."""
-    plt = import_matplotlib_pyplot()
+    import matplotlib.pyplot as plt
 
     if text and len(handles) > 1:
         # Create a blank handle that's not visible, the
@@ -1184,7 +1177,7 @@ def _legend_add_subtitle(handles, labels, text):
 
 def _adjust_legend_subtitles(legend):
     """Make invisible-handle "subtitles" entries look more like titles."""
-    plt = import_matplotlib_pyplot()
+    import matplotlib.pyplot as plt
 
     # Legend title not in rcParams until 3.0
     font_size = plt.rcParams.get("legend.title_fontsize", None)
@@ -1640,7 +1633,7 @@ class _Normalize(Sequence):
         >>> aa.format(1)
         '3.0'
         """
-        plt = import_matplotlib_pyplot()
+        import matplotlib.pyplot as plt
 
         def _func(x: Any, pos: None | Any = None):
             return f"{self._lookup_arr([x])[0]}"


Test patch:
diff --git a/xarray/tests/test_backends.py b/xarray/tests/test_backends.py
index ea088b36..5bd51709 100644
--- a/xarray/tests/test_backends.py
+++ b/xarray/tests/test_backends.py
@@ -54,7 +54,6 @@ from xarray.core import indexing
 from xarray.core.options import set_options
 from xarray.core.pycompat import array_type
 from xarray.tests import (
-    arm_xfail,
     assert_allclose,
     assert_array_equal,
     assert_equal,
@@ -526,7 +525,6 @@ class DatasetIOBase:
             assert_identical(expected, actual)
             assert actual["x"].encoding["_Encoding"] == "ascii"
 
-    @arm_xfail
     def test_roundtrip_numpy_datetime_data(self) -> None:
         times = pd.to_datetime(["2000-01-01", "2000-01-02", "NaT"])
         expected = Dataset({"t": ("t", times), "t0": times[0]})
diff --git a/xarray/tests/test_coding_times.py b/xarray/tests/test_coding_times.py
index 580de878..079e432b 100644
--- a/xarray/tests/test_coding_times.py
+++ b/xarray/tests/test_coding_times.py
@@ -20,6 +20,7 @@ from xarray import (
 )
 from xarray.coding.times import (
     _encode_datetime_with_cftime,
+    _numpy_to_netcdf_timeunit,
     _should_cftime_be_used,
     cftime_to_nptime,
     decode_cf_datetime,
@@ -29,7 +30,7 @@ from xarray.coding.times import (
 from xarray.coding.variables import SerializationWarning
 from xarray.conventions import _update_bounds_attributes, cf_encoder
 from xarray.core.common import contains_cftime_datetimes
-from xarray.testing import assert_equal, assert_identical
+from xarray.testing import assert_allclose, assert_equal, assert_identical
 from xarray.tests import (
     FirstElementAccessibleArray,
     arm_xfail,
@@ -110,6 +111,7 @@ def _all_cftime_date_types():
 
 @requires_cftime
 @pytest.mark.filterwarnings("ignore:Ambiguous reference date string")
+@pytest.mark.filterwarnings("ignore:Times can't be serialized faithfully")
 @pytest.mark.parametrize(["num_dates", "units", "calendar"], _CF_DATETIME_TESTS)
 def test_cf_datetime(num_dates, units, calendar) -> None:
     import cftime
@@ -567,6 +569,7 @@ def test_infer_cftime_datetime_units(calendar, date_args, expected) -> None:
     assert expected == coding.times.infer_datetime_units(dates)
 
 
+@pytest.mark.filterwarnings("ignore:Timedeltas can't be serialized faithfully")
 @pytest.mark.parametrize(
     ["timedeltas", "units", "numbers"],
     [
@@ -576,10 +579,10 @@ def test_infer_cftime_datetime_units(calendar, date_args, expected) -> None:
         ("1ms", "milliseconds", np.int64(1)),
         ("1us", "microseconds", np.int64(1)),
         ("1ns", "nanoseconds", np.int64(1)),
-        (["NaT", "0s", "1s"], None, [np.nan, 0, 1]),
+        (["NaT", "0s", "1s"], None, [np.iinfo(np.int64).min, 0, 1]),
         (["30m", "60m"], "hours", [0.5, 1.0]),
-        ("NaT", "days", np.nan),
-        (["NaT", "NaT"], "days", [np.nan, np.nan]),
+        ("NaT", "days", np.iinfo(np.int64).min),
+        (["NaT", "NaT"], "days", [np.iinfo(np.int64).min, np.iinfo(np.int64).min]),
     ],
 )
 def test_cf_timedelta(timedeltas, units, numbers) -> None:
@@ -1020,6 +1023,7 @@ def test_decode_ambiguous_time_warns(calendar) -> None:
     np.testing.assert_array_equal(result, expected)
 
 
+@pytest.mark.filterwarnings("ignore:Times can't be serialized faithfully")
 @pytest.mark.parametrize("encoding_units", FREQUENCIES_TO_ENCODING_UNITS.values())
 @pytest.mark.parametrize("freq", FREQUENCIES_TO_ENCODING_UNITS.keys())
 @pytest.mark.parametrize("date_range", [pd.date_range, cftime_range])
@@ -1191,3 +1195,147 @@ def test_contains_cftime_lazy() -> None:
     )
     array = FirstElementAccessibleArray(times)
     assert _contains_cftime_datetimes(array)
+
+
+@pytest.mark.parametrize(
+    "timestr, timeunit, dtype, fill_value, use_encoding",
+    [
+        ("1677-09-21T00:12:43.145224193", "ns", np.int64, 20, True),
+        ("1970-09-21T00:12:44.145224808", "ns", np.float64, 1e30, True),
+        (
+            "1677-09-21T00:12:43.145225216",
+            "ns",
+            np.float64,
+            -9.223372036854776e18,
+            True,
+        ),
+        ("1677-09-21T00:12:43.145224193", "ns", np.int64, None, False),
+        ("1677-09-21T00:12:43.145225", "us", np.int64, None, False),
+        ("1970-01-01T00:00:01.000001", "us", np.int64, None, False),
+    ],
+)
+def test_roundtrip_datetime64_nanosecond_precision(
+    timestr: str,
+    timeunit: str,
+    dtype: np.typing.DTypeLike,
+    fill_value: int | float | None,
+    use_encoding: bool,
+) -> None:
+    # test for GH7817
+    time = np.datetime64(timestr, timeunit)
+    times = [np.datetime64("1970-01-01T00:00:00", timeunit), np.datetime64("NaT"), time]
+
+    if use_encoding:
+        encoding = dict(dtype=dtype, _FillValue=fill_value)
+    else:
+        encoding = {}
+
+    var = Variable(["time"], times, encoding=encoding)
+    assert var.dtype == np.dtype("<M8[ns]")
+
+    encoded_var = conventions.encode_cf_variable(var)
+    assert (
+        encoded_var.attrs["units"]
+        == f"{_numpy_to_netcdf_timeunit(timeunit)} since 1970-01-01 00:00:00"
+    )
+    assert encoded_var.attrs["calendar"] == "proleptic_gregorian"
+    assert encoded_var.data.dtype == dtype
+
+    decoded_var = conventions.decode_cf_variable("foo", encoded_var)
+    assert decoded_var.dtype == np.dtype("<M8[ns]")
+    assert (
+        decoded_var.encoding["units"]
+        == f"{_numpy_to_netcdf_timeunit(timeunit)} since 1970-01-01 00:00:00"
+    )
+    assert decoded_var.encoding["dtype"] == dtype
+    assert decoded_var.encoding["calendar"] == "proleptic_gregorian"
+    assert_identical(var, decoded_var)
+
+
+def test_roundtrip_datetime64_nanosecond_precision_warning() -> None:
+    # test warning if times can't be serialized faithfully
+    times = [
+        np.datetime64("1970-01-01T00:01:00", "ns"),
+        np.datetime64("NaT"),
+        np.datetime64("1970-01-02T00:01:00", "ns"),
+    ]
+    units = "days since 1970-01-10T01:01:00"
+    needed_units = "hours"
+    encoding = dict(_FillValue=20, units=units)
+    var = Variable(["time"], times, encoding=encoding)
+    wmsg = (
+        f"Times can't be serialized faithfully with requested units {units!r}. "
+        f"Resolution of {needed_units!r} needed. "
+    )
+    with pytest.warns(UserWarning, match=wmsg):
+        encoded_var = conventions.encode_cf_variable(var)
+
+    decoded_var = conventions.decode_cf_variable("foo", encoded_var)
+    assert_identical(var, decoded_var)
+
+
+@pytest.mark.parametrize(
+    "dtype, fill_value",
+    [(np.int64, 20), (np.int64, np.iinfo(np.int64).min), (np.float64, 1e30)],
+)
+def test_roundtrip_timedelta64_nanosecond_precision(
+    dtype: np.typing.DTypeLike, fill_value: int | float
+) -> None:
+    # test for GH7942
+    one_day = np.timedelta64(1, "ns")
+    nat = np.timedelta64("nat", "ns")
+    timedelta_values = (np.arange(5) * one_day).astype("timedelta64[ns]")
+    timedelta_values[2] = nat
+    timedelta_values[4] = nat
+
+    encoding = dict(dtype=dtype, _FillValue=fill_value)
+    var = Variable(["time"], timedelta_values, encoding=encoding)
+
+    encoded_var = conventions.encode_cf_variable(var)
+    decoded_var = conventions.decode_cf_variable("foo", encoded_var)
+
+    assert_identical(var, decoded_var)
+
+
+def test_roundtrip_timedelta64_nanosecond_precision_warning() -> None:
+    # test warning if timedeltas can't be serialized faithfully
+    one_day = np.timedelta64(1, "D")
+    nat = np.timedelta64("nat", "ns")
+    timedelta_values = (np.arange(5) * one_day).astype("timedelta64[ns]")
+    timedelta_values[2] = nat
+    timedelta_values[4] = np.timedelta64(12, "h").astype("timedelta64[ns]")
+
+    units = "days"
+    needed_units = "hours"
+    wmsg = (
+        f"Timedeltas can't be serialized faithfully with requested units {units!r}. "
+        f"Resolution of {needed_units!r} needed. "
+    )
+    encoding = dict(_FillValue=20, units=units)
+    var = Variable(["time"], timedelta_values, encoding=encoding)
+    with pytest.warns(UserWarning, match=wmsg):
+        encoded_var = conventions.encode_cf_variable(var)
+    decoded_var = conventions.decode_cf_variable("foo", encoded_var)
+    assert_allclose(var, decoded_var)
+
+
+def test_roundtrip_float_times() -> None:
+    fill_value = 20.0
+    times = [np.datetime64("2000-01-01 12:00:00", "ns"), np.datetime64("NaT", "ns")]
+
+    units = "days since 2000-01-01"
+    var = Variable(
+        ["time"],
+        times,
+        encoding=dict(dtype=np.float64, _FillValue=fill_value, units=units),
+    )
+
+    encoded_var = conventions.encode_cf_variable(var)
+    np.testing.assert_array_equal(encoded_var, np.array([0.5, 20.0]))
+    assert encoded_var.attrs["units"] == units
+    assert encoded_var.attrs["_FillValue"] == fill_value
+
+    decoded_var = conventions.decode_cf_variable("foo", encoded_var)
+    assert_identical(var, decoded_var)
+    assert decoded_var.encoding["units"] == units
+    assert decoded_var.encoding["_FillValue"] == fill_value
diff --git a/xarray/tests/test_computation.py b/xarray/tests/test_computation.py
index 372998a6..052672ef 100644
--- a/xarray/tests/test_computation.py
+++ b/xarray/tests/test_computation.py
@@ -257,6 +257,168 @@ def test_apply_two_outputs() -> None:
     assert_identical(out1, dataset)
 
 
+def test_apply_missing_dims() -> None:
+    ## Single arg
+
+    def add_one(a, core_dims, on_missing_core_dim):
+        return apply_ufunc(
+            lambda x: x + 1,
+            a,
+            input_core_dims=core_dims,
+            output_core_dims=core_dims,
+            on_missing_core_dim=on_missing_core_dim,
+        )
+
+    array = np.arange(6).reshape(2, 3)
+    variable = xr.Variable(["x", "y"], array)
+    variable_no_y = xr.Variable(["x", "z"], array)
+
+    ds = xr.Dataset({"x_y": variable, "x_z": variable_no_y})
+
+    # Check the standard stuff works OK
+    assert_identical(
+        add_one(ds[["x_y"]], core_dims=[["y"]], on_missing_core_dim="raise"),
+        ds[["x_y"]] + 1,
+    )
+
+    # `raise` — should raise on a missing dim
+    with pytest.raises(ValueError):
+        add_one(ds, core_dims=[["y"]], on_missing_core_dim="raise"),
+
+    # `drop` — should drop the var with the missing dim
+    assert_identical(
+        add_one(ds, core_dims=[["y"]], on_missing_core_dim="drop"),
+        (ds + 1).drop_vars("x_z"),
+    )
+
+    # `copy` — should not add one to the missing with `copy`
+    copy_result = add_one(ds, core_dims=[["y"]], on_missing_core_dim="copy")
+    assert_identical(copy_result["x_y"], (ds + 1)["x_y"])
+    assert_identical(copy_result["x_z"], ds["x_z"])
+
+    ## Multiple args
+
+    def sum_add(a, b, core_dims, on_missing_core_dim):
+        return apply_ufunc(
+            lambda a, b, axis=None: a.sum(axis) + b.sum(axis),
+            a,
+            b,
+            input_core_dims=core_dims,
+            on_missing_core_dim=on_missing_core_dim,
+        )
+
+    # Check the standard stuff works OK
+    assert_identical(
+        sum_add(
+            ds[["x_y"]],
+            ds[["x_y"]],
+            core_dims=[["x", "y"], ["x", "y"]],
+            on_missing_core_dim="raise",
+        ),
+        ds[["x_y"]].sum() * 2,
+    )
+
+    # `raise` — should raise on a missing dim
+    with pytest.raises(
+        ValueError,
+        match=r".*Missing core dims \{'y'\} from arg number 1 on a variable named `x_z`:\n.*<xarray.Variable \(x: 2, z: ",
+    ):
+        sum_add(
+            ds[["x_z"]],
+            ds[["x_z"]],
+            core_dims=[["x", "y"], ["x", "z"]],
+            on_missing_core_dim="raise",
+        )
+
+    # `raise` on a missing dim on a non-first arg
+    with pytest.raises(
+        ValueError,
+        match=r".*Missing core dims \{'y'\} from arg number 2 on a variable named `x_z`:\n.*<xarray.Variable \(x: 2, z: ",
+    ):
+        sum_add(
+            ds[["x_z"]],
+            ds[["x_z"]],
+            core_dims=[["x", "z"], ["x", "y"]],
+            on_missing_core_dim="raise",
+        )
+
+    # `drop` — should drop the var with the missing dim
+    assert_identical(
+        sum_add(
+            ds[["x_z"]],
+            ds[["x_z"]],
+            core_dims=[["x", "y"], ["x", "y"]],
+            on_missing_core_dim="drop",
+        ),
+        ds[[]],
+    )
+
+    # `copy` — should drop the var with the missing dim
+    assert_identical(
+        sum_add(
+            ds[["x_z"]],
+            ds[["x_z"]],
+            core_dims=[["x", "y"], ["x", "y"]],
+            on_missing_core_dim="copy",
+        ),
+        ds[["x_z"]],
+    )
+
+    ## Multiple vars per arg
+    assert_identical(
+        sum_add(
+            ds[["x_y", "x_y"]],
+            ds[["x_y", "x_y"]],
+            core_dims=[["x", "y"], ["x", "y"]],
+            on_missing_core_dim="raise",
+        ),
+        ds[["x_y", "x_y"]].sum() * 2,
+    )
+
+    assert_identical(
+        sum_add(
+            ds,
+            ds,
+            core_dims=[["x", "y"], ["x", "y"]],
+            on_missing_core_dim="drop",
+        ),
+        ds[["x_y"]].sum() * 2,
+    )
+
+    assert_identical(
+        sum_add(
+            # The first one has the wrong dims — using `z` for the `x_t` var
+            ds.assign(x_t=ds["x_z"])[["x_y", "x_t"]],
+            ds.assign(x_t=ds["x_y"])[["x_y", "x_t"]],
+            core_dims=[["x", "y"], ["x", "y"]],
+            on_missing_core_dim="drop",
+        ),
+        ds[["x_y"]].sum() * 2,
+    )
+
+    assert_identical(
+        sum_add(
+            ds.assign(x_t=ds["x_y"])[["x_y", "x_t"]],
+            # The second one has the wrong dims — using `z` for the `x_t` var (seems
+            # duplicative but this was a bug in the initial impl...)
+            ds.assign(x_t=ds["x_z"])[["x_y", "x_t"]],
+            core_dims=[["x", "y"], ["x", "y"]],
+            on_missing_core_dim="drop",
+        ),
+        ds[["x_y"]].sum() * 2,
+    )
+
+    assert_identical(
+        sum_add(
+            ds.assign(x_t=ds["x_y"])[["x_y", "x_t"]],
+            ds.assign(x_t=ds["x_z"])[["x_y", "x_t"]],
+            core_dims=[["x", "y"], ["x", "y"]],
+            on_missing_core_dim="copy",
+        ),
+        ds.drop_vars("x_z").assign(x_y=30, x_t=ds["x_y"]),
+    )
+
+
 @requires_dask
 def test_apply_dask_parallelized_two_outputs() -> None:
     data_array = xr.DataArray([[0, 1, 2], [1, 2, 3]], dims=("x", "y"))
diff --git a/xarray/tests/test_plot.py b/xarray/tests/test_plot.py
index a6c15a86..6cc061a7 100644
--- a/xarray/tests/test_plot.py
+++ b/xarray/tests/test_plot.py
@@ -43,6 +43,7 @@ has_nc_time_axis = module_available("nc_time_axis")
 # import mpl and change the backend before other mpl imports
 try:
     import matplotlib as mpl
+    import matplotlib.dates
     import matplotlib.pyplot as plt
     import mpl_toolkits
 except ImportError:
@@ -428,6 +429,7 @@ class TestPlot(PlotTestCase):
         ]:
             p = a.plot.pcolormesh(x=x, y=y)
             v = p.get_paths()[0].vertices
+            assert isinstance(v, np.ndarray)
 
             # Check all vertices are different, except last vertex which should be the
             # same as the first
@@ -447,7 +449,7 @@ class TestPlot(PlotTestCase):
     def test_contourf_cmap_set(self) -> None:
         a = DataArray(easy_array((4, 4)), dims=["z", "time"])
 
-        cmap = mpl.cm.viridis
+        cmap_expected = mpl.colormaps["viridis"]
 
         # use copy to ensure cmap is not changed by contourf()
         # Set vmin and vmax so that _build_discrete_colormap is called with
@@ -457,55 +459,59 @@ class TestPlot(PlotTestCase):
         # extend='neither' (but if extend='neither' the under and over values
         # would not be used because the data would all be within the plotted
         # range)
-        pl = a.plot.contourf(cmap=copy(cmap), vmin=0.1, vmax=0.9)
+        pl = a.plot.contourf(cmap=copy(cmap_expected), vmin=0.1, vmax=0.9)
 
         # check the set_bad color
+        cmap = pl.cmap
+        assert cmap is not None
         assert_array_equal(
-            pl.cmap(np.ma.masked_invalid([np.nan]))[0],
             cmap(np.ma.masked_invalid([np.nan]))[0],
+            cmap_expected(np.ma.masked_invalid([np.nan]))[0],
         )
 
         # check the set_under color
-        assert pl.cmap(-np.inf) == cmap(-np.inf)
+        assert cmap(-np.inf) == cmap_expected(-np.inf)
 
         # check the set_over color
-        assert pl.cmap(np.inf) == cmap(np.inf)
+        assert cmap(np.inf) == cmap_expected(np.inf)
 
     def test_contourf_cmap_set_with_bad_under_over(self) -> None:
         a = DataArray(easy_array((4, 4)), dims=["z", "time"])
 
         # make a copy here because we want a local cmap that we will modify.
-        cmap = copy(mpl.cm.viridis)
+        cmap_expected = copy(mpl.colormaps["viridis"])
 
-        cmap.set_bad("w")
+        cmap_expected.set_bad("w")
         # check we actually changed the set_bad color
         assert np.all(
-            cmap(np.ma.masked_invalid([np.nan]))[0]
-            != mpl.cm.viridis(np.ma.masked_invalid([np.nan]))[0]
+            cmap_expected(np.ma.masked_invalid([np.nan]))[0]
+            != mpl.colormaps["viridis"](np.ma.masked_invalid([np.nan]))[0]
         )
 
-        cmap.set_under("r")
+        cmap_expected.set_under("r")
         # check we actually changed the set_under color
-        assert cmap(-np.inf) != mpl.cm.viridis(-np.inf)
+        assert cmap_expected(-np.inf) != mpl.colormaps["viridis"](-np.inf)
 
-        cmap.set_over("g")
+        cmap_expected.set_over("g")
         # check we actually changed the set_over color
-        assert cmap(np.inf) != mpl.cm.viridis(-np.inf)
+        assert cmap_expected(np.inf) != mpl.colormaps["viridis"](-np.inf)
 
         # copy to ensure cmap is not changed by contourf()
-        pl = a.plot.contourf(cmap=copy(cmap))
+        pl = a.plot.contourf(cmap=copy(cmap_expected))
+        cmap = pl.cmap
+        assert cmap is not None
 
         # check the set_bad color has been kept
         assert_array_equal(
-            pl.cmap(np.ma.masked_invalid([np.nan]))[0],
             cmap(np.ma.masked_invalid([np.nan]))[0],
+            cmap_expected(np.ma.masked_invalid([np.nan]))[0],
         )
 
         # check the set_under color has been kept
-        assert pl.cmap(-np.inf) == cmap(-np.inf)
+        assert cmap(-np.inf) == cmap_expected(-np.inf)
 
         # check the set_over color has been kept
-        assert pl.cmap(np.inf) == cmap(np.inf)
+        assert cmap(np.inf) == cmap_expected(np.inf)
 
     def test3d(self) -> None:
         self.darray.plot()
@@ -723,6 +729,13 @@ class TestPlot(PlotTestCase):
         expected = "dim_0_bins_center [m]"
         assert actual == expected
 
+    def test_multiplot_over_length_one_dim(self) -> None:
+        a = easy_array((3, 1, 1, 1))
+        d = DataArray(a, dims=("x", "col", "row", "hue"))
+        d.plot(col="col")
+        d.plot(row="row")
+        d.plot(hue="hue")
+
 
 class TestPlot1D(PlotTestCase):
     @pytest.fixture(autouse=True)
@@ -838,19 +851,25 @@ class TestPlotStep(PlotTestCase):
         """Test step plot with intervals."""
         bins = [-1, 0, 1, 2]
         self.darray.groupby_bins("dim_0", bins).mean(...).plot.step()
-        assert len(plt.gca().lines[0].get_xdata()) == ((len(bins) - 1) * 2)
+        line = plt.gca().lines[0]
+        assert isinstance(line, mpl.lines.Line2D)
+        assert len(np.asarray(line.get_xdata())) == ((len(bins) - 1) * 2)
 
     def test_coord_with_interval_step_x(self) -> None:
         """Test step plot with intervals explicitly on x axis."""
         bins = [-1, 0, 1, 2]
         self.darray.groupby_bins("dim_0", bins).mean(...).plot.step(x="dim_0_bins")
-        assert len(plt.gca().lines[0].get_xdata()) == ((len(bins) - 1) * 2)
+        line = plt.gca().lines[0]
+        assert isinstance(line, mpl.lines.Line2D)
+        assert len(np.asarray(line.get_xdata())) == ((len(bins) - 1) * 2)
 
     def test_coord_with_interval_step_y(self) -> None:
         """Test step plot with intervals explicitly on y axis."""
         bins = [-1, 0, 1, 2]
         self.darray.groupby_bins("dim_0", bins).mean(...).plot.step(y="dim_0_bins")
-        assert len(plt.gca().lines[0].get_xdata()) == ((len(bins) - 1) * 2)
+        line = plt.gca().lines[0]
+        assert isinstance(line, mpl.lines.Line2D)
+        assert len(np.asarray(line.get_xdata())) == ((len(bins) - 1) * 2)
 
     def test_coord_with_interval_step_x_and_y_raises_valueeerror(self) -> None:
         """Test that step plot with intervals both on x and y axes raises an error."""
@@ -890,8 +909,11 @@ class TestPlotHistogram(PlotTestCase):
         self.pass_in_axis(self.darray.plot.hist)
 
     def test_primitive_returned(self) -> None:
-        h = self.darray.plot.hist()
-        assert isinstance(h[-1][0], mpl.patches.Rectangle)
+        n, bins, patches = self.darray.plot.hist()
+        assert isinstance(n, np.ndarray)
+        assert isinstance(bins, np.ndarray)
+        assert isinstance(patches, mpl.container.BarContainer)
+        assert isinstance(patches[0], mpl.patches.Rectangle)
 
     @pytest.mark.slow
     def test_plot_nans(self) -> None:
@@ -935,9 +957,9 @@ class TestDetermineCmapParams:
             assert cmap_params["cmap"] == "magma"
 
     def test_cmap_sequential_explicit_option(self) -> None:
-        with xr.set_options(cmap_sequential=mpl.cm.magma):
+        with xr.set_options(cmap_sequential=mpl.colormaps["magma"]):
             cmap_params = _determine_cmap_params(self.data)
-            assert cmap_params["cmap"] == mpl.cm.magma
+            assert cmap_params["cmap"] == mpl.colormaps["magma"]
 
     def test_cmap_divergent_option(self) -> None:
         with xr.set_options(cmap_divergent="magma"):
@@ -1177,7 +1199,7 @@ class TestDiscreteColorMap:
     def test_discrete_colormap_int_levels(self) -> None:
         for extend, levels, vmin, vmax, cmap in [
             ("neither", 7, None, None, None),
-            ("neither", 7, None, 20, mpl.cm.RdBu),
+            ("neither", 7, None, 20, mpl.colormaps["RdBu"]),
             ("both", 7, 4, 8, None),
             ("min", 10, 4, 15, None),
         ]:
@@ -1727,8 +1749,8 @@ class TestContour(Common2dMixin, PlotTestCase):
     # matplotlib cmap.colors gives an rgbA ndarray
     # when seaborn is used, instead we get an rgb tuple
     @staticmethod
-    def _color_as_tuple(c):
-        return tuple(c[:3])
+    def _color_as_tuple(c: Any) -> tuple[Any, Any, Any]:
+        return c[0], c[1], c[2]
 
     def test_colors(self) -> None:
         # with single color, we don't want rgb array
@@ -1750,10 +1772,16 @@ class TestContour(Common2dMixin, PlotTestCase):
         # https://github.com/pydata/xarray/issues/3284
         levels = np.array([-0.5, 0.0, 0.5, 1.0])
         artist = self.darray.plot.contour(levels=levels, colors=["k", "r", "w", "b"])
-        assert self._color_as_tuple(artist.cmap.colors[1]) == (1.0, 0.0, 0.0)
-        assert self._color_as_tuple(artist.cmap.colors[2]) == (1.0, 1.0, 1.0)
+        cmap = artist.cmap
+        assert isinstance(cmap, mpl.colors.ListedColormap)
+        colors = cmap.colors
+        assert isinstance(colors, list)
+
+        assert self._color_as_tuple(colors[1]) == (1.0, 0.0, 0.0)
+        assert self._color_as_tuple(colors[2]) == (1.0, 1.0, 1.0)
         # the last color is now under "over"
-        assert self._color_as_tuple(artist.cmap._rgba_over) == (0.0, 0.0, 1.0)
+        assert hasattr(cmap, "_rgba_over")
+        assert self._color_as_tuple(cmap._rgba_over) == (0.0, 0.0, 1.0)
 
     def test_cmap_and_color_both(self) -> None:
         with pytest.raises(ValueError):
@@ -1805,7 +1833,9 @@ class TestPcolormesh(Common2dMixin, PlotTestCase):
         artist = self.plotmethod(x="x2d", y="y2d", ax=ax)
         assert isinstance(artist, mpl.collections.QuadMesh)
         # Let cartopy handle the axis limits and artist size
-        assert artist.get_array().size <= self.darray.size
+        arr = artist.get_array()
+        assert arr is not None
+        assert arr.size <= self.darray.size
 
 
 class TestPcolormeshLogscale(PlotTestCase):
@@ -1956,6 +1986,7 @@ class TestImshow(Common2dMixin, PlotTestCase):
     ) -> None:
         da = DataArray(easy_array((5, 5, 3), start=-0.6, stop=1.4))
         arr = da.plot.imshow(vmin=vmin, vmax=vmax, robust=robust).get_array()
+        assert arr is not None
         assert 0 <= arr.min() <= arr.max() <= 1
 
     def test_normalize_rgb_one_arg_error(self) -> None:
@@ -1972,7 +2003,10 @@ class TestImshow(Common2dMixin, PlotTestCase):
         da = DataArray(np.arange(75, dtype="uint8").reshape((5, 5, 3)))
         _, ax = plt.subplots()
         out = da.plot.imshow(ax=ax).get_array()
-        assert out.dtype == np.uint8
+        assert out is not None
+        dtype = out.dtype
+        assert dtype is not None
+        assert dtype == np.uint8
         assert (out[..., :3] == da.values).all()  # Compare without added alpha
 
     @pytest.mark.filterwarnings("ignore:Several dimensions of this array")
@@ -2007,6 +2041,7 @@ class TestSurface(Common2dMixin, PlotTestCase):
         self.plotmethod(x="x2d", y="y2d")
         # make sure labels came out ok
         ax = plt.gca()
+        assert isinstance(ax, mpl_toolkits.mplot3d.axes3d.Axes3D)
         assert "x2d" == ax.get_xlabel()
         assert "y2d" == ax.get_ylabel()
         assert f"{self.darray.long_name} [{self.darray.units}]" == ax.get_zlabel()
@@ -2129,6 +2164,7 @@ class TestFacetGrid(PlotTestCase):
         self.g.map_dataarray(xplt.imshow, "x", "y")
 
         for image in plt.gcf().findobj(mpl.image.AxesImage):
+            assert isinstance(image, mpl.image.AxesImage)
             clim = np.array(image.get_clim())
             assert np.allclose(expected, clim)
 
@@ -2139,7 +2175,9 @@ class TestFacetGrid(PlotTestCase):
         fg: xplt.FacetGrid = ds.plot.scatter(x="a", y="a", row="x", hue="a")
         cbar = fg.cbar
         assert cbar is not None
+        assert hasattr(cbar, "vmin")
         assert cbar.vmin == 0
+        assert hasattr(cbar, "vmax")
         assert cbar.vmax == 3
 
     @pytest.mark.slow
@@ -2206,6 +2244,7 @@ class TestFacetGrid(PlotTestCase):
         self.g.map_dataarray(xplt.imshow, "x", "y", vmin=vmin, vmax=vmax)
 
         for image in plt.gcf().findobj(mpl.image.AxesImage):
+            assert isinstance(image, mpl.image.AxesImage)
             clim = np.array(image.get_clim())
             assert np.allclose(expected, clim)
 
@@ -2222,6 +2261,7 @@ class TestFacetGrid(PlotTestCase):
         norm = mpl.colors.SymLogNorm(0.1)
         self.g.map_dataarray(xplt.imshow, "x", "y", norm=norm)
         for image in plt.gcf().findobj(mpl.image.AxesImage):
+            assert isinstance(image, mpl.image.AxesImage)
             assert image.norm is norm
 
     @pytest.mark.slow
@@ -2759,15 +2799,20 @@ class TestDatasetScatterPlots(PlotTestCase):
         ds2 = self.ds.copy()
         ds2["hue"] = ["a", "b", "c", "d"]
         pc = ds2.plot.scatter(x="A", y="B", markersize="hue")
+        axes = pc.axes
+        assert axes is not None
         # should make a discrete legend
-        assert pc.axes.legend_ is not None
+        assert hasattr(axes, "legend_")
+        assert axes.legend_ is not None
 
     def test_legend_labels(self) -> None:
         # regression test for #4126: incorrect legend labels
         ds2 = self.ds.copy()
         ds2["hue"] = ["a", "a", "b", "b"]
         pc = ds2.plot.scatter(x="A", y="B", markersize="hue")
-        actual = [t.get_text() for t in pc.axes.get_legend().texts]
+        axes = pc.axes
+        assert axes is not None
+        actual = [t.get_text() for t in axes.get_legend().texts]
         expected = ["hue", "a", "b"]
         assert actual == expected
 
@@ -2788,7 +2833,9 @@ class TestDatasetScatterPlots(PlotTestCase):
 
     def test_add_legend_by_default(self) -> None:
         sc = self.ds.plot.scatter(x="A", y="B", hue="hue")
-        assert len(sc.figure.axes) == 2
+        fig = sc.figure
+        assert fig is not None
+        assert len(fig.axes) == 2
 
 
 class TestDatetimePlot(PlotTestCase):
@@ -2841,6 +2888,7 @@ class TestDatetimePlot(PlotTestCase):
 
         p = da.plot.pcolormesh()
         ax = p.axes
+        assert ax is not None
 
         # Make sure only mpl converters are used, use type() so only
         # mpl.dates.AutoDateLocator passes and no other subclasses: