import ccdproc
from astropy.nddata import CCDData
import numpy as np
import matplotlib.pyplot as plt
from astropy.modeling.functional_models import Gaussian2D
from astropy.utils.misc import NumpyRNGContext
from scipy.ndimage import uniform_filter

# Create some source signal
source = Gaussian2D(60, 70, 70, 20, 25)
data = source(*np.mgrid[0:250, 0:250])

# and another one
source = Gaussian2D(70, 150, 180, 15, 15)
data += source(*np.mgrid[0:250, 0:250])

# create some random signals
with NumpyRNGContext(1234):
    noise = np.random.default_rng().exponential(40, (250, 250))
    # remove low signal
    noise[noise < 100] = 0
    data += noise

# create a CCD object based on the data
ccd = CCDData(data, unit='adu')

# Create some plots
fig, (ax1, ax2, ax3) = plt.subplots(1, 3)
ax1.set_title('Unprocessed')
ax1.imshow(ccd, origin='lower', interpolation='none', cmap=plt.cm.gray)
ax2.set_title('Mean filtered')
ax2.imshow(uniform_filter(ccd.data, 5), origin='lower', interpolation='none', cmap=plt.cm.gray)
ax3.set_title('Median filtered')
ax3.imshow(ccdproc.median_filter(ccd, 5), origin='lower', interpolation='none', cmap=plt.cm.gray)
plt.tight_layout()
plt.show()