From e125dbcfab154f8845a1e31a6a936d8e7bc76cfb Mon Sep 17 00:00:00 2001 From: Michael Dulude Date: Thu, 11 Jan 2024 08:52:06 -0500 Subject: [PATCH] Add WFC3 notebook 'Correcting_for_Scattered_Light_in_IR_Exposures_Using_calwf3_to_Mask_Bad_Reads.ipynb' (#107) * PEP 8 corrections, updating pinned version of ginga in requirements.txt, adding a pre-requirements.sh * added crds to requirements.txt, one last PEP8 correction to script * Update ima_visualization_and_differencing.py added an empty line at the end of the script... * Some changes based on notebook review * Update links to hst notebooks repository and fix get data * Fix style errors --------- Co-authored-by: annierose3 Co-authored-by: annierose3 <112646956+annierose3@users.noreply.github.com> Co-authored-by: Hatice Karatay --- _config.yml | 1 - _toc.yml | 4 +- ...sures_Using_calwf3_to_Mask_Bad_Reads.ipynb | 172 +++++++-------- .../ima_visualization_and_differencing.py | 206 +++++++++--------- .../pre-requirements.sh | 1 + .../requirements.txt | 3 +- 6 files changed, 195 insertions(+), 192 deletions(-) create mode 100644 notebooks/WFC3/ir_scattered_light_calwf3_corrections/pre-requirements.sh diff --git a/_config.yml b/_config.yml index 834af357d..938404b26 100644 --- a/_config.yml +++ b/_config.yml @@ -52,5 +52,4 @@ exclude_patterns: [notebooks/DrizzlePac/align_mosaics/align_mosaics.ipynb, notebooks/WFC3/dash/dash.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/zeropoints/zeropoints.ipynb] diff --git a/_toc.yml b/_toc.yml index ee7023e90..87a784f18 100644 --- a/_toc.yml +++ b/_toc.yml @@ -69,7 +69,9 @@ parts: # - file: notebooks/WFC3/flux_conversion_tool/flux_conversion_tool.ipynb - 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/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/persistence/wfc3_ir_persistence.ipynb - file: notebooks/WFC3/tvb_flattenramp/TVB_flattenramp_notebook.ipynb diff --git a/notebooks/WFC3/ir_scattered_light_calwf3_corrections/Correcting_for_Scattered_Light_in_IR_Exposures_Using_calwf3_to_Mask_Bad_Reads.ipynb b/notebooks/WFC3/ir_scattered_light_calwf3_corrections/Correcting_for_Scattered_Light_in_IR_Exposures_Using_calwf3_to_Mask_Bad_Reads.ipynb index 666dff9a5..85d37744a 100644 --- a/notebooks/WFC3/ir_scattered_light_calwf3_corrections/Correcting_for_Scattered_Light_in_IR_Exposures_Using_calwf3_to_Mask_Bad_Reads.ipynb +++ b/notebooks/WFC3/ir_scattered_light_calwf3_corrections/Correcting_for_Scattered_Light_in_IR_Exposures_Using_calwf3_to_Mask_Bad_Reads.ipynb @@ -22,7 +22,7 @@ "- Compare the original DRZ to the reprocessed DRZ image.\n", "\n", "A second method (manually subtracting bad reads from the final IMA read) can \n", - "be found in the notebook [Correcting for Scattered Light in WFC3/IR Exposures: Manually Subtracting Bad Reads (O'Connor 2023)](https://github.com/spacetelescope/WFC3Library/tree/master/notebooks/ir_scattered_light_manual_corrections). This provides a method for removing anomalies such as satellite trails which appear in the middle of an IMA exposure. One caveat of the second method is that it does not perform the 'ramp fitting' step and therefore the calibrated FLT products will still contain cosmic rays. \n", + "be found in the notebook [Correcting for Scattered Light in WFC3/IR Exposures: Manually Subtracting Bad Reads (O'Connor 2023)](https://github.com/spacetelescope/hst_notebooks/tree/main/notebooks/WFC3/ir_scattered_light_manual_corrections). This provides a method for removing anomalies such as satellite trails which appear in the middle of an IMA exposure. One caveat of the second method is that it does not perform the 'ramp fitting' step and therefore the calibrated FLT products will still contain cosmic rays. \n", "\n", "
\n", "\n", @@ -68,7 +68,7 @@ "have a reduced total exposure time, given the rejection of some number of reads, and \n", "therefore a lower signal-to-noise ratio.\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 in 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/blob/main/notebooks/WFC3/ir_ima_visualization/IR_IMA_Visualization_with_an_Example_of_Time_Variable_Background.ipynb) for a walkthrough of how to identify a TVB in due to scattered light. " ] }, { @@ -117,7 +117,7 @@ "\n", "import ima_visualization_and_differencing as diff\n", "\n", - "%matplotlib inline\n" + "%matplotlib inline" ] }, { @@ -142,16 +142,15 @@ "metadata": {}, "outputs": [], "source": [ - "OBS_ID = 'ICQTBB020' #Earth-limb at the start \n", - "#OBS_ID = 'ICXT27020' #Earth-limb at the end \n", + "OBS_ID = 'ICQTBB020' # Earth-limb at the start \n", + "# OBS_ID = 'ICXT27020' # Earth-limb at the end \n", "data_list = Observations.query_criteria(obs_id=OBS_ID)\n", "\n", "file_id = \"icqtbbbxq\"\n", - "#file_id = 'icxt27hkq'\n", - "Observations.download_products(data_list['obsid'], project='CALWF3', obs_id=file_id, \n", - " mrp_only=False, productSubGroupDescription=['RAW', 'IMA', 'FLT'])\n", - "\n", - "\n" + "# file_id = 'icxt27hkq'\n", + "Observations.download_products(data_list['obsid'], project='CALWF3', \n", + " obs_id=file_id, mrp_only=False, \n", + " productSubGroupDescription=['RAW', 'IMA', 'FLT'])" ] }, { @@ -175,7 +174,7 @@ "\n", "raw_file = f'mastDownload/HST/{file_id}/{file_id}_raw.fits'\n", "\n", - "remove_files_list=glob.glob(f'./{file_id}_*.fits')\n", + "remove_files_list = glob.glob(f'./{file_id}_*.fits')\n", "for rm_file in remove_files_list:\n", " os.remove(rm_file)\n", "\n", @@ -193,7 +192,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "In this section, we show how to identify the reads impacted by TVB 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 in WFC3/IR images. \n", + "In this section, we show how to identify the reads impacted by TVB 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/blob/main/notebooks/WFC3/ir_ima_visualization/IR_IMA_Visualization_with_an_Example_of_Time_Variable_Background.ipynb) notebook, which includes a more comprehensive walkthrough of identifying time variable background in WFC3/IR images. \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", @@ -208,33 +207,35 @@ "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", "path, filename = os.path.split(ima_filepath)\n", "cube, integ_time = diff.read_wfc3(ima_filepath)\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", - "#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", + "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(\n", + " diff_cube, \n", + " lhs_region=lhs_region, \n", + " rhs_region=rhs_region))\n", "\n", - "plt.rc('xtick', labelsize = 20) \n", - "plt.rc('ytick', labelsize = 20) \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", "\n", - "plt.ylim(0.5,2.5)\n", - "_ = plt.title(filename)\n", - "\n" + "plt.ylim(0.5, 2.5)\n", + "_ = plt.title(filename)" ] }, { @@ -259,14 +260,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", + "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)" + "diff.plot_ima_difference_subplots(ima_filepath, \n", + " difference_method='instantaneous', \n", + " lhs_region=lhs_region, \n", + " rhs_region=rhs_region)" ] }, { @@ -338,8 +342,7 @@ "source": [ "raw_file = f'{file_id}_raw.fits'\n", "\n", - "\n", - "print( f\"Querying CRDS for the reference file associated with {raw_file}.\")\n", + "print(f\"Querying CRDS for the reference file associated with {raw_file}.\")\n", "command_line_input = 'crds bestrefs --files {:} --sync-references=1 --update-bestrefs'.format(raw_file)\n", "os.system(command_line_input)" ] @@ -357,7 +360,7 @@ "source": [ "As discussed in the [introduction](#intro) to this notebook, the accuracy of the ramp fit performed during the CRCORR cosmic-ray rejection step of the pipeline determines the quality of the calibrated WFC3/IR FLT data products. Given that a time variable background can compromise the quality of the ramp fit, observations affected by TVB will likely result in poor-quality calibrated FLT images (see the Appendix of [WFC3 ISR 2021-01](https://www.stsci.edu/files/live/sites/www/files/home/hst/instrumentation/wfc3/documentation/instrument-science-reports-isrs/_documents/2021/2021-01.pdf)).\n", "\n", - "To address poorly calibrated FLT images where some reads are affected by scattered light TVB (as in our example observation), we can remove these reads and re-run `calwf3` to produce cleaner FLT images. 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 poorly calibrated FLT images where some reads are affected by scattered light TVB (as in our example observation), we can remove these reads and re-run `calwf3` to produce cleaner FLT images. 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", "\n", "The following reprocessing example is replacing section [7.10.1 of the WFC3 Data Handbook](https://hst-docs.stsci.edu/wfc3dhb/chapter-7-wfc3-ir-sources-of-error/7-10-time-variable-background). " ] @@ -393,11 +396,11 @@ " if os.path.exists(filepath):\n", " os.remove(filepath)\n", "\n", - "reads = np.arange(11,16) #numpy arange creates an array including the start value (11) but excluding the stop value (16), so the array is actually 11-15. \n", + "reads = np.arange(11, 16) # numpy arange creates an array including the start value (11) but excluding the stop value (16), so the array is actually 11-15. \n", "for read in reads:\n", " fits.setval(raw_filepath, extver=read, extname='DQ', keyword='pixvalue', value=1024)\n", "\n", - "calwf3(raw_filepath)\n" + "calwf3(raw_filepath)" ] }, { @@ -436,9 +439,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "scrolled": false - }, + "metadata": {}, "outputs": [], "source": [ "image_new = fits.getdata(reprocessed_flt)\n", @@ -449,25 +450,23 @@ "rows = 1\n", "columns = 2\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 FLT image\", fontsize=20)\n", "im1 = plt.imshow(image_new, vmin=0.8, vmax=1.3, origin='lower', cmap='Greys_r')\n", - "ax1.tick_params(axis='both', labelsize = 10)\n", + "ax1.tick_params(axis='both', labelsize=10)\n", "cbar1 = plt.colorbar(im1, ax=ax1)\n", "cbar1.ax.tick_params(labelsize=10)\n", "\n", - "\n", - "ax2 = fig.add_subplot(rows,columns,2)\n", + "ax2 = fig.add_subplot(rows, columns, 2)\n", "ax2.set_title(\"Original FLT image\", fontsize=20)\n", - "im2 = plt.imshow(image_old, vmin=0.8, vmax=1.3, origin='lower', cmap='Greys_r')\n", + "im2 = plt.imshow(image_old, 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", - "\n" + "plt.rc('ytick', labelsize=10) " ] }, { @@ -478,7 +477,7 @@ "\n", "Finally, note that the reprocessed FLT product now includes a larger sky background in science pixels corresponding to IR \"blobs\", regions of reduced sensitivity due particulate matter on the channel selection mechanism (CSM) mirror ([section 7.5 of the WFC3 Data Handbook](https://hst-docs.stsci.edu/wfc3dhb/chapter-7-wfc3-ir-sources-of-error/7-5-blobs)). \n", "\n", - "The method for correcting WFC3/IR images described in the notebook [Correcting for Scattered Light in WFC3/IR Exposures: Manually Subtracting Bad Reads (O'Connor 2023)](https://github.com/spacetelescope/WFC3Library/tree/master/notebooks/ir_scattered_light_manual_corrections) provides FLT products without blobs, but that do include cosmic rays.\n", + "The method for correcting WFC3/IR images described in the notebook [Correcting for Scattered Light in WFC3/IR Exposures: Manually Subtracting Bad Reads (O'Connor 2023)](https://github.com/spacetelescope/hst_notebooks/tree/main/notebooks/WFC3/ir_scattered_light_manual_corrections) provides FLT products without blobs, but that do include cosmic rays.\n", "\n", "We update the FLT header keywords \"EXPTIME\" (in the primary header) and \"SAMPTIME\" (in the science header) to reflect the new total exposure time." ] @@ -491,32 +490,28 @@ }, "outputs": [], "source": [ - "image_new = fits.open(reprocessed_flt, mode ='update')\n", - "\n", - "ima_orig = fits.open(original_ima)\n", + "with fits.open(reprocessed_flt, mode='update') as image_new, fits.open(original_ima) as ima_orig:\n", " \n", - "hdr = ima_orig[0].header\n", - "NSAMP = hdr['NSAMP']\n", - "hdr1 = ima_orig[1].header\n", - "integ_time = np.zeros(shape=(NSAMP))\n", - "for i in range(1, NSAMP+1):\n", + " hdr = ima_orig[0].header\n", + " NSAMP = hdr['NSAMP']\n", + " hdr1 = ima_orig[1].header\n", + " integ_time = np.zeros(shape=(NSAMP))\n", + " for i in range(1, NSAMP+1):\n", " integ_time[i-1] = ima_orig[('TIME', i)].header['PIXVALUE']\n", - "integ_time = integ_time[::-1]\n", "\n", - "dt = np.diff(integ_time)\n", + " integ_time = integ_time[::-1]\n", "\n", - "final_time = integ_time[-1]\n", + " dt = np.diff(integ_time)\n", "\n", - "if (len(reads) > 0):\n", - " for read in reads:\n", - " index = len(integ_time)-read-1 #because the reads are stored in reverse order\n", - " final_time -= dt[index]\n", - " \n", - " \n", - "print(f'The final exposure time after reprocessing is {final_time}.')\n", - "image_new[0].header['EXPTIME'] = final_time\n", - "image_new.flush()\n", - "image_new.close()" + " final_time = integ_time[-1]\n", + "\n", + " if (len(reads) > 0):\n", + " for read in reads:\n", + " index = len(integ_time)-read-1 # because the reads are stored in reverse order\n", + " final_time -= dt[index]\n", + "\n", + " print(f'The final exposure time after reprocessing is {final_time}.')\n", + " image_new[0].header['EXPTIME'] = final_time" ] }, { @@ -543,8 +538,10 @@ "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'])\n" + "Observations.download_products(data_list['obsid'], \n", + " project='CALWF3',\n", + " mrp_only=False, \n", + " productSubGroupDescription=['FLT', 'DRZ'])" ] }, { @@ -554,7 +551,7 @@ "outputs": [], "source": [ "nominal_file_ids = [\"icqtbbc0q\", \"icqtbbbrq\", \"icqtbbbtq\"]\n", - "#nominal_file_ids = [\"icxt27hoq\"]\n", + "# nominal_file_ids = [\"icxt27hoq\"]\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", @@ -592,7 +589,7 @@ "outputs": [], "source": [ "astrodrizzle.AstroDrizzle('icqtbb*flt.fits', output='f140w', preserve=False, build=False, context=False, skymethod='match', driz_separate=False, median=False, blot=False, driz_cr=False, final_bits='16')\n", - "#astrodrizzle.AstroDrizzle('icxt27*flt.fits', output='f105w', preserve=False, build=False, context=False, skymethod='match', driz_separate=False, median=False, blot=False, driz_cr=False, final_bits='16')" + "# astrodrizzle.AstroDrizzle('icxt27*flt.fits', output='f105w', preserve=False, build=False, context=False, skymethod='match', driz_separate=False, median=False, blot=False, driz_cr=False, final_bits='16')" ] }, { @@ -610,27 +607,22 @@ }, "outputs": [], "source": [ - "DRZ_image = fits.open(\"f140w_drz_sci.fits\")\n", - "Orig_DRZ = fits.open('mastDownload/HST/icqtbb020/icqtbb020_drz.fits')\n", + "with fits.open(\"f140w_drz_sci.fits\") as DRZ_image, fits.open('mastDownload/HST/icqtbb020/icqtbb020_drz.fits') as Orig_DRZ:\n", "\n", - "#DRZ_image = fits.getdata(\"f105w_drz_sci.fits\")\n", - "#Orig_DRZ = fits.getdata('mastDownload/HST/icxt27020/icxt27020_drz.fits')\n", + " fig = plt.figure(figsize=(20, 10))\n", + " rows = 1\n", + " columns = 2\n", "\n", - "fig = plt.figure(figsize=(20,10))\n", - "rows = 1\n", - "columns = 2\n", + " ax1 = fig.add_subplot(rows, columns, 1)\n", + " ax1.set_title(\"Reprocessed DRZ Image\", fontsize=20)\n", + " vmin, vmax = zscale(Orig_DRZ[1].data)\n", + " im1 = plt.imshow(DRZ_image[0].data, vmin=vmin, vmax=vmax, origin='lower', cmap='Greys_r')\n", + " _ = plt.colorbar()\n", "\n", - "\n", - "ax1 = fig.add_subplot(rows, columns, 1)\n", - "ax1.set_title(\"Reprocessed DRZ Image\", fontsize=20)\n", - "vmin,vmax = zscale(Orig_DRZ[1].data)\n", - "im1 = plt.imshow(DRZ_image[0].data, vmin=vmin, vmax=vmax, origin='lower', cmap='Greys_r')\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[1].data, vmin=vmin, vmax=vmax, origin='lower', cmap='Greys_r')\n", - "_= plt.colorbar()\n" + " ax2 = fig.add_subplot(rows, columns, 2)\n", + " ax2.set_title(\"Original Pipeline DRZ Image\", fontsize=20)\n", + " im2 = plt.imshow(Orig_DRZ[1].data, vmin=vmin, vmax=vmax, origin='lower', cmap='Greys_r')\n", + " _ = plt.colorbar()" ] }, { @@ -666,7 +658,7 @@ "\n", "**Author:** Anne O'Connor, Jennifer Mack, Annalisa Calamida, Harish Khandrika -- WFC3 Instrument\n", "\n", - "**Updated On:** 2023-05-25\n", + "**Updated On:** 2023-12-26\n", "\n", "## Citations \n", "\n", @@ -716,9 +708,9 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.13" + "version": "3.11.6" } }, "nbformat": 4, - "nbformat_minor": 2 + "nbformat_minor": 4 } diff --git a/notebooks/WFC3/ir_scattered_light_calwf3_corrections/ima_visualization_and_differencing.py b/notebooks/WFC3/ir_scattered_light_calwf3_corrections/ima_visualization_and_differencing.py index 62f6ee235..d43e744b7 100644 --- a/notebooks/WFC3/ir_scattered_light_calwf3_corrections/ima_visualization_and_differencing.py +++ b/notebooks/WFC3/ir_scattered_light_calwf3_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) + 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_calwf3_corrections/pre-requirements.sh b/notebooks/WFC3/ir_scattered_light_calwf3_corrections/pre-requirements.sh new file mode 100644 index 000000000..23a5dffae --- /dev/null +++ b/notebooks/WFC3/ir_scattered_light_calwf3_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_calwf3_corrections/requirements.txt b/notebooks/WFC3/ir_scattered_light_calwf3_corrections/requirements.txt index 2814a9434..a742650f7 100644 --- a/notebooks/WFC3/ir_scattered_light_calwf3_corrections/requirements.txt +++ b/notebooks/WFC3/ir_scattered_light_calwf3_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