You can reshape the array to a four dimensional array that reflects the desired block structure, and then sum along both axes within each block. Example:

>>> a = np.arange(24).reshape(4, 6)
>>> a
array([[ 0,  1,  2,  3,  4,  5],
       [ 6,  7,  8,  9, 10, 11],
       [12, 13, 14, 15, 16, 17],
       [18, 19, 20, 21, 22, 23]])
>>> a.reshape(2, 2, 2, 3).sum(3).sum(1)
array([[ 24,  42],
       [ 96, 114]])

If a has the shape m, n, the reshape should have the form

a.reshape(m_bins, m // m_bins, n_bins, n // n_bins)

It might be surprising, but summing up some values in the matrix isn't an easy task. This and this answers of mine give some insights, so I'm not going to repeat in much detail, but to get the best performance one must utilize the cache and SIMD/pipeline-nature floating operations on modern CPUs in the best possible way.

Numpy gets all of the above right and it is quite hard to beat. It is not impossible, but one must really be very intimate with low level optimization to be successful - and one should not expect much of improvement. Naive implementation won't be able to beat Numpy at all.

Here are my tries using cython and numba, where all of the speed-up comes from parallelization.

Let's start with the baseline-algorithm:

def bin2d(a,K):
    m_bins = a.shape[0]//K
    n_bins = a.shape[1]//K
    return a.reshape(m_bins, K, n_bins, K).sum(3).sum(1)

and measure it's performance:

import numpy as np
N,K=2000,50
a=np.arange(N*N, dtype=np.float64).reshape(N,N)
%timeit bin2d(a,K)
# 2.76 ms ± 107 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Cython:

Here is my implementation with Cython, which uses OpenMP for parallelization. To keep things simple I perform summation in-place, thus passed array will be changed (which is not the case for numpy's version):

%%cython -a  -c=/openmp --link-args=/openmp -c=/arch:AVX512
import numpy as np
import cython
from cython.parallel import prange
cdef extern from *:
    """
    //assumes direct memory and row-major-order (i.e. rows are continuous)
    double calc_bin(double *ptr, int N, int y_offset){
       double *row = ptr;
       for(int y=1;y<N;y++){
           row+=y_offset; //next row
           //there is no dependencies, so the summation can be done in parallel
           for(int x=0;x<N;x++){
               ptr[x]+=row[x];
           }
       }
       double res=0.0;
       for(int x=0;x<N;x++){
           //could be made slightly faster (summation is not in parallel), but is it needed?
           res+=ptr[x];
       }
       return res;
    }
    """
    double calc_bin(double *ptr, int N, int y_offset) nogil


@cython.boundscheck(False)
@cython.wraparound(False)
def cy_bin2d_parallel(double[:, ::1] a, int K):
    cdef int y_offset = a.shape[0]
    cdef int m_bins = a.shape[0]//K
    cdef int n_bins = a.shape[1]//K
    cdef double[:,:] res = np.empty((m_bins, n_bins), dtype=np.float64)
    cdef int i,j,k
    for k in prange(m_bins*n_bins, nogil=True):
        i = k//m_bins
        j = k%m_bins
        res[i,j] = calc_bin(&a[i*K, j*K], K, y_offset)
    return res.base

And now (after verifying that the results are the same):

a=np.arange(N*N, dtype=np.float64).reshape(N,N)
%timeit cy_bin2d_parallel(a,K)
# 1.51 ms ± 25.5 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
# without parallelization: 2.93 ms

which is about 30% faster. Using -c=/arch:AVX512 (otherwise compiled only with /Ox with MSVC) as proposed by @max9111, made the single-threaded version about 20% faster but brought only about 5% for the parallel version.

Numba:

There is the same algorithm compiled with numba (which often can beat cython due to better performance of the clang-compiler) - but the result is slightly slower than cython but beats numpy by about 20%:

import numba as nb
@nb.njit(parallel=True)
def nb_bin2d_parallel(a, K):
    m_bins = a.shape[0]//K
    n_bins = a.shape[1]//K
    res = np.zeros((m_bins, n_bins), dtype=np.float64)

    for k in nb.prange(m_bins*n_bins):
        i = k//m_bins
        j = k%m_bins
        for y in range(i*K+1, (i+1)*K):
            for x in range(j*K, (j+1)*K):
                a[i*K, x] += a[y,x]
        s=0.0
        for x in range(j*K, (j+1)*K):
            s+=a[i*K, x]
        res[i,j] = s
    return res

and now:

a=np.arange(N*N, dtype=np.float64).reshape(N,N)
%timeit nb_bin2d_parallel(a,K)
# 1.98 ms ± 162 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
# without parallelization: 5.8 ms

In a nutshell: I guess it is possible to beat the above, but there is no free lunch anymore as numpy does its job pretty well. The most potential is probably in the parallelization, but it is limited due to the memory-bandwidth-bound nature of the problem (and probably one should use a smarter strategy as mine - for 50*50 summation one possible still can see an overhead for creation/management of the threads).

There is another faster (at least for sizes of about 2000) attempt with Cython, which perform the summation not for small 50-element parts, but for the whole row thus reducing the overhead (but may have more cache misses once the row is larger than 1-2k):

%%cython -a --verbose -c=/openmp --link-args=/openmp -c=/arch:AVX512
import numpy as np
import cython
from cython.parallel import prange
cdef extern from *:
    """
    void calc_bin_row(double *ptr, int N, int y_offset, double* out){
       double *row = ptr;
       for(int y=1;y<N;y++){
           row+=y_offset; //next row
           for(int x=0;x<y_offset;x++){
               ptr[x]+=row[x];
           }
       }
       double res=0.0;
       int i=0;
       int k=0;
       for(int x=0;x<y_offset;x++){//could be made faster, but is it needed?
           res+=ptr[x];
           k++;
           if(k==N){
              k=0;
              out[i]=res;
              i++;
              res=0.0;
           }
       }
    }
    """
    void calc_bin_row(double *ptr, int N, int y_offset, double* out) nogil


@cython.boundscheck(False)
@cython.wraparound(False)
def cy_bin2d_parallel_rowise(double[:, ::1] a, int K):
    cdef int y_offset = a.shape[1]
    cdef int m_bins = a.shape[0]//K
    cdef int n_bins = a.shape[1]//K
    cdef double[:,:] res = np.empty((m_bins, n_bins), dtype=np.float64)
    cdef int i,j,k
    for k in prange(0, y_offset, K, nogil=True):
        calc_bin_row(&a[k, 0], K, y_offset, &res[k//K, 0])
    return res.base

which is already faster singlethreaded (2ms) and about 60% (1.27 ms ± 50.8 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)) faster with parallelization.

You can use a higher dimensional view of your array and take the average along the extra dimensions:

In [12]: a = np.arange(36).reshape(6, 6)

In [13]: a
Out[13]: 
array([[ 0,  1,  2,  3,  4,  5],
       [ 6,  7,  8,  9, 10, 11],
       [12, 13, 14, 15, 16, 17],
       [18, 19, 20, 21, 22, 23],
       [24, 25, 26, 27, 28, 29],
       [30, 31, 32, 33, 34, 35]])

In [14]: a_view = a.reshape(3, 2, 3, 2)

In [15]: a_view.mean(axis=3).mean(axis=1)
Out[15]: 
array([[  3.5,   5.5,   7.5],
       [ 15.5,  17.5,  19.5],
       [ 27.5,  29.5,  31.5]])

In general, if you want bins of shape (a, b) for an array of (rows, cols), your reshaping of it should be .reshape(rows // a, a, cols // b, b). Note also that the order of the .mean is important, e.g. a_view.mean(axis=1).mean(axis=3) will raise an error, because a_view.mean(axis=1) only has three dimensions, although a_view.mean(axis=1).mean(axis=2) will work fine, but it makes it harder to understand what is going on.

As is, the above code only works if you can fit an integer number of bins inside your array, i.e. if a divides rows and b divides cols. There are ways to deal with other cases, but you will have to define the behavior you want then.

The problem is that bincount isn't always returning the same shaped objects, in particular when values are missing. For example:

>>> m = np.array([[0,0,1],[1,1,0],[1,1,1]])
>>> np.apply_along_axis(np.bincount, 1, m)
array([[2, 1],
       [1, 2],
       [0, 3]])
>>> [np.bincount(m[i]) for i in range(m.shape[1])]
[array([2, 1]), array([1, 2]), array([0, 3])]

works, but:

>>> m = np.array([[0,0,0],[1,1,0],[1,1,0]])
>>> m
array([[0, 0, 0],
       [1, 1, 0],
       [1, 1, 0]])
>>> [np.bincount(m[i]) for i in range(m.shape[1])]
[array([3]), array([1, 2]), array([1, 2])]
>>> np.apply_along_axis(np.bincount, 1, m)
Traceback (most recent call last):
  File "<ipython-input-49-72e06e26a718>", line 1, in <module>
    np.apply_along_axis(np.bincount, 1, m)
  File "/usr/local/lib/python2.7/dist-packages/numpy/lib/shape_base.py", line 117, in apply_along_axis
    outarr[tuple(i.tolist())] = res
ValueError: could not broadcast input array from shape (2) into shape (1)

won't.

You could use the minlength parameter and pass it using a lambda or partial or something:

>>> np.apply_along_axis(lambda x: np.bincount(x, minlength=2), axis=1, arr=m)
array([[3, 0],
       [1, 2],
       [1, 2]])

How about this:

import numpy as np
import skimage.measure

a = np.arange(36).reshape(6, 6)

b = skimage.measure.block_reduce(a, (2,2), np.mean)

output:

a = 
[[ 0  1  2  3  4  5]
 [ 6  7  8  9 10 11]
 [12 13 14 15 16 17]
 [18 19 20 21 22 23]
 [24 25 26 27 28 29]
 [30 31 32 33 34 35]]
b = 
[[ 3.5  5.5  7.5]
 [15.5 17.5 19.5]
 [27.5 29.5 31.5]]

But instead of my 2d example, you can do that for a block size of (1, 10, 10, 10) of your data.

values, xbins, ybins = np.histogram2d(x=a[:,0], y=a[:,1]) gives the actual number of points of each bin into values. Note that many matplotlib functions index first by y, so you might need values.T depending on the use case.

Here is a visualization showing how the values can be used.

import matplotlib.pyplot as plt
import matplotlib.patheffects as path_effects
import numpy as np

x = np.linspace(-0.212, 0.233, 50)
y = x * 0.5 - 0.01

hist, xbins, ybins = np.histogram2d(x=x, y=y, bins=(np.arange(-0.25, 0.25001, 0.02), np.arange(-0.15, 0.15001, 0.02)))

fig, ax = plt.subplots(figsize=(11, 6))
for i in range(len(xbins) - 1):
    for j in range(len(ybins) - 1):
        text = ax.text((xbins[i] + xbins[i + 1]) / 2, (ybins[j] + ybins[j + 1]) / 2, f"{hist[i, j]:.0f}",
                       color='cornflowerblue', size=16, ha='center', va='center')
        text.set_path_effects([path_effects.Stroke(linewidth=3, foreground='white', alpha=0.6), path_effects.Normal()])
ax.plot(x, y, '-ro')
ax.set_xlim(xbins.min(), xbins.max())
ax.set_ylim(ybins.min(), ybins.max())
ax.set_xticks(xbins + 0.0001, minor=True)
ax.set_yticks(ybins + 0.0001, minor=True)
ax.grid(which='minor', color='dodgerblue', ls='--')
ax.set_aspect(1)
plt.show()

Just use reshape and then mean(axis=1).

As the simplest possible example:

import numpy as np

data = np.array([4,2,5,6,7,5,4,3,5,7])

print data.reshape(-1, 2).mean(axis=1)

More generally, we'd need to do something like this to drop the last bin when it's not an even multiple:

import numpy as np

width=3
data = np.array([4,2,5,6,7,5,4,3,5,7])

result = data[:(data.size // width) * width].reshape(-1, width).mean(axis=1)

print result