def detect_outliers(input_array, thr, coord=None, k=None, verbose=False):
    """
    Detect outliers in a (multivariate and georeferenced) dataset.

    Assume a (multivariate) Gaussian distribution and detect outliers based on
    the number of standard deviations from the mean.

    If spatial information is provided through coordinates, the outlier
    detection can be localized by considering only the k-nearest neighbours
    when computing the local mean and standard deviation.

    Parameters
    ----------
    input_array: array_like
        Array of shape (n) or (n, m), where *n* is the number of samples and
        *m* the number of variables. If *m* > 1, the Mahalanobis distance
        is used.
        All values in ``input_array`` are required to have finite values.
    thr: float
        The number of standard deviations from the mean used to define an outlier.
    coord: array_like or None, optional
        Array of shape (n, d) containing the coordinates of the input data into
        a space of *d* dimensions.
        Passing ``coord`` requires that ``k`` is not None.
    k: int or None, optional
        The number of nearest neighbours used to localize the outlier
        detection. If set to None (the default), it employs all the data points (global
        detection). Setting ``k`` requires that ``coord`` is not None.
    verbose: bool, optional
        Print out information.

    Returns
    -------
    out: array_like
        A 1-D boolean array of shape (n) with True values indicating the outliers
        detected in ``input_array``.
    """

    input_array = np.copy(input_array)

    if np.any(~np.isfinite(input_array)):
        raise ValueError("input_array contains non-finite values")

    if input_array.ndim == 1:
        nsamples = input_array.size
        nvar = 1
    elif input_array.ndim == 2:
        nsamples = input_array.shape[0]
        nvar = input_array.shape[1]
    else:
        raise ValueError(
            f"input_array must have 1 (n) or 2 dimensions (n, m), "
            f"but it has {input_array.ndim}"
        )

    if nsamples < 2:
        return np.zeros(nsamples, dtype=bool)

    if coord is not None and k is not None:
        coord = np.copy(coord)
        if coord.ndim == 1:
            coord = coord[:, None]

        elif coord.ndim > 2:
            raise ValueError(
                "coord must have 2 dimensions (n, d)," f"but it has {coord.ndim}"
            )

        if coord.shape[0] != nsamples:
            raise ValueError(
                "the number of samples in input_array does not match the "
                f"number of coordinates {nsamples}!={coord.shape[0]}"
            )

        k = np.min((nsamples, k + 1))

    # global

    if k is None or coord is None:
        if nvar == 1:
            # univariate
            zdata = np.abs(input_array - np.mean(input_array)) / np.std(input_array)
            outliers = zdata > thr
        else:
            # multivariate (mahalanobis distance)
            zdata = input_array - np.mean(input_array, axis=0)
            V = np.cov(zdata.T)
            try:
                VI = np.linalg.inv(V)
                MD = np.sqrt(np.dot(np.dot(zdata, VI), zdata.T).diagonal())
            except np.linalg.LinAlgError as err:
                warnings.warn(f"{err} during outlier detection")
                MD = np.zeros(nsamples)
            outliers = MD > thr

    # local
    else:
        tree = scipy.spatial.cKDTree(coord)
        __, inds = tree.query(coord, k=k)
        outliers = np.empty(shape=0, dtype=bool)
        for i in range(inds.shape[0]):
            if nvar == 1:
                # univariate
                thisdata = input_array[i]
                neighbours = input_array[inds[i, 1:]]
                thiszdata = np.abs(thisdata - np.mean(neighbours)) / np.std(neighbours)
                outliers = np.append(outliers, thiszdata > thr)
            else:
                # multivariate (mahalanobis distance)
                thisdata = input_array[i, :]
                neighbours = input_array[inds[i, 1:], :].copy()
                thiszdata = thisdata - np.mean(neighbours, axis=0)
                neighbours = neighbours - np.mean(neighbours, axis=0)
                V = np.cov(neighbours.T)
                try:
                    VI = np.linalg.inv(V)
                    MD = np.sqrt(np.dot(np.dot(thiszdata, VI), thiszdata.T))
                except np.linalg.LinAlgError as err:
                    warnings.warn(f"{err} during outlier detection")
                    MD = 0
                outliers = np.append(outliers, MD > thr)

    if verbose:
        print(f"--- {np.sum(outliers)} outliers detected ---")

    return outliers
Detect outliers in a (multivariate and georeferenced) dataset.

Assume a (multivariate) Gaussian distribution and detect outliers based on
the number of standard deviations from the mean.

If spatial information is provided through coordinates, the outlier
detection can be localized by considering only the k-nearest neighbours
when computing the local mean and standard deviation.

Parameters
----------
input_array: array_like
Array of shape (n) or (n, m), where *n* is the number of samples and
*m* the number of variables. If *m* > 1, the Mahalanobis distance
is used.
All values in ``input_array`` are required to have finite values.
thr: float
The number of standard deviations from the mean used to define an outlier.
coord: array_like or None, optional
Array of shape (n, d) containing the coordinates of the input data into
a space of *d* dimensions.
Passing ``coord`` requires that ``k`` is not None.
k: int or None, optional
The number of nearest neighbours used to localize the outlier
detection. If set to None (the default), it employs all the data points (global
detection). Setting ``k`` requires that ``coord`` is not None.
verbose: bool, optional
Print out information.

Returns
-------
out: array_like
A 1-D boolean array of shape (n) with True values indicating the outliers
detected in ``input_array``.
pysteps/utils/cleansing.py