[Numpy-discussion] contiguous regions

[John Hunter]

I frequently want to break a 1D array into regions above and below
some threshold, identifying all such subslices where the contiguous
elements are above the threshold.  I have two related implementations
below to illustrate what I am after.  The first "crossings" is rather
naive in that it doesn't handle the case where an element is equal to
the threshold (assuming zero for the threshold in the examples below).
 The second is correct (I think) but is pure python.  Has anyone got a
nifty numpy solution for this?

import numpy as np
import matplotlib.pyplot as plt
t = np.arange(0.0123, 2, 0.05)
s = np.sin(2*np.pi*t)

def crossings(x):
    """
    return a list of (above, ind0, ind1).  ind0 and ind1 are regions
    such that the slice x[ind0:ind1]>0 when above is True and
    x[ind0:ind1]<0 when above is False
    """
    N = len(x)
    crossings = x[:-1]*x[1:]<0
    ind = np.nonzero(crossings)[0]+1
    lastind = 0
    data = []
    for i in range(len(ind)):
        above = x[lastind]>0
        thisind = ind[i]
        data.append((above, lastind, thisind))
        lastind = thisind

    # put the one past the end index if not already in
    if len(data) and data[-1]!=N-1:
        data.append((not data[-1][0], thisind, N))
    return data

def contiguous_regions(mask):
    """
    return a list of (ind0, ind1) such that mask[ind0:ind1].all() is
    True and we cover all such regions
    """

    in_region = None
    boundaries = []
    for i, val in enumerate(mask):
        if in_region is None and val:
            in_region = i
        elif in_region is not None and not val:
            boundaries.append((in_region, i))
            in_region = None

    if in_region is not None:
        boundaries.append((in_region, i+1))
    return boundaries



fig = plt.figure()
ax = fig.add_subplot(111)
ax.set_title('using crossings')

ax.plot(t, s, 'o')
ax.axhline(0)


for above, ind0, ind1 in crossings(s):
    if above: color='green'
    else: color = 'red'
    tslice = t[ind0:ind1]
    ax.axvspan(tslice[0], tslice[-1], facecolor=color, alpha=0.5)

fig = plt.figure()
ax = fig.add_subplot(111)
ax.set_title('using contiguous regions')
ax.plot(t, s, 'o')
ax.axhline(0)

for ind0, ind1 in contiguous_regions(s>0):
    tslice = t[ind0:ind1]
    ax.axvspan(tslice[0], tslice[-1], facecolor='green', alpha=0.5)

for ind0, ind1 in contiguous_regions(s<0):
    tslice = t[ind0:ind1]
    ax.axvspan(tslice[0], tslice[-1], facecolor='red', alpha=0.5)


plt.show()