def banddepth(X, thr=None, norm=False):
    """
    Compute the modified band depth (Lopez-Pintado and Romo, 2009) for a
    k-member ensemble data set.

    Implementation of the exact fast algorithm for computing the modified band
    depth as described in Sun et al (2012).

    Parameters
    ----------
    X: array_like
        Array of shape (k, m, ...) representing an ensemble of *k* members
        (i.e., samples) with shape (m, ...).
    thr: float
        Optional threshold for excluding pixels that have no samples equal or
        above the **thr** value.

    Returns
    -------
    out: array_like
        Array of shape *k* containing the (normalized) band depth values for
        each ensemble member.

    References
    ----------
    Lopez-Pintado, Sara, and Juan Romo. 2009. "On the Concept of Depth for
    Functional Data." Journal of the American Statistical Association
    104 (486): 718–34. https://doi.org/10.1198/jasa.2009.0108.

    Sun, Ying, Marc G. Genton, and Douglas W. Nychka. 2012. "Exact Fast
    Computation of Band Depth for Large Functional Datasets: How Quickly Can
    One Million Curves Be Ranked?" Stat 1 (1): 68–74.
    https://doi.org/10.1002/sta4.8.
    """

    # mask invalid pixels
    if thr is None:
        thr = np.nanmin(X)
    mask = np.logical_and(np.all(np.isfinite(X), axis=0), np.any(X >= thr, axis=0))

    n = X.shape[0]
    p = np.sum(mask)
    depth = np.zeros(n)

    # assign ranks
    b = np.random.random((n, p))
    order = np.lexsort((b, X[:, mask]), axis=0)  # random rank for ties
    rank = order.argsort(axis=0) + 1

    # compute band depth
    nabove = n - rank
    nbelow = rank - 1
    match = nabove * nbelow
    nchoose2 = comb(n, 2)
    proportion = np.sum(match, axis=1) / p
    depth = (proportion + n - 1) / nchoose2

    # normalize depth between 0 and 1
    if norm:
        depth = (depth - depth.min()) / (depth.max() - depth.min())

    return depth
Compute the modified band depth (Lopez-Pintado and Romo, 2009) for a
k-member ensemble data set.

Implementation of the exact fast algorithm for computing the modified band
depth as described in Sun et al (2012).

Parameters
----------
X: array_like
Array of shape (k, m, ...) representing an ensemble of *k* members
(i.e., samples) with shape (m, ...).
thr: float
Optional threshold for excluding pixels that have no samples equal or
above the **thr** value.

Returns
-------
out: array_like
Array of shape *k* containing the (normalized) band depth values for
each ensemble member.

References
----------
Lopez-Pintado, Sara, and Juan Romo. 2009. "On the Concept of Depth for
Functional Data." Journal of the American Statistical Association
104 (486): 718–34. https://doi.org/10.1198/jasa.2009.0108.

Sun, Ying, Marc G. Genton, and Douglas W. Nychka. 2012. "Exact Fast
Computation of Band Depth for Large Functional Datasets: How Quickly Can
One Million Curves Be Ranked?" Stat 1 (1): 68–74.
https://doi.org/10.1002/sta4.8.
pysteps/postprocessing/ensemblestats.py