From ed090679c15e99e2e8e3aa90f2b00022f9e4328e Mon Sep 17 00:00:00 2001 From: Camilla Pacifici Date: Tue, 17 Oct 2023 13:49:50 -0400 Subject: [PATCH] Some updates to 00 --- .../00_Optimal_extraction.ipynb | 55 +++++++++---------- 1 file changed, 25 insertions(+), 30 deletions(-) diff --git a/notebooks/NIRISS_WFSS_postpipeline/00_Optimal_extraction.ipynb b/notebooks/NIRISS_WFSS_postpipeline/00_Optimal_extraction.ipynb index 1315fb83c..54952a06c 100755 --- a/notebooks/NIRISS_WFSS_postpipeline/00_Optimal_extraction.ipynb +++ b/notebooks/NIRISS_WFSS_postpipeline/00_Optimal_extraction.ipynb @@ -53,13 +53,8 @@ "from scipy.ndimage import rotate\n", "from scipy.optimize import curve_fit\n", "\n", - "\n", - "from astropy.convolution import Gaussian2DKernel\n", "from astropy.io import fits\n", - "from astropy.stats import gaussian_fwhm_to_sigma\n", - "from astropy.table import QTable\n", "import astropy.units as u\n", - "from astropy.visualization import make_lupton_rgb, SqrtStretch, ImageNormalize, simple_norm\n", "import astropy.wcs as wcs\n", "from astropy.io import ascii\n", "\n", @@ -128,7 +123,7 @@ "DIR_DATA = './pipeline_products/'\n", "\n", "# Output directory;\n", - "DIR_OUT = './output/'\n", + "DIR_OUT = './output/'\n", "if not os.path.exists(DIR_OUT):\n", " os.mkdir(DIR_OUT)\n", "\n", @@ -137,11 +132,11 @@ "\n", "# Image array from direct image. This is for optimal extraction and masking.\n", "# This image should already be sky-subtracted; otherwise, you will encounter a wrong result with optimal extraction.\n", - "infile = '%sl3_nis_%s_i2d_skysub.fits'%(DIR_DATA,filt_det)\n", + "infile= '{}13_nis_{}_i2d_skysub.fits'.format(DIR_DATA, filt_det)\n", "hdu = fits.open(infile)\n", "\n", "# This is just for error array;\n", - "infile = '%sl3_nis_%s_i2d.fits'%(DIR_DATA,filt_det)\n", + "infile = '{}l3_nis_{}_i2d.fits'.format(DIR_DATA, filt_det)\n", "hdu_err = fits.open(infile)\n", "\n", "data = hdu[0].data\n", @@ -152,7 +147,7 @@ "\n", "# Segmentation map;\n", "# This can be prepared by running Photutils, if the pipeline does not generate one.\n", - "segfile = '%sl3_nis_%s_i2d_seg.fits'%(DIR_DATA, filt_det)\n", + "segfile = '{}l3_nis_{}_i2d_seg.fits'.format(DIR_DATA, filt_det)\n", "seghdu = fits.open(segfile)\n", "segdata = seghdu[0].data" ] @@ -165,8 +160,8 @@ "source": [ "# Load catalog from image level3;\n", "# to obtain source position in pixel coordinate.\n", - "catfile = '%sl3_nis_%s_cat.ecsv'%(DIR_DATA, filt_det)\n", - "fd = ascii.read('%s'%catfile)" + "catfile = '{}l3_nis_{}_cat.ecsv'.format(DIR_DATA, filt_det)\n", + "fd = ascii.read(catfile)" ] }, { @@ -201,15 +196,15 @@ "magzp = 25.0 # magnitude zeropoint in the catalog.\n", "\n", "# Read catalog;\n", - "fd_input = ascii.read('%ssources_extend_01.cat'%(DIR_DATA))\n", + "fd_input = ascii.read('{}sources_extend_01.cat'.format(DIR_DATA))\n", "ra_input = fd_input['x_or_RA']\n", "dec_input = fd_input['y_or_Dec']\n", "\n", "# Header;\n", - "fw = open('%sl3_nis_flux.cat'%(DIR_OUT), 'w')\n", + "fw = open('{}l3_nis_flux.cat'.format(DIR_OUT), 'w')\n", "fw.write('# id')\n", "for ff in range(len(filts)):\n", - " fw.write(' F%d E%d'%(eazy_filts[ff], eazy_filts[ff]))\n", + " fw.write(' F{} E{}'.format(eazy_filts[ff], eazy_filts[ff]))\n", "fw.write('\\n')\n", "\n", "# Contents;\n", @@ -220,9 +215,9 @@ " \n", " for ff in range(len(filts)):\n", " if ff == 0:\n", - " fw.write('%d'%(fd['id'][ii]))\n", + " fw.write('{}'.format(fd['id'][ii]))\n", "\n", - " mag = fd_input['niriss_%s_magnitude'%filts[ff]][iix]\n", + " mag = fd_input['niriss_{}_magnitude'.format(filts[ff])][iix]\n", " flux_nu = 10**((mag-magzp)/(-2.5))\n", "\n", " # Currently, the catalog does not provide proper error;\n", @@ -230,7 +225,7 @@ " \n", " flux_err_nu = flux_nu * 0.05\n", "\n", - " fw.write(' %.5e %.5e'%(flux_nu, flux_err_nu))\n", + " fw.write(f' {flux_nu:.5e} {flux_err_nu:.5e}')\n", "\n", " fw.write('\\n')\n", "fw.close()" @@ -262,7 +257,7 @@ "# Zero-indexed number for dither --- the test data here has two dither positions, so 0 or 1.\n", "ndither = 0\n", "\n", - "file_2d = '%sl3_nis_%s_%s_s%s_cal.fits'%(DIR_DATA, filt, grism, id)\n", + "file_2d = '{}l3_nis_{}_{}_{}s_cal.fits'.format(DIR_DATA, filt, grism, id)\n", "hdu_2d = fits.open(file_2d)\n", "\n", "# Align grism direction\n", @@ -306,10 +301,10 @@ "# Get light profile of the source;\n", "\n", "# Again, y is for cross-dispersion, and x is for dispersion directions.\n", - "y2d,x2d = data_2d.shape[:]\n", + "y2d, x2d = data_2d.shape[:]\n", "\n", "# Cut out segmentation map;\n", - "iix = np.where(fd['id']==int(id))[0][0]\n", + "iix = np.where(fd['id'] == int(id))[0][0]\n", "\n", "# Target position from image 3 catalog;\n", "ycen = fd['ycentroid'][iix].value\n", @@ -361,9 +356,9 @@ "outputs": [], "source": [ "for ii in range(sci_rot.shape[1]):\n", - " flux_tmp = sci_rot[:,ii]\n", - " xx_tmp = np.arange(0, len(sci_rot[:,ii]), 1)\n", - " plt.plot(xx_tmp, flux_tmp, label='x=%d'%ii)\n", + " flux_tmp = sci_rot[:, ii]\n", + " xx_tmp = np.arange(0, len(sci_rot[:, ii]), 1)\n", + " plt.plot(xx_tmp, flux_tmp, label='x={}'.format(ii))\n", "plt.legend(loc=0, fontsize=8)" ] }, @@ -376,10 +371,10 @@ "# Sum along x (disperse) direction\n", "flux_y = np.zeros(len(sci_rot[:,0]), 'float')\n", "for ii in range(sci_rot.shape[0]):\n", - " flux_y[ii] = np.sum(sci_rot[ii,:])\n", + " flux_y[ii] = np.sum(sci_rot[ii, :])\n", "\n", "# Sky subtraction, if needed.\n", - "#sky = np.mean([flux_y[0], flux_y[-1]])\n", + "# sky = np.mean([flux_y[0], flux_y[-1]])\n", "\n", "# Normalize;\n", "flux_y[:] /= flux_y.sum()\n", @@ -415,15 +410,15 @@ "wave_disp1 = np.zeros(x2d, 'float')\n", " \n", "for ii in range(x2d): # Wavelength direction.\n", - " mask_tmp = (dq_2d[:,ii] == 0) & (err_2d[:,ii]>0)\n", + " mask_tmp = (dq_2d[:, ii] == 0) & (err_2d[:, ii] > 0)\n", "\n", " # Sum within a box;\n", - " flux_disp1[ii] = np.sum(data_2d[:,ii][mask_tmp]) \n", - " err_disp1[ii] = np.sqrt(np.sum(err_2d[:,ii][mask_tmp]**2)) \n", - " wave_disp1[ii] = wave_2d[0,ii]\n", + " flux_disp1[ii] = np.sum(data_2d[:, ii][mask_tmp]) \n", + " err_disp1[ii] = np.sqrt(np.sum(err_2d[:, ii][mask_tmp]**2)) \n", + " wave_disp1[ii] = wave_2d[0, ii]\n", "\n", "plt.errorbar(wave_disp1, flux_disp1, yerr=err_disp1)\n", - "plt.xlim(1.7,2.3)" + "plt.xlim(1.7, 2.3)" ] }, {