The best approach I can think of for this problem is to fit a spline to the array, take the derivative, and then find all local maxima. These local maxima should represent the boundaries of peaks, which I think is what you are after. My approach:
from scipy import signal
from scipy import interpolate
import numpy as np
from numpy import linspace
x = [0,0,0,0,4,5,6,6,4,0,0,0,0,0,0,2,0,0,0,6,4,5,6,0,0,0,0,0]
s = interpolate.UnivariateSpline( linspace(0,len(x)-1,len(x)), np.array(x) )
ds = s.derivative()
slope_down_begin_points = [ p for p in signal.find_peaks_cwt( vector = [ -ds(v) for v in range(len(x)) ], widths = np.array([2]) ) if x[p-1] >= 1 ]
slope_up_begin_points = [ p for p in signal.find_peaks_cwt( vector = [ ds(v) for v in range(len(x)) ], widths = np.array([2]) ) if x[p+1] >= 1 ]
slope_up_begin_points + slope_down_begin_points
>> [4, 9, 16, 19, 23]
16
is included in this approach because it is a little micro-peak of its own, if you fiddle with the find_peaks_cwt
/UnivariateSpline
parameters you should be able to filter it out..