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.

I don't think you need cython for this, I think you're looking for numpy.bincount. Here is an example:

import numpy as np
d = np.random.random(10**5)
numbins = 10

bins = np.linspace(d.min(), d.max(), numbins+1)
# This line is not necessary, but without it the smallest bin only has 1 value.
bins = bins[1:]
digitized = bins.searchsorted(d)

bin_means = (np.bincount(digitized, weights=d, minlength=numbins) /
             np.bincount(digitized, minlength=numbins))

Update

Lets take a second to discuss why the above code is faster than the code in your question and why cython will (probably) not help much in this case. In your code when you do [digitized == i] for i in range(numbins)], you're doing numbins passes over the digitized array. If you're familiar with big O notation, that's O(n * m). On the other hand bincount does something a little different. Bincount is equivelent, more or less, to:

def bincount(digitized, Weights):
   out = zeros(digitized.max() + 1)
   for i, w = zip(digitized, Weights):
       out[i] += w
   return out

It has 1 pass (well 2 passes if you count the max) over digitized so it has complexity O(n). Also bincount is already written in C and compiled so it already has very little overhead and is very fast. Cython is most helpful when you have code which has a lot of interpreter and type-check overhead so that declaring types and compiling the code removes that overhead. Hope that's helpful.

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.

If you are already importing scipy.spatial.distance, might as well use pdist. And then you're just making a 2d histogram. Use np.histogram2d.

def binDists2d(dat, f1 = 'euclidean', f2 = 'cosine'):
    dist0 = d.pdist(dat, f1)
    dist1 = d.pdist(dat, f2)
    rng = np.array([[0, 1], [0, 1]])
    return np.histogram2d(dist0, dist1, bins = 10, range = rng)

pdist only returns the upper triangular elements. If you want to do this manually, use np.triu_indices, which you could use to generate the distances if scipy is unavailable.

def cosdist(u, v):
    return 1 - u.dot(v) / (np.linalg.norm(u) * np.linlg.norm(v))

def binDists2d(dat, f0 = lambda u, v: np.linalg.norm(u - v), f1 = cosdist):
    i, j = np.triu_indices(dat.shape[0], 1)
    dist0 = f0(dat[i], dat[j])
    dist1 = f1(dat[i], dat[j])
    rng = np.array([[0, 1], [0, 1]])
    return np.histogram2d(dist0, dist1, bins = 10, range = rng)  

EDIT: Less memory-hungry version:

def binDists2d(dat, f0, f1, n = 1, bins = 10, rng = np.array([[0, 1], [0, 1]])):
    i_, j_ = np.triu_indices(dat.shape[0], 1)
    out = np.zeros((bins, bins))
    i_, j_ = np.array_split(i_, n), np.array_split(j_, n)
    for k, (i, j) in enumerate(zip(i_, j_)):
        dist0 = f0(dat[i], dat[j])
        dist1 = f1(dat[i], dat[j])
        out += np.histogram2d(dist0, dist1, bins = bins, range = rng)
        print(str(k) + " of " + str(n) + "completed")
    return out

So you are looking for median over strided reshape:

import numpy as np
a = np.arange(24).reshape(4,6)

def median_binner(a,bin_x,bin_y):
    m,n = np.shape(a)
    strided_reshape = np.lib.stride_tricks.as_strided(a,shape=(bin_x,bin_y,m//bin_x,n//bin_y),strides = a.itemsize*np.array([(m / bin_x) * n, (n / bin_y), n, 1]))
    return np.array([np.median(col) for row in strided_reshape for col in row]).reshape(bin_x,bin_y)



print "Original Matrix:"
print a
print "\n"
bin_tester1 = median_binner(a,2,3)
print "2x3 median bin :"
print bin_tester1
print "\n"
bin_tester2 = median_binner(a,2,2)
print "2x2 median bin :"
print bin_tester2

result:

Original Matrix:
[[ 0  1  2  3  4  5]
 [ 6  7  8  9 10 11]
 [12 13 14 15 16 17]
 [18 19 20 21 22 23]]


2x3 median bin :
[[  3.5   5.5   7.5]
 [ 15.5  17.5  19.5]]


2x2 median bin :
[[  4.   7.]
 [ 16.  19.]]

Read this in order to completely understand the following line in the code:

strided_reshape = np.lib.stride_tricks.as_strided(a,shape=(bin_x,bin_y,m//bin_x,n//bin_y),strides = a.itemsize*np.array([(m / bin_x) * n, (n / bin_y), n, 1])) .

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()

This is what I came up with, I hope it's useful. It's not necessarily cleaner or better than using a loop, but maybe it'll get you started toward something better.

import numpy as np
bins_x, bins_y = 1., 1.
x = np.array([1,1,2,2,3,3,3])
y = np.array([1,1,2,2,3,3,3])
w = np.array([1,2,3,4,5,6,7], 'float')

# You can get a bin number for each point like this
x = (x // bins_x).astype('int')
y = (y // bins_y).astype('int')
shape = [x.max()+1, y.max()+1]
bin = np.ravel_multi_index([x, y], shape)

# You could get the mean by doing something like:
mean = np.bincount(bin, w) / np.bincount(bin)

# Median is a bit harder
order = bin.argsort()
bin = bin[order]
w = w[order]
edges = (bin[1:] != bin[:-1]).nonzero()[0] + 1
med_index = (np.r_[0, edges] + np.r_[edges, len(w)]) // 2
median = w[med_index]

# But that's not quite right, so maybe
median2 = [np.median(i) for i in np.split(w, edges)]

Also take a look at numpy.histogram2d

You may use scipy.stats.binned_statistic to get the mean of the data in each bin. The bins would best be created via numpy.logspace. You may then plot those means e.g. as horiziontal lines spanning the bin width or as scatter at the mean position.

import numpy as np; np.random.seed(42)
from scipy.stats import binned_statistic
import matplotlib.pyplot as plt

x = np.logspace(0,5,300)
y = np.logspace(0,5,300)+np.random.rand(300)*1.e3


fig, ax = plt.subplots()
ax.scatter(x,y, s=9)

s, edges, _ = binned_statistic(x,y, statistic='mean', bins=np.logspace(0,5,6))

ys = np.repeat(s,2)
xs = np.repeat(edges,2)[1:-1]
ax.hlines(s,edges[:-1],edges[1:], color="crimson", )

for e in edges:
    ax.axvline(e, color="grey", linestyle="--")

ax.scatter(edges[:-1]+np.diff(edges)/2, s, c="limegreen", zorder=3)

ax.set_xscale("log")
ax.set_yscale("log")
plt.show()

You obviously can't do this if you want the bins to be separated at the same $x_i$ and $y_i$ values. This is easy to see if you want to have, say, 4 bins and have 10 data points at $(0,0)$ and 10 points at $(1,1)$.

But you can use a $kd$ tree data structure in which you recursively subdivide each bin so that it contains the same number of points.