code: def block_stats(x,y,z,ds,stat='median',bins=None):
      
      import scipy.stats as stats
      extent = ds_extent(ds)
      #[[xmin, xmax], [ymin, ymax]]
      range = [[extent[0], extent[2]], [extent[1], extent[3]]]
      if bins is None:
          bins = (ds.RasterXSize, ds.RasterYSize)
      if stat == 'max':
          stat = np.max
      elif stat == 'min':
          stat = np.min
      #block_count, xedges, yedges, bin = stats.binned_statistic_2d(x,y,z,'count',bins,range)
      block_stat, xedges, yedges, bin = stats.binned_statistic_2d(x,y,z,stat,bins,range)
      #Get valid blocks
      #if (stat == 'median') or (stat == 'mean'):
      if stat in ('median', 'mean', np.max, np.min):
          idx = ~np.isnan(block_stat)
      else:
          idx = (block_stat != 0)
      idx_idx = idx.nonzero()
      #Cell centers
      res = [(xedges[1] - xedges[0]), (yedges[1] - yedges[0])]
      out_x = xedges[:-1]+res[0]/2.0
      out_y = yedges[:-1]+res[1]/2.0
      out_x = out_x[idx_idx[0]]
      out_y = out_y[idx_idx[1]]
      out_z = block_stat[idx]
      return out_x, out_y, out_z
docs: Compute points on a regular grid (matching input GDAL Dataset) from scattered point data using specified statistic
  
      Wrapper for  scipy.stats.binned_statistic_2d
      
      Note: this is very fast for mean, std, count, but bignificantly slower for median