From 8dde431fd8ff165425e52431298b29cd0172603d Mon Sep 17 00:00:00 2001 From: Michael Dulude Date: Tue, 9 Jan 2024 14:51:27 -0500 Subject: [PATCH] Add WFC3 notebook 'Correcting_for_Scattered_Light_in_IR_Exposures_by_Manually_Subtracting_Bad_Reads.ipynb' (#108) * updated _toc.yml and _config.yml files * Correcting_for_Scattered_Light_in_IR_Exposures_by_Manually_Subtracting_Bad_Reads.ipynb: cleared notebook outputs. * ima_visualization_and_differencing.py: added trailing line. * Minor tweaks to index.md * Update index.md * Updated README.md * Update ci_html_build.yml * Changes made to requirements file to succesfully execute notebooks. (#116) * PEP 8 fixes and updating requirements.txt to the version of ginga that works for this notebook * adding pre-requirements.sh * last few PEP 8 corrections * notebook PEP8 corrections and adding crds to requirements.txt * Link to WFC3 notebooks in the hst-notebooks instead * Refactor file operations to use context managers --------- Co-authored-by: annierose3 Co-authored-by: annierose3 <112646956+annierose3@users.noreply.github.com> Co-authored-by: Hatice Karatay --- _config.yml | 1 - _toc.yml | 2 +- ...es_by_Manually_Subtracting_Bad_Reads.ipynb | 492 ++++++++---------- .../ima_visualization_and_differencing.py | 206 ++++---- .../pre-requirements.sh | 1 + .../requirements.txt | 3 +- 6 files changed, 341 insertions(+), 364 deletions(-) create mode 100644 notebooks/WFC3/ir_scattered_light_manual_corrections/pre-requirements.sh diff --git a/_config.yml b/_config.yml index bed71fe8f..834af357d 100644 --- a/_config.yml +++ b/_config.yml @@ -53,5 +53,4 @@ exclude_patterns: [notebooks/DrizzlePac/align_mosaics/align_mosaics.ipynb, notebooks/WFC3/filter_transformations/filter_transformations.ipynb, notebooks/WFC3/flux_conversion_tool/flux_conversion_tool.ipynb, notebooks/WFC3/ir_scattered_light_calwf3_corrections/Correcting_for_Scattered_Light_in_IR_Exposures_Using_calwf3_to_Mask_Bad_Reads.ipynb, - notebooks/WFC3/ir_scattered_light_manual_corrections/Correcting_for_Scattered_Light_in_IR_Exposures_by_Manually_Subtracting_Bad_Reads.ipynb, notebooks/WFC3/zeropoints/zeropoints.ipynb] diff --git a/_toc.yml b/_toc.yml index 3503588df..ee7023e90 100644 --- a/_toc.yml +++ b/_toc.yml @@ -70,7 +70,7 @@ parts: - file: notebooks/WFC3/image_displayer_analyzer/wfc3_image_displayer_analyzer.ipynb - file: notebooks/WFC3/ir_ima_visualization/IR_IMA_Visualization_with_an_Example_of_Time_Variable_Background.ipynb # - file: notebooks/WFC3/ir_scattered_light_calwf3_corrections/Correcting_for_Scattered_Light_in_IR_Exposures_Using_calwf3_to_Mask_Bad_Reads.ipynb -# - file: notebooks/WFC3/ir_scattered_light_manual_corrections/Correcting_for_Scattered_Light_in_IR_Exposures_by_Manually_Subtracting_Bad_Reads.ipynb + - file: notebooks/WFC3/ir_scattered_light_manual_corrections/Correcting_for_Scattered_Light_in_IR_Exposures_by_Manually_Subtracting_Bad_Reads.ipynb - file: notebooks/WFC3/persistence/wfc3_ir_persistence.ipynb - file: notebooks/WFC3/tvb_flattenramp/TVB_flattenramp_notebook.ipynb - file: notebooks/WFC3/photometry_examples/phot_examples.ipynb diff --git a/notebooks/WFC3/ir_scattered_light_manual_corrections/Correcting_for_Scattered_Light_in_IR_Exposures_by_Manually_Subtracting_Bad_Reads.ipynb b/notebooks/WFC3/ir_scattered_light_manual_corrections/Correcting_for_Scattered_Light_in_IR_Exposures_by_Manually_Subtracting_Bad_Reads.ipynb index 7b05eae8a..5c7159d55 100644 --- a/notebooks/WFC3/ir_scattered_light_manual_corrections/Correcting_for_Scattered_Light_in_IR_Exposures_by_Manually_Subtracting_Bad_Reads.ipynb +++ b/notebooks/WFC3/ir_scattered_light_manual_corrections/Correcting_for_Scattered_Light_in_IR_Exposures_by_Manually_Subtracting_Bad_Reads.ipynb @@ -19,7 +19,7 @@ "- Compare the original FLT to the reprocessed FLT image.\n", "- Compare the original DRZ to the reprocessed DRZ image.\n", "\n", - "A second method for correcting TVB (masking bad reads in the RAW image and re-running `calwf3`) can be found in the notebook [Correcting for Scattered Light in WFC3/IR Exposures: Using `calwf3` to Mask Bad Reads (O'Connor 2023)](https://github.com/spacetelescope/WFC3Library/tree/master/notebooks/ir_scattered_light_calwf3_corrections). This method performs the 'ramp \n", + "A second method for correcting TVB (masking bad reads in the RAW image and re-running `calwf3`) can be found in the notebook [Correcting for Scattered Light in WFC3/IR Exposures: Using `calwf3` to Mask Bad Reads (O'Connor 2023)](https://github.com/spacetelescope/hst_notebooks/tree/main/notebooks/WFC3/ir_scattered_light_calwf3_corrections). This method performs the 'ramp \n", "fitting' step in calwf3 and therefore removes cosmic rays from the final image. It only works, \n", "however, when TVB occurs at the beginning or end of the IR exposure and may leave sky \n", "residuals in regions flagged as IR 'blobs'.\n", @@ -68,7 +68,7 @@ "This new FLT will have a reduced total exposure time (and signal-to-noise) given the rejection \n", "of some number of reads but will now have a flat background and an improved ramp fit. \n", "\n", - "Please see the notebook [WFC3/IR IMA Visualization with An Example of Time Variable Background (O'Connor 2023)](https://github.com/spacetelescope/WFC3Library/tree/master/notebooks/ir_ima_visualization) for a walkthrough of how to identify a TVB due to scattered light. " + "Please see the notebook [WFC3/IR IMA Visualization with An Example of Time Variable Background (O'Connor 2023)](https://github.com/spacetelescope/hst_notebooks/tree/main/notebooks/WFC3/ir_ima_visualization) for a walkthrough of how to identify a TVB due to scattered light. " ] }, { @@ -82,7 +82,6 @@ "We import:\n", "\n", "- *os* for setting environment variables\n", - "- *glob* for finding lists of files\n", "- *shutil* for managing directories\n", "- *numpy* for handling array functions\n", "- *matplotlib.pyplot* for plotting data\n", @@ -94,8 +93,7 @@ "- *drizzlepac* for combining images with AstroDrizzle\n", "\n", "We import the following modules:\n", - "- *ima_visualization_and_differencing* to take the difference between reads, plot the ramp, and visualize the difference in images\n", - "\n" + "- *ima_visualization_and_differencing* to take the difference between reads, plot the ramp, and visualize the difference in images" ] }, { @@ -105,7 +103,6 @@ "outputs": [], "source": [ "import os \n", - "import glob\n", "import shutil \n", "import numpy as np\n", "from matplotlib import pyplot as plt\n", @@ -118,7 +115,7 @@ "\n", "import ima_visualization_and_differencing as diff\n", "\n", - "%matplotlib inline\n" + "%matplotlib inline" ] }, { @@ -147,8 +144,11 @@ "data_list = Observations.query_criteria(obs_id=OBS_ID)\n", "\n", "file_id = \"icqtbbbxq\"\n", - "Observations.download_products(data_list['obsid'], project = 'CALWF3', obs_id = file_id, \n", - " mrp_only = False, productSubGroupDescription = ['RAW','IMA','FLT'])" + "Observations.download_products(data_list['obsid'], \n", + " project='CALWF3', \n", + " obs_id=file_id, \n", + " mrp_only=False, \n", + " productSubGroupDescription=['RAW', 'IMA', 'FLT'])" ] }, { @@ -167,8 +167,8 @@ "if not os.path.exists('orig/'):\n", " os.mkdir('orig/')\n", "\n", - "shutil.copy(f'mastDownload/HST/{file_id}/{file_id}_ima.fits',f'orig/{file_id}_ima.fits')\n", - "shutil.copy(f'mastDownload/HST/{file_id}/{file_id}_flt.fits',f'orig/{file_id}_flt.fits') \n", + "shutil.copy(f'mastDownload/HST/{file_id}/{file_id}_ima.fits', f'orig/{file_id}_ima.fits')\n", + "shutil.copy(f'mastDownload/HST/{file_id}/{file_id}_flt.fits', f'orig/{file_id}_flt.fits') \n", "\n", "raw_file = f'mastDownload/HST/{file_id}/{file_id}_raw.fits'\n", "\n", @@ -186,7 +186,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "In this section, we show how to identify the reads impacted by the scattered light by examining the difference in count rate between reads. This section was taken from the [WFC3/IR IMA Visualization with An Example of Time Variable Background (O'Connor 2023)](https://github.com/spacetelescope/WFC3Library/tree/master/notebooks/ir_ima_visualization) notebook, which includes a more comprehensive walkthrough of identifying time variable background. \n", + "In this section, we show how to identify the reads impacted by the scattered light by examining the difference in count rate between reads. This section was taken from the [WFC3/IR IMA Visualization with An Example of Time Variable Background (O'Connor 2023)](https://github.com/spacetelescope/hst_notebooks/tree/main/notebooks/WFC3/ir_ima_visualization) notebook, which includes a more comprehensive walkthrough of identifying time variable background. \n", "\n", "Here we implement a technique to examine the count rate difference between consecutive reads. In this case, we first convert from count rate (electrons/second) back to counts (electrons) before taking the difference, as shown in equation 3 from [WFC3 ISR 2018-05](https://www.stsci.edu/files/live/sites/www/files/home/hst/instrumentation/wfc3/documentation/instrument-science-reports-isrs/_documents/2018/WFC3-2018-05.pdf).\n", "\n", @@ -201,7 +201,7 @@ "metadata": {}, "outputs": [], "source": [ - "fig = plt.figure(figsize = (20, 10))\n", + "fig = plt.figure(figsize=(20, 10))\n", "\n", "ima_filepath = f'orig/{file_id}_ima.fits'\n", "\n", @@ -209,26 +209,27 @@ "cube, integ_time = diff.read_wfc3(ima_filepath)\n", "\n", "\n", - "lhs_region = {\"x0\":50,\"x1\":250,\"y0\":100,\"y1\":900}\n", - "rhs_region = {\"x0\":700,\"x1\":900,\"y0\":100,\"y1\":900}\n", + "lhs_region = {\"x0\": 50, \"x1\": 250, \"y0\": 100, \"y1\": 900}\n", + "rhs_region = {\"x0\": 700, \"x1\": 900, \"y0\": 100, \"y1\": 900}\n", "\n", - "#Please use a limit that makes sense for your own data, when running your images through this notebook.\n", + "# Please use a limit that makes sense for your own data, when running your images through this notebook.\n", "cube[np.abs(cube) > 3] = np.nan\n", "\n", - "diff_cube = diff.compute_diff_imas(cube, integ_time, diff_method = \"instantaneous\")\n", - "median_diff_fullframe, median_diff_lhs, median_diff_rhs = diff.get_median_fullframe_lhs_rhs(diff_cube, lhs_region = lhs_region, rhs_region = rhs_region)\n", - "\n", - "plt.rc('xtick', labelsize = 20) \n", - "plt.rc('ytick', labelsize = 20) \n", + "diff_cube = diff.compute_diff_imas(cube, integ_time, diff_method=\"instantaneous\")\n", + "median_diff_fullframe, median_diff_lhs, median_diff_rhs = (\n", + " diff.get_median_fullframe_lhs_rhs(diff_cube, \n", + " lhs_region=lhs_region, \n", + " rhs_region=rhs_region))\n", + "plt.rc('xtick', labelsize=20) \n", + "plt.rc('ytick', labelsize=20) \n", "plt.rcParams.update({'font.size': 30})\n", "plt.rcParams.update({'lines.markersize': 15})\n", "\n", - "diff.plot_ramp(ima_filepath, integ_time, median_diff_fullframe, median_diff_lhs, median_diff_rhs)\n", + "diff.plot_ramp(ima_filepath, integ_time, median_diff_fullframe, \n", + " median_diff_lhs, median_diff_rhs)\n", "\n", - "plt.ylim(0.5,2.5)\n", - "_ = plt.title(filename)\n", - "\n", - "\n" + "plt.ylim(0.5, 2.5)\n", + "_ = plt.title(filename)" ] }, { @@ -253,15 +254,17 @@ "metadata": {}, "outputs": [], "source": [ - "fig = plt.figure(figsize = (20, 10))\n", + "fig = plt.figure(figsize=(20, 10))\n", "\n", "ima_filepath = f'orig/{file_id}_ima.fits'\n", "\n", + "lhs_region = {\"x0\": 50, \"x1\": 250, \"y0\": 100, \"y1\": 900}\n", + "rhs_region = {\"x0\": 700, \"x1\": 900, \"y0\": 100, \"y1\": 900}\n", "\n", - "lhs_region = {\"x0\":50,\"x1\":250,\"y0\":100,\"y1\":900}\n", - "rhs_region = {\"x0\":700,\"x1\":900,\"y0\":100,\"y1\":900}\n", - "\n", - "diff.plot_ima_difference_subplots(ima_filepath, difference_method='instantaneous', lhs_region=lhs_region, rhs_region=rhs_region)\n" + "diff.plot_ima_difference_subplots(ima_filepath, \n", + " difference_method='instantaneous',\n", + " lhs_region=lhs_region, \n", + " rhs_region=rhs_region)" ] }, { @@ -326,17 +329,14 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "scrolled": false - }, + "metadata": {}, "outputs": [], "source": [ "raw_filepath = f'{file_id}_raw.fits'\n", "\n", - "\n", - "print( f\"Querying CRDS for the reference file associated with {raw_filepath}.\")\n", + "print(f\"Querying CRDS for the reference file associated with {raw_filepath}.\")\n", "command_line_input = 'crds bestrefs --files {:} --sync-references=1 --update-bestrefs'.format(raw_file)\n", - "os.system(command_line_input);" + "os.system(command_line_input)" ] }, { @@ -350,7 +350,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "To address anomalies such as scattered light, satellite trails, and more in WFC3/IR images we can remove the individual affected reads and correct the science and error arrays accordingly. For our example image, we choose to exclude reads where the ratio of background signal is greater than 1.1 e-/s (see the notebook [WFC3/IR IMA Visualization with An Example of Time Variable Background (O'Connor 2023)](https://github.com/spacetelescope/WFC3Library/tree/master/notebooks/ir_ima_visualization) for a more complete demonstration of how we find this ratio).\n" + "To address anomalies such as scattered light, satellite trails, and more in WFC3/IR images we can remove the individual affected reads and correct the science and error arrays accordingly. For our example image, we choose to exclude reads where the ratio of background signal is greater than 1.1 e-/s (see the notebook [WFC3/IR IMA Visualization with An Example of Time Variable Background (O'Connor 2023)](https://github.com/spacetelescope/hst_notebooks/tree/main/notebooks/WFC3/ir_ima_visualization) for a more complete demonstration of how we find this ratio).\n" ] }, { @@ -393,11 +393,10 @@ "metadata": {}, "outputs": [], "source": [ + "bad_reads = [11, 12, 13, 14, 15] # note, we leave the \"zero read\" (read number 16)\n", "\n", - "bad_reads = [11,12,13,14,15] #note, we leave the \"zero read\" (read number 16)\n", - "\n", - "#### Remove any existing ima or flt products from working directory or calwf3 will not run \n", - "for ext in ['flt','ima']:\n", + "# Remove any existing ima or flt products from working directory or calwf3 will not run \n", + "for ext in ['flt', 'ima']:\n", " if os.path.exists(raw_filepath.replace('raw', ext)):\n", " os.remove(raw_filepath.replace('raw', ext))" ] @@ -405,12 +404,10 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "scrolled": false - }, + "metadata": {}, "outputs": [], "source": [ - " ##run CALWF3\n", + "# run CALWF3\n", "calwf3(raw_filepath)" ] }, @@ -439,12 +436,10 @@ "outputs": [], "source": [ "ima_filepath_new = f'{file_id}_ima.fits'\n", - "ima = fits.open(ima_filepath_new)\n", - "\n", - "cube, integ_time = diff.read_wfc3(ima_filepath_new)\n", - "\n", - "dark_file = ima[0].header['DARKFILE'].replace('iref$', os.getenv('iref'))\n", - "dark_counts, dark_time = diff.read_wfc3(dark_file)" + "with fits.open(ima_filepath_new) as ima:\n", + " cube, integ_time = diff.read_wfc3(ima_filepath_new)\n", + " dark_file = ima[0].header['DARKFILE'].replace('iref$', os.getenv('iref'))\n", + " dark_counts, dark_time = diff.read_wfc3(dark_file)" ] }, { @@ -461,8 +456,9 @@ "outputs": [], "source": [ "cube_counts = np.zeros(np.shape(cube))\n", + "\n", "for img in range(cube.shape[2]):\n", - " cube_counts[:,:,img] = cube[:,:,img]*integ_time[img]" + " cube_counts[:, :, img] = cube[:, :, img]*integ_time[img]" ] }, { @@ -478,8 +474,8 @@ "metadata": {}, "outputs": [], "source": [ - "cube_diff = np.diff(cube_counts, axis = 2)\n", - "dark_diff = np.diff(dark_counts, axis = 2)\n", + "cube_diff = np.diff(cube_counts, axis=2)\n", + "dark_diff = np.diff(dark_counts, axis=2)\n", "dt = np.diff(integ_time)" ] }, @@ -496,19 +492,17 @@ "metadata": {}, "outputs": [], "source": [ - "final_dark = dark_counts[:,:,-1] \n", + "final_dark = dark_counts[:, :, -1] \n", "\n", - "final_sci = cube_counts[:,:,-1]\n", + "final_sci = cube_counts[:, :, -1]\n", "final_time = np.ones((1024, 1024))*integ_time[-1]\n", "\n", - "\n", "if (len(bad_reads) > 0):\n", " for read in bad_reads:\n", " index = cube_counts.shape[2]-read-1\n", - " final_sci -= cube_diff[:,:,index]\n", - " final_dark -= dark_diff[:,:,index]\n", - " final_time -= dt[index]\n", - " " + " final_sci -= cube_diff[:, :, index]\n", + " final_dark -= dark_diff[:, :, index]\n", + " final_time -= dt[index]" ] }, { @@ -572,56 +566,51 @@ "metadata": {}, "outputs": [], "source": [ - "\n", - "##getting variance terms\n", - "#### Readnoise in 4 amps\n", - "RN = np.zeros((1024,1024))\n", - "RN[512: ,0:512] += ima[0].header['READNSEA']\n", - "RN[0:512,0:512] += ima[0].header['READNSEB']\n", + "# Readnoise in 4 amps\n", + "RN = np.zeros((1024, 1024))\n", + "RN[512:, 0:512] += ima[0].header['READNSEA']\n", + "RN[0:512, 0:512] += ima[0].header['READNSEB']\n", "RN[0:512, 512:] += ima[0].header['READNSEC']\n", - "RN[512: , 512:] += ima[0].header['READNSED']\n", - "\n", + "RN[512:, 512:] += ima[0].header['READNSED']\n", "\n", - "#### Gain in 4 amps\n", - "gain_2D = np.zeros((1024,1024))\n", - "gain_2D[512: ,0:512] += ima[0].header['ATODGNA']\n", - "gain_2D[0:512,0:512] += ima[0].header['ATODGNB']\n", + "# Gain in 4 amps\n", + "gain_2D = np.zeros((1024, 1024))\n", + "gain_2D[512:, 0:512] += ima[0].header['ATODGNA']\n", + "gain_2D[0:512, 0:512] += ima[0].header['ATODGNB']\n", "gain_2D[0:512, 512:] += ima[0].header['ATODGNC']\n", - "gain_2D[512: , 512:] += ima[0].header['ATODGND']\n", + "gain_2D[512:, 512:] += ima[0].header['ATODGND']\n", "\n", - "#dark image\n", - "dark_im = fits.open(dark_file)\n", + "# Dark image\n", + "with fits.open(dark_file) as dark_im:\n", "\n", - "##And grabbing the flats\n", - "pflat_file = ima[0].header['PFLTFILE'].replace('iref$', os.getenv('iref'))\n", - "pflat_im = fits.open(pflat_file)\n", - "pflat = pflat_im[1].data\n", - "\n", - "dflat_file = ima[0].header['DFLTFILE'].replace('iref$', os.getenv('iref'))\n", - "dflat_im = fits.open(dflat_file)\n", - "dflat = dflat_im[1].data\n", - "\n", - "\n", - "#Computing the final error\n", - "\n", - "#### Variance terms\n", - "\n", - "## poisson error: sqrt(signal), flux = final_sci, dark is converted to electrons\n", - "signal = (final_sci*(pflat*dflat)+final_dark*gain_2D)\n", - "## flat errors\n", - "dflat_err=dflat_im['ERR'].data\n", - "pflat_err=pflat_im['ERR'].data\n", - "#dark error\n", - "dark_err = (dark_im['ERR'].data)*gain_2D \n", - "##final\n", - "final_err = np.sqrt(RN**2 + signal + (dark_err**2) + (pflat_err*final_sci*dflat)**2\\\n", - " + (dflat_err*final_sci*pflat)**2)\n", + " # Grabbing the flats\n", + " pflat_file = ima[0].header['PFLTFILE'].replace('iref$', os.getenv('iref'))\n", + " dflat_file = ima[0].header['DFLTFILE'].replace('iref$', os.getenv('iref'))\n", " \n", - "final_err /= (pflat*dflat*final_time)\n", - "final_err[np.isnan(final_err)] = 0\n", - "\n", - "#and finally the final science image (converted back to count rate)\n", - "final_sci /= final_time " + " with fits.open(pflat_file) as pflat_im, fits.open(dflat_file) as dflat_im:\n", + " pflat = pflat_im[1].data\n", + " dflat = dflat_im[1].data\n", + " \n", + " # Computing the final error\n", + " \n", + " # Poisson error: sqrt(signal), flux = final_sci, dark is converted to electrons\n", + " signal = (final_sci*(pflat*dflat)+final_dark*gain_2D)\n", + " \n", + " # Flat errors\n", + " dflat_err = dflat_im['ERR'].data\n", + " pflat_err = pflat_im['ERR'].data\n", + " \n", + " # Dark error\n", + " dark_err = (dark_im['ERR'].data)*gain_2D \n", + " \n", + " # Final error term\n", + " final_err = np.sqrt(RN**2 + signal + (dark_err**2) + (pflat_err*final_sci*dflat)**2\n", + " + (dflat_err*final_sci*pflat)**2) \n", + " final_err /= (pflat*dflat*final_time)\n", + " final_err[np.isnan(final_err)] = 0\n", + " \n", + " # Finally the final science image is converted back to count rate\n", + " final_sci /= final_time " ] }, { @@ -645,20 +634,16 @@ "outputs": [], "source": [ "flt_filepath_new = f'{file_id}_flt.fits'\n", - "flt = fits.open(flt_filepath_new, mode='update')\n", + "with fits.open(flt_filepath_new, mode='update') as flt:\n", " \n", - "#Updating the flt data\n", - "\n", - "flt['SCI'].data = final_sci[5:-5,5:-5]\n", - "flt['ERR'].data = final_err[5:-5,5:-5]\n", - "\n", - "\n", - "# Updating the FLT header\n", - "flt[0].header['IMA2FLT'] = (1, 'FLT extracted from IMA file')\n", - "flt[0].header['EXPTIME'] = np.max(final_time)\n", - "flt[0].header['NPOP'] = (len(bad_reads), 'Number of reads popped from the sequence')\n", - "\n", - "flt.flush()" + " # Updating the flt data\n", + " flt['SCI'].data = final_sci[5:-5, 5:-5]\n", + " flt['ERR'].data = final_err[5:-5, 5:-5]\n", + " \n", + " # Updating the FLT header\n", + " flt[0].header['IMA2FLT'] = (1, 'FLT extracted from IMA file')\n", + " flt[0].header['EXPTIME'] = np.max(final_time)\n", + " flt[0].header['NPOP'] = (len(bad_reads), 'Number of reads popped from the sequence')" ] }, { @@ -683,7 +668,7 @@ }, "outputs": [], "source": [ - "def remove_reads(raw_filepath, bad_reads = []):\n", + "def remove_reads(raw_filepath, bad_reads=[]):\n", " '''\n", " From the final IMA read, subtract data of reads affected by anomalies.\n", " Compute a corrected science image and error image, and place them in the FLT product. \n", @@ -700,111 +685,96 @@ " flt_filepath = raw_filepath.replace('raw', 'flt')\n", " ima_filepath = raw_filepath.replace('raw', 'ima')\n", " \n", - " #### Remove existing products or calwf3 will die\n", + " # Remove existing products or calwf3 will die\n", " for filepath in [flt_filepath, ima_filepath]:\n", " if os.path.exists(filepath):\n", " os.remove(filepath)\n", " \n", - " #### Run calwf3\n", + " # Run calwf3\n", " calwf3(raw_filepath)\n", " \n", - " #take calwf3 produced IMA, get rid of bad reads \n", - " ima = fits.open(ima_filepath)\n", - "\n", - " cube, integ_time = diff.read_wfc3(ima_filepath)\n", - "\n", - " dark_file = ima[0].header['DARKFILE'].replace('iref$', os.getenv('iref'))\n", - " dark_counts, dark_time = diff.read_wfc3(dark_file)\n", - "\n", - "\n", - " cube_counts = np.zeros(np.shape(cube))\n", - " for img in range(cube.shape[2]):\n", - " cube_counts[:,:,img] = cube[:,:,img]*integ_time[img]\n", - "\n", - " cube_diff = np.diff(cube_counts, axis = 2)\n", - " dark_diff = np.diff(dark_counts, axis = 2)\n", - " dt = np.diff(integ_time)\n", - "\n", - " final_dark = dark_counts[:,:,-1] \n", - "\n", - " final_sci = cube_counts[:,:,-1]\n", - " final_time = np.ones((1024, 1024))*integ_time[-1]\n", - "\n", + " # Take calwf3 produced IMA, get rid of bad reads \n", + " with fits.open(ima_filepath) as ima:\n", "\n", - " if (len(bad_reads) > 0):\n", - " for read in bad_reads:\n", - " index = cube_counts.shape[2]-read-1\n", - " final_sci -= cube_diff[:,:,index]\n", - " final_dark -= dark_diff[:,:,index]\n", - " final_time -= dt[index]\n", - "\n", - "\n", - "\n", - " ##getting variance terms\n", - " #### Readnoise in 4 amps\n", - " RN = np.zeros((1024,1024))\n", - " RN[512: ,0:512] += ima[0].header['READNSEA']\n", - " RN[0:512,0:512] += ima[0].header['READNSEB']\n", - " RN[0:512, 512:] += ima[0].header['READNSEC']\n", - " RN[512: , 512:] += ima[0].header['READNSED']\n", - "\n", - "\n", - " #### Gain in 4 amps\n", - " gain_2D = np.zeros((1024,1024))\n", - " gain_2D[512: ,0:512] += ima[0].header['ATODGNA']\n", - " gain_2D[0:512,0:512] += ima[0].header['ATODGNB']\n", - " gain_2D[0:512, 512:] += ima[0].header['ATODGNC']\n", - " gain_2D[512: , 512:] += ima[0].header['ATODGND']\n", - "\n", - " #dark image\n", - " dark_im = fits.open(dark_file)\n", - "\n", - " ##And grabbing the flats\n", - " pflat_file = ima[0].header['PFLTFILE'].replace('iref$', os.getenv('iref'))\n", - " pflat_im = fits.open(pflat_file)\n", - " pflat = pflat_im[1].data\n", - "\n", - " dflat_file = ima[0].header['DFLTFILE'].replace('iref$', os.getenv('iref'))\n", - " dflat_im = fits.open(dflat_file)\n", - " dflat = dflat_im[1].data\n", - "\n", - "\n", - " #Computing the final error\n", - "\n", - " #### Variance terms\n", - "\n", - " ## poisson error: sqrt(signal), flux = final_sci, dark is converted to electrons\n", - " signal = (final_sci*(pflat*dflat)+final_dark*gain_2D)\n", - " ## flat errors\n", - " dflat_err=dflat_im['ERR'].data\n", - " pflat_err=pflat_im['ERR'].data\n", - " #dark error\n", - " dark_err = (dark_im['ERR'].data)*gain_2D \n", - " ##final\n", - " final_err = np.sqrt(RN**2 + signal + (dark_err**2) + (pflat_err*final_sci*dflat)**2\\\n", - " + (dflat_err*final_sci*pflat)**2)\n", - "\n", - " final_err /= (pflat*dflat*final_time)\n", - " final_err[np.isnan(final_err)] = 0\n", - "\n", - " #and finally the final science image (converted back to count rate)\n", - " final_sci /= final_time \n", - "\n", - "\n", - " #Updating the flt data\n", + " cube, integ_time = diff.read_wfc3(ima_filepath)\n", " \n", - " flt = fits.open(flt_filepath, mode='update')\n", - "\n", - " flt['SCI'].data = final_sci[5:-5,5:-5]\n", - " flt['ERR'].data = final_err[5:-5,5:-5]\n", - "\n", - "\n", - " # Updating the FLT header\n", - " flt[0].header['IMA2FLT'] = (1, 'FLT extracted from IMA file')\n", - " flt[0].header['EXPTIME'] = np.max(final_time)\n", - " flt[0].header['NPOP'] = (len(bad_reads), 'Number of reads popped from the sequence')\n", - "\n", - " flt.flush()\n" + " dark_file = ima[0].header['DARKFILE'].replace('iref$', os.getenv('iref'))\n", + " dark_counts, dark_time = diff.read_wfc3(dark_file)\n", + " \n", + " cube_counts = np.zeros(np.shape(cube))\n", + " for img in range(cube.shape[2]):\n", + " cube_counts[:, :, img] = cube[:, :, img]*integ_time[img]\n", + " \n", + " cube_diff = np.diff(cube_counts, axis=2)\n", + " dark_diff = np.diff(dark_counts, axis=2)\n", + " dt = np.diff(integ_time)\n", + " \n", + " final_dark = dark_counts[:, :, -1] \n", + " final_sci = cube_counts[:, :, -1]\n", + " final_time = np.ones((1024, 1024))*integ_time[-1]\n", + " \n", + " if (len(bad_reads) > 0):\n", + " for read in bad_reads:\n", + " index = cube_counts.shape[2]-read-1\n", + " final_sci -= cube_diff[:, :, index]\n", + " final_dark -= dark_diff[:, :, index]\n", + " final_time -= dt[index]\n", + " \n", + " # Variance terms\n", + " # Readnoise in 4 amps\n", + " RN = np.zeros((1024, 1024))\n", + " RN[512:, 0:512] += ima[0].header['READNSEA']\n", + " RN[0:512, 0:512] += ima[0].header['READNSEB']\n", + " RN[0:512, 512:] += ima[0].header['READNSEC']\n", + " RN[512:, 512:] += ima[0].header['READNSED']\n", + " \n", + " # Gain in 4 amps\n", + " gain_2D = np.zeros((1024, 1024))\n", + " gain_2D[512:, 0:512] += ima[0].header['ATODGNA']\n", + " gain_2D[0:512, 0:512] += ima[0].header['ATODGNB']\n", + " gain_2D[0:512, 512:] += ima[0].header['ATODGNC']\n", + " gain_2D[512:, 512:] += ima[0].header['ATODGND']\n", + " \n", + " # Dark image\n", + " with fits.open(dark_file) as dark_im:\n", + " \n", + " # Grabbing the flats\n", + " pflat_file = ima[0].header['PFLTFILE'].replace('iref$', os.getenv('iref'))\n", + " dflat_file = ima[0].header['DFLTFILE'].replace('iref$', os.getenv('iref'))\n", + " with fits.open(pflat_file) as pflat_im, fits.open(dflat_file) as dflat_im:\n", + " pflat = pflat_im[1].data\n", + " dflat = dflat_im[1].data\n", + " \n", + " # Computing the final error\n", + " \n", + " # Poisson error: sqrt(signal), flux = final_sci, dark is converted to electrons\n", + " signal = (final_sci*(pflat*dflat)+final_dark*gain_2D)\n", + " # Flat errors\n", + " dflat_err = dflat_im['ERR'].data\n", + " pflat_err = pflat_im['ERR'].data\n", + " # Dark error\n", + " dark_err = (dark_im['ERR'].data)*gain_2D \n", + " \n", + " # Final error term\n", + " final_err = np.sqrt(RN**2 + signal + (dark_err**2) + (pflat_err*final_sci*dflat)**2\n", + " + (dflat_err*final_sci*pflat)**2)\n", + " \n", + " final_err /= (pflat*dflat*final_time)\n", + " final_err[np.isnan(final_err)] = 0\n", + " \n", + " # Finally the final science image (converted back to count rate)\n", + " final_sci /= final_time \n", + " \n", + " # Updating the flt data\n", + " with fits.open(flt_filepath, mode='update') as flt:\n", + " \n", + " flt['SCI'].data = final_sci[5:-5, 5:-5]\n", + " flt['ERR'].data = final_err[5:-5, 5:-5]\n", + " \n", + " # Updating the FLT header\n", + " flt[0].header['IMA2FLT'] = (1, 'FLT extracted from IMA file')\n", + " flt[0].header['EXPTIME'] = np.max(final_time)\n", + " flt[0].header['NPOP'] = (len(bad_reads), 'Number of reads popped from the sequence')" ] }, { @@ -829,7 +799,7 @@ "source": [ "reprocessed_flt = f'{file_id}_flt.fits'\n", "\n", - "original_flt = f'orig/{file_id}_flt.fits'\n" + "original_flt = f'orig/{file_id}_flt.fits'" ] }, { @@ -842,36 +812,34 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "scrolled": false - }, + "metadata": {}, "outputs": [], "source": [ - "image_new = fits.open(reprocessed_flt)\n", - "image_old = fits.open(original_flt)\n", - "\n", - "fig = plt.figure(figsize = (20, 7))\n", - "fig\n", - "rows = 1\n", - "columns = 2\n", - "\n", - "#add the total exptime in the title \n", - "ax1 = fig.add_subplot(rows, columns, 1)\n", - "ax1.set_title(\"Reprocessed FLT image\", fontsize = 20)\n", - "im1 = plt.imshow(image_new[\"SCI\", 1].data, vmin = 0.8, vmax = 1.3, origin = 'lower', cmap = 'Greys_r')\n", - "ax1.tick_params(axis = 'both',labelsize = 10)\n", - "cbar1 = plt.colorbar(im1, ax = ax1)\n", - "cbar1.ax.tick_params(labelsize = 10)\n", - "\n", - "ax2 = fig.add_subplot(rows, columns, 2)\n", - "ax2.set_title(\"Original FLT image\", fontsize = 20)\n", - "im2=plt.imshow(image_old[\"SCI\", 1].data, vmin = 0.8, vmax = 1.3, origin = 'lower', cmap = 'Greys_r')\n", - "ax2.tick_params(axis = 'both', labelsize = 10)\n", - "cbar2 = plt.colorbar(im2, ax = ax2)\n", - "cbar2.ax.tick_params(labelsize = 10)\n", - "\n", - "plt.rc('xtick', labelsize = 10) \n", - "plt.rc('ytick', labelsize = 10) \n" + "with fits.open(reprocessed_flt) as image_new, fits.open(original_flt) as image_old:\n", + " fig = plt.figure(figsize=(20, 7))\n", + " fig\n", + " rows = 1\n", + " columns = 2\n", + " \n", + " # Add the total exptime in the title \n", + " ax1 = fig.add_subplot(rows, columns, 1)\n", + " ax1.set_title(\"Reprocessed FLT image\", fontsize=20)\n", + " im1 = plt.imshow(image_new[\"SCI\", 1].data, vmin=0.8, vmax=1.3, \n", + " origin='lower', cmap='Greys_r')\n", + " ax1.tick_params(axis='both', labelsize=10)\n", + " cbar1 = plt.colorbar(im1, ax=ax1)\n", + " cbar1.ax.tick_params(labelsize=10)\n", + " \n", + " ax2 = fig.add_subplot(rows, columns, 2)\n", + " ax2.set_title(\"Original FLT image\", fontsize=20)\n", + " im2 = plt.imshow(image_old[\"SCI\", 1].data, vmin=0.8, vmax=1.3, \n", + " origin='lower', cmap='Greys_r')\n", + " ax2.tick_params(axis='both', labelsize=10)\n", + " cbar2 = plt.colorbar(im2, ax=ax2)\n", + " cbar2.ax.tick_params(labelsize=10)\n", + " \n", + " plt.rc('xtick', labelsize=10) \n", + " plt.rc('ytick', labelsize=10) " ] }, { @@ -909,8 +877,8 @@ "source": [ "data_list = Observations.query_criteria(obs_id=OBS_ID)\n", "\n", - "Observations.download_products(data_list['obsid'], project='CALWF3', \n", - " mrp_only=False, productSubGroupDescription=['FLT','DRZ'])" + "Observations.download_products(data_list['obsid'], project='CALWF3', \n", + " mrp_only=False, productSubGroupDescription=['FLT', 'DRZ'])" ] }, { @@ -923,7 +891,8 @@ "\n", "nominal_list = []\n", "for nominal_file_id in nominal_file_ids:\n", - " shutil.copy(f'mastDownload/HST/{nominal_file_id}/{nominal_file_id}_flt.fits', f'{nominal_file_id}_flt.fits')\n", + " shutil.copy(f'mastDownload/HST/{nominal_file_id}/{nominal_file_id}_flt.fits', \n", + " f'{nominal_file_id}_flt.fits')\n", " nominal_list.append(f'{nominal_file_id}_flt.fits')\n", "print(nominal_list)" ] @@ -938,9 +907,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "scrolled": false - }, + "metadata": {}, "outputs": [], "source": [ "updatewcs.updatewcs(reprocessed_flt, use_db=True)" @@ -961,7 +928,9 @@ }, "outputs": [], "source": [ - "astrodrizzle.AstroDrizzle('*flt.fits', output='f140w', mdriztab=True, preserve=False, build=False, context=False)" + "astrodrizzle.AstroDrizzle('*flt.fits', output='f140w', \n", + " mdriztab=True, preserve=False, \n", + " build=False, context=False)" ] }, { @@ -980,22 +949,21 @@ "DRZ_image = fits.getdata(\"f140w_drz.fits\")\n", "Orig_DRZ = fits.getdata('mastDownload/HST/icqtbb020/icqtbb020_drz.fits')\n", "\n", - "fig = plt.figure(figsize=(20,10))\n", + "fig = plt.figure(figsize=(20, 10))\n", "rows = 1\n", "columns = 2\n", "\n", - "\n", - "#add the total exptime in the title \n", + "# Add the total exptime in the title \n", "ax1 = fig.add_subplot(rows, columns, 1)\n", "ax1.set_title(\"Reprocessed DRZ Image\", fontsize=20)\n", - "vmin,vmax = zscale(Orig_DRZ)\n", + "vmin, vmax = zscale(Orig_DRZ)\n", "im1 = plt.imshow(DRZ_image, vmin=vmin, vmax=vmax, origin='lower', cmap='Greys_r')\n", - "_= plt.colorbar()\n", + "_ = plt.colorbar()\n", "\n", "ax2 = fig.add_subplot(rows, columns, 2)\n", "ax2.set_title(\"Original Pipeline DRZ Image\", fontsize=20)\n", "im2 = plt.imshow(Orig_DRZ, vmin=vmin, vmax=vmax, origin='lower', cmap='Greys_r')\n", - "_= plt.colorbar()" + "_ = plt.colorbar()" ] }, { @@ -1031,7 +999,7 @@ "\n", "**Author:** Anne O'Connor, Jennifer Mack, Annalisa Calamida, Harish Khandrika -- WFC3 Instrument\n", "\n", - "**Updated On:** 2023-05-16\n", + "**Updated On:** 2023-11-13\n", "\n", "## Citations \n", "\n", @@ -1081,9 +1049,9 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.0" + "version": "3.11.6" } }, "nbformat": 4, - "nbformat_minor": 2 + "nbformat_minor": 4 } diff --git a/notebooks/WFC3/ir_scattered_light_manual_corrections/ima_visualization_and_differencing.py b/notebooks/WFC3/ir_scattered_light_manual_corrections/ima_visualization_and_differencing.py index f242bedba..d43e744b7 100644 --- a/notebooks/WFC3/ir_scattered_light_manual_corrections/ima_visualization_and_differencing.py +++ b/notebooks/WFC3/ir_scattered_light_manual_corrections/ima_visualization_and_differencing.py @@ -5,6 +5,7 @@ from ginga.util.zscale import zscale import matplotlib.patheffects as path_effects + def read_wfc3(filename): ''' Read a full-frame IR image and return the datacube plus integration times for each read. @@ -24,26 +25,23 @@ def read_wfc3(filename): integ_time : array-like Integration times associated with the datacube in ascending order. ''' - with fits.open(filename) as f: hdr = f[0].header NSAMP = hdr['NSAMP'] hdr1 = f[1].header - cube = np.zeros((hdr1['NAXIS1'],hdr1['NAXIS2'],NSAMP), dtype = float) - integ_time = np.zeros(shape = (NSAMP)) + cube = np.zeros((hdr1['NAXIS1'], hdr1['NAXIS2'], NSAMP), dtype=float) + integ_time = np.zeros(shape=NSAMP) for i in range(1, NSAMP+1): - cube[:,:,i-1] = f[('SCI', i)].data + cube[:, :, i-1] = f[('SCI', i)].data integ_time[i-1] = f[('TIME', i)].header['PIXVALUE'] - cube = cube[:,:,::-1] + cube = cube[:, :, ::-1] integ_time = integ_time[::-1] return cube, integ_time - def compute_diff_imas(cube, integ_time, diff_method): - ''' Compute the difference in signal between reads of a WFC3 IR IMA file. @@ -66,10 +64,9 @@ def compute_diff_imas(cube, integ_time, diff_method): 1024x1024x(NSAMP-1) datacube of the differebce between IR IMA reads in ascending time order, where NSAMP is the number of samples taken. ''' - if diff_method == 'instantaneous': ima_j = cube[:, :, 1:] - ima_j_1 = cube[:,:,0:-1] + ima_j_1 = cube[:, :, 0:-1] t_0 = integ_time[0] t_j = integ_time[1:] t_j_1 = integ_time[0:-1] @@ -77,7 +74,7 @@ def compute_diff_imas(cube, integ_time, diff_method): diff = ((ima_j*(t_j-t_0))-(ima_j_1*(t_j_1-t_0)))/(t_j-t_j_1) elif diff_method == 'cumulative': - diff = cube[:,:,0:-1] - cube[:,:,1:] + diff = cube[:, :, 0:-1] - cube[:, :, 1:] else: # if an incorrect method is chosen raise an error raise ValueError(f"{diff_method} is an invalid method. The allowed methods are 'instantaneous' and 'cumulative'.") @@ -86,7 +83,6 @@ def compute_diff_imas(cube, integ_time, diff_method): def get_median_fullframe_lhs_rhs(cube, lhs_region, rhs_region): - ''' Compute the median in the full-frame image, the user-defined left side region, and the user-defined right side region. @@ -113,19 +109,19 @@ def get_median_fullframe_lhs_rhs(cube, lhs_region, rhs_region): median_rhs : array of floats The median signal of the right side of each read. ''' - - - median_full_frame = np.nanmedian(cube[5:-5,5:-5,:], axis = (0,1)) + median_full_frame = np.nanmedian(cube[5:-5, 5:-5, :], + axis=(0, 1)) median_lhs = np.nanmedian(cube[lhs_region['y0']:lhs_region['y1'], - lhs_region['x0']:lhs_region['x1'],:], axis = (0,1)) + lhs_region['x0']:lhs_region['x1'], :], + axis=(0, 1)) median_rhs = np.nanmedian(cube[rhs_region['y0']:rhs_region['y1'], - rhs_region['x0']:rhs_region['x1'],:], axis = (0,1)) - + rhs_region['x0']:rhs_region['x1'], :], + axis=(0, 1)) return median_full_frame, median_lhs, median_rhs + def get_std_fullframe_lhs_rhs(cube, lhs_region, rhs_region): - ''' Compute the standard deviation of the signal in the full-frame image, the user-defined left side region, and the user-defined right side region. @@ -154,20 +150,19 @@ def get_std_fullframe_lhs_rhs(cube, lhs_region, rhs_region): standard_dev_rhs : array of floats The standard deviation of the signal of the right side of each read. ''' - - - standard_dev_fullframe = np.nanstd(cube[5:-5,5:-5,:], axis = (0,1)) + standard_dev_fullframe = np.nanstd(cube[5:-5, 5:-5, :], + axis=(0, 1)) standard_dev_lhs = np.nanstd(cube[lhs_region['y0']:lhs_region['y1'], - lhs_region['x0']:lhs_region['x1'],:], axis = (0,1)) + lhs_region['x0']:lhs_region['x1'], :], + axis=(0, 1)) standard_dev_rhs = np.nanstd(cube[rhs_region['y0']:rhs_region['y1'], - rhs_region['x0']:rhs_region['x1'],:], axis = (0,1)) - - + rhs_region['x0']:rhs_region['x1'], :], + axis=(0, 1)) + return standard_dev_fullframe, standard_dev_lhs, standard_dev_rhs - + def plot_ramp(ima, integ_time, median_diff_fullframe, median_diff_lhs, median_diff_rhs): - ''' Plots the signal accumulation ramp of an IMA image. Each point is the median signal (in e-/s) of the difference between subsequent reads. The median signal difference is plotted for the full @@ -190,19 +185,22 @@ def plot_ramp(ima, integ_time, median_diff_fullframe, median_diff_lhs, median_di median_diff_rhs: array-like The median difference in signal between the right side of each read. ''' - - plt.plot(integ_time[2:], median_diff_fullframe[1:], 's', markersize = 25, label = 'Full Frame', color = 'black') - plt.plot(integ_time[2:], median_diff_lhs[1:], '<', markersize = 20, label = 'LHS', color = 'orange') - plt.plot(integ_time[2:], median_diff_rhs[1:], '>', markersize = 20, label = 'RHS', color = 'green') + plt.plot(integ_time[2:], median_diff_fullframe[1:], 's', markersize=25, + label='Full Frame', color='black') + plt.plot(integ_time[2:], median_diff_lhs[1:], '<', markersize=20, + label='LHS', color='orange') + plt.plot(integ_time[2:], median_diff_rhs[1:], '>', markersize=20, + label='RHS', color='green') ax = plt.gca() - for spine in ['top', 'bottom', 'left', 'right']: ax.spines[spine].set_visible(False) + for spine in ['top', 'bottom', 'left', 'right']: + ax.spines[spine].set_visible(False) plt.grid() plt.xlabel('SAMPTIME [s]') - plt.ylabel('$\mu$ [e-/s]') - plt.legend(loc = 0) + plt.ylabel(r'$\mu$ [e-/s]') + plt.legend(loc=0) plt.title(ima) - - + + def panel_plot(cube, integ_time, median_diff_full_frame, median_diff_lhs, median_diff_rhs, standard_dev_fullframe, standard_dev_lhs, standard_dev_rhs, diff_method): ''' @@ -245,7 +243,7 @@ def panel_plot(cube, integ_time, median_diff_full_frame, median_diff_lhs, median fig: figure object Panel plot with subplots showing the difference between subsequent IMA reads. - Above each panel, we print the median difference $\mu$ in the count rate over the entire image. + Above each panel, we print the median difference mu in the count rate over the entire image. Below each panel, we list the IMSET difference, along with the time interval between the two IMSETs. The statistics in orange (on the left and right side of each panel) give the median rate and standard deviation of each side of the image, respectively. The value in green 'delta' is the @@ -253,55 +251,52 @@ def panel_plot(cube, integ_time, median_diff_full_frame, median_diff_lhs, median The value in white "Ratio" gives the ratio of the median difference in orange for the left versus the right side. ''' - - - xlabel_list = ["SCI[16-15]","SCI[15-14]","SCI[14-13]","SCI[13-12]","SCI[12-11]", - "SCI[11-10]","SCI[10-9]","SCI[9-8]","SCI[8-7]","SCI[[7-6]]","SCI[6-5]", - "SCI[5-4]","SCI[4-3]","SCI[3-2]","SCI[2-1]"] + xlabel_list = ["SCI[16-15]", "SCI[15-14]", "SCI[14-13]", "SCI[13-12]", "SCI[12-11]", + "SCI[11-10]", "SCI[10-9]", "SCI[9-8]", "SCI[8-7]", "SCI[[7-6]]", "SCI[6-5]", + "SCI[5-4]", "SCI[4-3]", "SCI[3-2]", "SCI[2-1]"] fig, axarr = plt.subplots(4, 4) fig.set_size_inches(40, 40) fig.set_dpi(40) itime = integ_time[0:-1] - integ_time[1:] - diff = compute_diff_imas(cube, integ_time, diff_method = diff_method) - - + diff = compute_diff_imas(cube, integ_time, diff_method=diff_method) + for i, ax in enumerate(axarr.reshape(-1)): - if (i < cube.shape[-1]-2): - i=i+1 + if i < cube.shape[-1]-2: + i += 1 - diff_i = diff[:,:,i] - vmin = 0 - vmax = 2 + diff_i = diff[:, :, i] + vmin, vmax = zscale(diff_i) im = ax.imshow(np.abs(diff_i), cmap='Greys_r', origin='lower', - vmin = vmin, vmax = vmax) - ax.set_title(f'$\mu = ${median_diff_full_frame[i]:.2f}±{standard_dev_fullframe[i]:.2f} e-/s', fontsize = 30) - - text = ax.text(50, 500, f'{median_diff_lhs[i]:.3f}\n±\n{standard_dev_lhs[i]:.3f}', color='Orange', fontsize=30) + vmin=vmin, vmax=vmax) + title = fr'$\mu = ${median_diff_full_frame[i]:.2f}±{standard_dev_fullframe[i]:.2f} e-/s' + ax.set_title(title, fontsize=30) + text_lhs = f'{median_diff_lhs[i]:.3f}\n±\n{standard_dev_lhs[i]:.3f}' + text = ax.text(50, 500, text_lhs, color='Orange', fontsize=30) text.set_path_effects([path_effects.Stroke(linewidth=15, foreground='black'), - path_effects.Normal()]) - text = ax.text(700, 500, f'{median_diff_rhs[i]:.3f}\n±\n{standard_dev_rhs[i]:.3f}', color='Orange', fontsize=30) + path_effects.Normal()]) + text_rhs = f'{median_diff_rhs[i]:.3f}\n±\n{standard_dev_rhs[i]:.3f}' + text = ax.text(700, 500, text_rhs, color='Orange', fontsize=30) text.set_path_effects([path_effects.Stroke(linewidth=15, foreground='black'), - path_effects.Normal()]) - text = ax.text(200, 900, f'Ratio = {median_diff_lhs[i]/median_diff_rhs[i]:.2f}', color='White', fontsize=30) + path_effects.Normal()]) + text_ratio = f'Ratio = {median_diff_lhs[i]/median_diff_rhs[i]:.2f}' + text = ax.text(200, 900, text_ratio, color='White', fontsize=30) text.set_path_effects([path_effects.Stroke(linewidth=15, foreground='black'), - path_effects.Normal()]) - text = ax.text(300, 300, f'$\Delta = ${median_diff_lhs[i]-median_diff_rhs[i]:.2f}', color='#32CD32', fontsize=30) + path_effects.Normal()]) + text_delta = fr'$\Delta = ${median_diff_lhs[i]-median_diff_rhs[i]:.2f}' + text = ax.text(300, 300, text_delta, color='#32CD32', fontsize=30) text.set_path_effects([path_effects.Stroke(linewidth=15, foreground='black'), - path_effects.Normal()]) + path_effects.Normal()]) - cbar = plt.colorbar(im, ax = ax) - cbar.ax.tick_params(labelsize = 20) - + cbar = plt.colorbar(im, ax=ax) + cbar.ax.tick_params(labelsize=20) ax.set_yticklabels([]) ax.set_xticklabels([]) - ax.set_xlabel(f'{xlabel_list[i]}, $\Delta t = ${np.abs(itime[i]):.2f} sec', fontsize = 30) - + ax.set_xlabel(fr'{xlabel_list[i]}, $\Delta t = ${np.abs(itime[i]):.2f} sec', fontsize=30) else: - ax.set_axis_off() return fig @@ -322,7 +317,6 @@ def plot_ima_subplots(ima_filename, vmin, vmax): vmax: float Maximum magnitude for scaling the data range that the colormap covers. ''' - path, filename = os.path.split(ima_filename) cube, integ_time = read_wfc3(ima_filename) @@ -330,23 +324,22 @@ def plot_ima_subplots(ima_filename, vmin, vmax): fig_panel1, axarr = plt.subplots(4, 4) fig_panel1.set_size_inches(40, 40) fig_panel1.set_dpi(40) - plt.rcParams.update({'font.size':40}) - itime = integ_time[0:-1] - integ_time[1:] - read_title=np.arange(16,0,-1) + plt.rcParams.update({'font.size': 40}) + read_title = np.arange(16, 0, -1) for i, ax in enumerate(axarr.reshape(-1)): - im = ax.imshow(cube[:,:,i], cmap = 'Greys_r', origin = 'lower', vmin = vmin , vmax = vmax) + im = ax.imshow(cube[:, :, i], cmap='Greys_r', origin='lower', vmin=vmin, vmax=vmax) - cbar=plt.colorbar(im, ax = ax) - cbar.ax.tick_params(labelsize = 20) - ax.set_title(f'SCI, {read_title[i]}', fontsize = 40) + cbar = plt.colorbar(im, ax=ax) + cbar.ax.tick_params(labelsize=20) + ax.set_title(f'SCI, {read_title[i]}', fontsize=40) ax.set_yticklabels([]) ax.set_xticklabels([]) - _=fig_panel1.suptitle(filename, fontsize = 40) - plt.subplots_adjust(bottom = 0.3, right = 0.9, top = 0.95) - - + _ = fig_panel1.suptitle(filename, fontsize=40) + plt.subplots_adjust(bottom=0.3, right=0.9, top=0.95) + + def plot_ramp_subplots(ima_files, difference_method, ylims, exclude_sources, lhs_region, rhs_region): ''' Build a simple figure with subplots of IMA accumulation ramps. @@ -373,40 +366,42 @@ def plot_ramp_subplots(ima_files, difference_method, ylims, exclude_sources, lhs rhs_region : dict The four corners (x0, x1, y0, y1) of the right hand region. ''' - - fig = plt.figure(figsize = (50, 20)) + fig = plt.figure(figsize=(50, 20)) fig rows = 1 columns = 2 subplot_titles = ['scattered', 'nominal'] - for i,ima in enumerate(ima_files): - + for i, ima in enumerate(ima_files): path, filename = os.path.split(ima) cube, integ_time = read_wfc3(ima) - if exclude_sources == True: + if exclude_sources is True: cube[np.abs(cube) > 3] = np.nan - diff_cube = compute_diff_imas(cube, integ_time, diff_method = difference_method) - median_diff_fullframe, median_diff_lhs, median_diff_rhs = get_median_fullframe_lhs_rhs(diff_cube, lhs_region = lhs_region, rhs_region = rhs_region) + diff_cube = compute_diff_imas(cube, integ_time, diff_method=difference_method) + + median_diff_fullframe, median_diff_lhs, median_diff_rhs = ( + get_median_fullframe_lhs_rhs(diff_cube, + lhs_region=lhs_region, + rhs_region=rhs_region)) ax = fig.add_subplot(rows, columns, i+1) plot_ramp(ima, integ_time, median_diff_fullframe, median_diff_lhs, median_diff_rhs) - ax.set_ylim(ylims[0],ylims[1]) + ax.set_ylim(ylims[0], ylims[1]) - ax.tick_params(axis = "x", labelsize = 30) - ax.tick_params(axis = "y", labelsize = 30) + ax.tick_params(axis="x", labelsize=30) + ax.tick_params(axis="y", labelsize=30) - _=ax.set_title(f'{filename}, {subplot_titles[i]}', fontsize=50) + _ = ax.set_title(f'{filename}, {subplot_titles[i]}', fontsize=50) def plot_ima_difference_subplots(ima_filename, difference_method, lhs_region, rhs_region): ''' Build a complex panel plot of the difference between individual IMA reads. - The median difference $\mu$ in the count rate over the entire image is printed above each panel. Below each panel, + The median difference mu in the count rate over the entire image is printed above each panel. Below each panel, The IMSET difference, along with the time interval between the two IMSETs, is printed below. The statistics in orange (on the left and right side of each panel) give the median rate and standard deviation of each side of the image, respectively. The value in green 'delta' is the @@ -430,17 +425,30 @@ def plot_ima_difference_subplots(ima_filename, difference_method, lhs_region, rh ''' - path,filename = os.path.split(ima_filename) + path, filename = os.path.split(ima_filename) cube, integ_time = read_wfc3(ima_filename) - median_fullframe, median_lhs, median_rhs = get_median_fullframe_lhs_rhs(cube, lhs_region = lhs_region, rhs_region = rhs_region) + median_fullframe, median_lhs, median_rhs = ( + get_median_fullframe_lhs_rhs(cube, + lhs_region=lhs_region, + rhs_region=rhs_region)) - diff_cube = compute_diff_imas(cube, integ_time, diff_method = difference_method) + diff_cube = compute_diff_imas(cube, integ_time, diff_method=difference_method) - median_diff_fullframe, median_diff_lhs, median_diff_rhs = get_median_fullframe_lhs_rhs(diff_cube, lhs_region = lhs_region, rhs_region = rhs_region) - standard_dev_fullframe, standard_dev_lhs, standard_dev_rhs = get_std_fullframe_lhs_rhs(diff_cube, lhs_region = lhs_region, rhs_region = rhs_region) + median_diff_fullframe, median_diff_lhs, median_diff_rhs = ( + get_median_fullframe_lhs_rhs(diff_cube, + lhs_region=lhs_region, + rhs_region=rhs_region)) + + standard_dev_fullframe, standard_dev_lhs, standard_dev_rhs = ( + get_std_fullframe_lhs_rhs(diff_cube, + lhs_region=lhs_region, + rhs_region=rhs_region)) - fig_0 = panel_plot(cube, integ_time, median_diff_fullframe, median_diff_lhs, median_diff_rhs, standard_dev_fullframe, standard_dev_lhs, standard_dev_rhs, diff_method = difference_method) - _=fig_0.suptitle(filename, fontsize = 40) - plt.subplots_adjust(bottom = 0.25, right = 0.9, top = 0.95) \ No newline at end of file + fig_0 = panel_plot(cube, integ_time, median_diff_fullframe, median_diff_lhs, + median_diff_rhs, standard_dev_fullframe, standard_dev_lhs, + standard_dev_rhs, diff_method=difference_method) + + _ = fig_0.suptitle(filename, fontsize=40) + plt.subplots_adjust(bottom=0.25, right=0.9, top=0.95) diff --git a/notebooks/WFC3/ir_scattered_light_manual_corrections/pre-requirements.sh b/notebooks/WFC3/ir_scattered_light_manual_corrections/pre-requirements.sh new file mode 100644 index 000000000..23a5dffae --- /dev/null +++ b/notebooks/WFC3/ir_scattered_light_manual_corrections/pre-requirements.sh @@ -0,0 +1 @@ +conda install --yes -c conda-forge hstcal \ No newline at end of file diff --git a/notebooks/WFC3/ir_scattered_light_manual_corrections/requirements.txt b/notebooks/WFC3/ir_scattered_light_manual_corrections/requirements.txt index 2814a9434..a742650f7 100644 --- a/notebooks/WFC3/ir_scattered_light_manual_corrections/requirements.txt +++ b/notebooks/WFC3/ir_scattered_light_manual_corrections/requirements.txt @@ -1,8 +1,9 @@ astropy==5.2.1 astroquery==0.4.6 drizzlepac==3.5.1 -ginga==4.1.1 +ginga==4.0.1 matplotlib==3.7.0 numpy==1.23.4 stwcs==1.7.2 wfc3tools==1.4.0 +crds==11.17.9