Calculate the background and RMS of an astronomical image

calc_bg.py

import numpy as np

def calc_bg(a):
    """ Estimate the background level and rms. """
    good = ~np.isnan(a)
    assert good.sum(), 'no good pixels!'
    # poor man's source detection...
    vmax = np.percentile(a[good], 80)
    c0 = a[good] < vmax
    temp = a[good][c0]
    bg = np.median(temp)
    # now find the rms in the background
    belowbg = temp[temp < bg]
    # remove lowest 2% to get rid of any outliers
    flo = np.percentile(belowbg, 2)
    belowbg = belowbg[belowbg > flo]
    rms = np.concatenate([belowbg, 2*bg - belowbg]).std()
    return bg, rms

Thanks, I wasn't able to find this for a long time.