Skip to content

Commit

Permalink
Some updates to 00
Browse files Browse the repository at this point in the history
  • Loading branch information
Camilla Pacifici committed Oct 17, 2023
1 parent 0e391fd commit ed09067
Showing 1 changed file with 25 additions and 30 deletions.
55 changes: 25 additions & 30 deletions notebooks/NIRISS_WFSS_postpipeline/00_Optimal_extraction.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand Down Expand Up @@ -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",
Expand All @@ -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",
Expand All @@ -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"
]
Expand All @@ -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)"
]
},
{
Expand Down Expand Up @@ -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",
Expand All @@ -220,17 +215,17 @@
" \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",
" # Assuming 5% error for flux.\n",
" \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()"
Expand Down Expand Up @@ -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",
Expand Down Expand Up @@ -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",
Expand Down Expand Up @@ -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)"
]
},
Expand All @@ -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",
Expand Down Expand Up @@ -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)"
]
},
{
Expand Down

0 comments on commit ed09067

Please sign in to comment.