Python - How to efficiently bin multiple dependent variables to a common independent variable?

Setup: I have ~30 parameters (dependent variables) measured simultaneously along a common time axis with ~1Hz resolution (independent variable). I need to calculate bin averages for all parameters, with a bin width of 10 seconds on the time axis. Potentially also other statistical parameters like standard deviation and median are desired. Binning should be lower limit included, upper excluded, with the last digit of the binned time axis ending to 5. I need to do this quite often (>3M entries total for one parameter set), so I look for an efficient solution.

Context: Binned data I need is an intermediate data product, meaning it's just read data --> binning --> write binned data. That makes me believe something 'high-level' like pandas might bring too much overhead. I've read some posts on the topic on Stackoverflow, e.g. here or here. Also, I know there is scipy.stats.binned_statistic (see here) but that would call binning of the time axis for every parameter, thus being inefficient imho.

What I've done so far:

# example input time axis t:
# note #1: there can be 'gaps' when data recording is interrupted
# note #2: t is guaranteed to be strictly increasing
t = np.concatenate((np.linspace(13, 24, 12), np.linspace(56, 100, 45)))

# 'binned' time axis:
tmin, tmax = np.rint(t), np.rint(t[-1])
t_binned = np.arange((tmin-tmin%10)+5, (tmax-tmax%10)+6, 10)
# array([ 15,  25,  35,  45,  55,  65,  75,  85,  95, 105])

# calculate bin starting indices and bins
binstarts = np.append(t_binned-5, t_binned[-1]+5)
# since t is strictly increasing, I use np.searchsorted which is slightly faster than np.digitized
bins = np.searchsorted(binstarts, t, side='right')

# remove elements from the binned time axis that would not map to any values:
t_binned = t_binned[np.bincount(bins-1).astype(bool)]
# array([ 15,  25,  55,  65,  75,  85,  95, 105])

now that I have the bins, I can apply it to any parameter v_i from my input data with parameters v_n:

v_i = np.random.rand(t.size)
v_i_binavg = []
for bin_no in np.unique(bins):
v_i_binavg.append(np.nanmean(v_i[bins == bin_no]))
# potentially also other calculations here, like stdev or median

This works in principle but feels a bit hacky and I think efficiency could be improved, e.g. by getting rid of the for loop in the last step. I found this to be a bottleneck in this approach (which is basically what I do now as well). Any suggestions?

p.s. for regular, 1D binning (avg or sum), I found this approach to be very efficient - however, that excludes the whole 'mapping to a time axis' problem.