Attempt to JWST Data Analysis

import ginga
import astropy.io.fits as fits

# Load the astronomy dataset
dataset = fits.open('dataset.fits')

# Create a Ginga viewer
viewer = ginga.ginga.Viewer()

# Add the dataset to the viewer
viewer.add_data(dataset)

# Drizzle the dataset
viewer.do('AstroDrizzle')

# Save the drizzled dataset
viewer.save('drizzled_dataset.fits')

import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits
from astropy.stats import gaussian_sigma_to_fwhm
from astropy.table import Table
from astropy.utils.misc import NumpyRNGContext
from astropy.visualization import imshow_norm, PercentileInterval, SqrtStretch
from photutils import BasicPSFPhotometry
from photutils import DBSCANGroup, MMMBackground
from photutils.datasets import make_noise_image
from webbpsf.gridded_library import display_psf_grid
from webbpsf.utils import to_griddedpsfmodel


def main():
    """
    Perform PSF-fitting photometry on a simulated starfield.
    """

    # Load the PSF model
    psfmodel = to_griddedpsfmodel(fits.open('https://stsci.box.com/shared/static/6h8wsr2ts0t24s79cxmnyk8bldt0vv3i.fits', ext=0))

    # Simulate a starfield
    shape = (512, 512)
    data = np.zeros(shape)
    nstars = 25
    flux_max = 200000.

    with NumpyRNGContext(12345):
        xx = np.random.uniform(low=0, high=shape[1], size=nstars)
        yy = np.random.uniform(low=0, high=shape[0], size=nstars)
        zz = np.random.uniform(low=0, high=flux_max, size=nstars)

    # Evaluate the PSF model at the star positions
    for xxi, yyi, zzi in zip(xx, yy, zz):
        x0 = np.int(np.floor(xxi - (psfmodel.data.shape[2] / psfmodel.oversampling) / 2.))
        y0 = np.int(np.floor(yyi - (psfmodel.data.shape[1] / psfmodel.oversampling) / 2.))
        data[y0:y0 + psfmodel.oversampling, x0:x0 + psfmodel.oversampling] += psfmodel.evaluate(
            x=x0, y=y0, flux=zzi
        )

    # Add noise to the image
    noise = make_noise_image(data.shape, 'gaussian', mean=0, stddev=2, random_state=123)
    data += noise

    # Perform PSF-fitting photometry
    daogroup = DBSCANGroup(2.0 * gaussian_sigma_to_fwhm(3.0))
    mmm_bkg = MMMBackground()
    fit_shape = (data.shape[0] // psfmodel.oversampling, data.shape[1] // psfmodel.oversampling)
    phot = BasicPSFPhotometry(daogroup, mmm_bkg, psfmodel, fit_shape, finder=None, aperture_radius=3.)
    tbl = phot(data, init_guesses={'x_0': xx, 'y_0': yy, 'flux_0': zz})

    # Print the results
    print(tbl)


if __name__ == '__main__':
    main()