import matplotlib.pyplot as plt
import numpy as np
from astropy.stats import SigmaClip
from astropy.visualization import simple_norm
from photutils.background import Background2D, MedianBackground
from photutils.datasets import make_100gaussians_image

data = make_100gaussians_image()
ny, nx = data.shape
y, x = np.mgrid[:ny, :nx]
gradient = x * y / 5000.0
data2 = data + gradient
sigma_clip = SigmaClip(sigma=3.0)
bkg_estimator = MedianBackground()
bkg = Background2D(data2, (15, 15), filter_size=(3, 3),
                   sigma_clip=sigma_clip, bkg_estimator=bkg_estimator)
data2_sub = data2 - bkg.background
norm = simple_norm(data2_sub, 'sqrt', percent=99.9)
fig, ax = plt.subplots()
ax.imshow(data2_sub, norm=norm, origin='lower')
ax.set_title('Background-subtracted Data')