code: def block_stats_grid(x,y,z,ds,stat='median'):
      
      mx, my, mz = block_stats(x,y,z,ds,stat)
      gt = ds.GetGeoTransform()
      pX, pY = mapToPixel(mx, my, gt)
      shape = (ds.RasterYSize, ds.RasterXSize)
      ndv = -9999.0
      a = np.full(shape, ndv)
      a[pY.astype('int'), pX.astype('int')] = mz
      return np.ma.masked_equal(a, ndv)
docs: Fill regular grid (matching input GDAL Dataset) from scattered point data using specified statistic