def polyfit(
    obj,
    dim: Hashable,
    deg: int,
    skipna: bool | None = None,
    rcond: np.floating[Any] | float | None = None,
    w: Hashable | Any = None,
    full: bool = False,
    cov: bool | Literal["unscaled"] = False,
):
    """
    Least squares polynomial fit.

    This replicates the behaviour of `numpy.polyfit` but differs by skipping
    invalid values when `skipna = True`.

    Parameters
    ----------
    obj : Dataset or DataArray
        Object to perform the polyfit on
    dim : hashable
        Coordinate along which to fit the polynomials.
    deg : int
        Degree of the fitting polynomial.
    skipna : bool or None, optional
        If True, removes all invalid values before fitting each 1D slices of the array.
        Default is True if data is stored in a dask.array or if there is any
        invalid values, False otherwise.
    rcond : float or None, optional
        Relative condition number to the fit.
    w : hashable or Any, optional
        Weights to apply to the y-coordinate of the sample points.
        Can be an array-like object or the name of a coordinate in the dataset.
    full : bool, default: False
        Whether to return the residuals, matrix rank and singular values in addition
        to the coefficients.
    cov : bool or "unscaled", default: False
        Whether to return to the covariance matrix in addition to the coefficients.
        The matrix is not scaled if `cov='unscaled'`.

    Returns
    -------
    Dataset
        A single dataset which contains (for each "var" in the input dataset):

        [var]_polyfit_coefficients
            The coefficients of the best fit for each variable in this dataset.
        [var]_polyfit_residuals
            The residuals of the least-square computation for each variable (only included if `full=True`)
            When the matrix rank is deficient, np.nan is returned.
        [dim]_matrix_rank
            The effective rank of the scaled Vandermonde coefficient matrix (only included if `full=True`)
            The rank is computed ignoring the NaN values that might be skipped.
        [dim]_singular_values
            The singular values of the scaled Vandermonde coefficient matrix (only included if `full=True`)
        [var]_polyfit_covariance
            The covariance matrix of the polynomial coefficient estimates (only included if `full=False` and `cov=True`)

    Warns
    -----
    RankWarning
        The rank of the coefficient matrix in the least-squares fit is deficient.
        The warning is not raised with in-memory (not dask) data and `full=True`.

    See Also
    --------
    numpy.polyfit
    numpy.polyval
    xarray.polyval
    """
    variables: dict[Hashable, Variable] = {}
    skipna_da = skipna

    x = np.asarray(_ensure_numeric(obj.coords[dim]).astype(np.float64))

    xname = f"{obj[dim].name}_"
    order = int(deg) + 1
    degree_coord_values = np.arange(order)[::-1]
    lhs = np.vander(x, order)

    if rcond is None:
        rcond = x.shape[0] * np.finfo(x.dtype).eps

    # Weights:
    if w is not None:
        if isinstance(w, Hashable):
            w = obj.coords[w]
        w = np.asarray(w)
        if w.ndim != 1:
            raise TypeError("Expected a 1-d array for weights.")
        if w.shape[0] != lhs.shape[0]:
            raise TypeError(f"Expected w and {dim} to have the same length")
        lhs *= w[:, np.newaxis]

    # Scaling
    scale = np.sqrt((lhs * lhs).sum(axis=0))
    lhs /= scale

    from xarray.core import utils

    degree_dim = utils.get_temp_dimname(obj.dims, "degree")

    rank = np.linalg.matrix_rank(lhs)

    if full:
        rank = Variable(dims=(), data=rank)
        variables[xname + "matrix_rank"] = rank
        _sing = np.linalg.svd(lhs, compute_uv=False)
        variables[xname + "singular_values"] = Variable(
            dims=(degree_dim,),
            data=np.concatenate([np.full((order - rank.data,), np.nan), _sing]),
        )

    # If we have a coordinate get its underlying dimension.
    (true_dim,) = obj.coords[dim].dims

    other_coords = {
        dim: obj._variables[dim]
        for dim in set(obj.dims) - {true_dim}
        if dim in obj._variables
    }
    present_dims: set[Hashable] = set()
    for name, var in obj._variables.items():
        if name in obj._coord_names or name in obj.dims:
            continue
        if true_dim not in var.dims:
            continue

        if is_duck_dask_array(var._data) and (rank != order or full or skipna is None):
            # Current algorithm with dask and skipna=False neither supports
            # deficient ranks nor does it output the "full" info (issue dask/dask#6516)
            skipna_da = True
        elif skipna is None:
            skipna_da = bool(np.any(var.isnull()))

        if var.ndim > 1:
            rhs = var.transpose(true_dim, ...)
            other_dims = rhs.dims[1:]
            scale_da = scale.reshape(-1, *((1,) * len(other_dims)))
        else:
            rhs = var
            scale_da = scale
            other_dims = ()

        present_dims.update(other_dims)
        if w is not None:
            rhs = rhs * w.reshape(-1, *((1,) * len(other_dims)))

        with warnings.catch_warnings():
            if full:  # Copy np.polyfit behavior
                warnings.simplefilter("ignore", RankWarning)
            else:  # Raise only once per variable
                warnings.simplefilter("once", RankWarning)

            coeffs, residuals = least_squares(
                lhs, rhs.data, rcond=rcond, skipna=skipna_da
            )

        from xarray.core.dataarray import _THIS_ARRAY

        if name is _THIS_ARRAY:
            # When polyfit is called on a DataArray, ensure the resulting
            # dataset is backwards compatible with previous behavior
            name = ""
        elif isinstance(name, str):
            name = f"{name}_"
        else:
            # For other non-string names
            name = ""

        variables[name + "polyfit_coefficients"] = Variable(
            data=coeffs / scale_da, dims=(degree_dim,) + other_dims
        )

        if full or (cov is True):
            variables[name + "polyfit_residuals"] = Variable(
                data=residuals if var.ndim > 1 else residuals.squeeze(),
                dims=other_dims,
            )

        fac: Variable | int
        if cov:
            Vbase = np.linalg.inv(np.dot(lhs.T, lhs))
            Vbase /= np.outer(scale, scale)
            if cov == "unscaled":
                fac = 1
            else:
                if x.shape[0] <= order:
                    raise ValueError(
                        "The number of data points must exceed order to scale the covariance matrix."
                    )
                fac = variables[name + "polyfit_residuals"] / (x.shape[0] - order)
            variables[name + "polyfit_covariance"] = (
                Variable(data=Vbase, dims=("cov_i", "cov_j")) * fac
            )

    return type(obj)(
        data_vars=variables,
        coords={
            degree_dim: degree_coord_values,
            **{
                name: coord
                for name, coord in other_coords.items()
                if name in present_dims
            },
        },
        attrs=obj.attrs.copy(),
    )

    Least squares polynomial fit.

    This replicates the behaviour of `numpy.polyfit` but differs by skipping
    invalid values when `skipna = True`.

    Parameters
    ----------
    obj : Dataset or DataArray
        Object to perform the polyfit on
    dim : hashable
        Coordinate along which to fit the polynomials.
    deg : int
        Degree of the fitting polynomial.
    skipna : bool or None, optional
        If True, removes all invalid values before fitting each 1D slices of the array.
        Default is True if data is stored in a dask.array or if there is any
        invalid values, False otherwise.
    rcond : float or None, optional
        Relative condition number to the fit.
    w : hashable or Any, optional
        Weights to apply to the y-coordinate of the sample points.
        Can be an array-like object or the name of a coordinate in the dataset.
    full : bool, default: False
        Whether to return the residuals, matrix rank and singular values in addition
        to the coefficients.
    cov : bool or "unscaled", default: False
        Whether to return to the covariance matrix in addition to the coefficients.
        The matrix is not scaled if `cov='unscaled'`.

    Returns
    -------
    Dataset
        A single dataset which contains (for each "var" in the input dataset):

        [var]_polyfit_coefficients
            The coefficients of the best fit for each variable in this dataset.
        [var]_polyfit_residuals
            The residuals of the least-square computation for each variable (only included if `full=True`)
            When the matrix rank is deficient, np.nan is returned.
        [dim]_matrix_rank
            The effective rank of the scaled Vandermonde coefficient matrix (only included if `full=True`)
            The rank is computed ignoring the NaN values that might be skipped.
        [dim]_singular_values
            The singular values of the scaled Vandermonde coefficient matrix (only included if `full=True`)
        [var]_polyfit_covariance
            The covariance matrix of the polynomial coefficient estimates (only included if `full=False` and `cov=True`)

    Warns
    -----
    RankWarning
        The rank of the coefficient matrix in the least-squares fit is deficient.
        The warning is not raised with in-memory (not dask) data and `full=True`.

    See Also
    --------
    numpy.polyfit
    numpy.polyval
    xarray.polyval
    
xarray/computation/fit.py