From 4e18ec8818f208d4ee60c32856a1822b7bea3156 Mon Sep 17 00:00:00 2001 From: Sierra Gomez Date: Fri, 12 Jul 2024 12:37:19 -0400 Subject: [PATCH 01/26] Add HASP Wavelength Adjustment notebook --- .../WavelengthAdjustment.ipynb | 2624 +++++++++++++++++ .../WavelengthAdjustment/pre-requirements.sh | 1 + .../WavelengthAdjustment/requirements.txt | 11 + 3 files changed, 2636 insertions(+) create mode 100644 notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb create mode 100755 notebooks/HASP/WavelengthAdjustment/pre-requirements.sh create mode 100644 notebooks/HASP/WavelengthAdjustment/requirements.txt diff --git a/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb b/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb new file mode 100644 index 000000000..8e8d7eb51 --- /dev/null +++ b/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb @@ -0,0 +1,2624 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "bdd3127b-f10f-4a3a-bdc3-8e96290db8a9", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "# Wavelength Adjustments for Running the Hubble Advanced Spectral Products (HASP) Script\n", + "\n", + "### This Notebook is designed to walk you through making wavelength adjustments to COS & STIS data when running the **Hubble Advanced Spectral Products (HASP)** coadd script.\n", + "\n", + "## Learning Goals:\n", + "\n", + "By the end of this tutorial, you will learn how to:\n", + "\n", + "- Determine if you must make wavelength adjustments based on your coadded product\n", + "\n", + "- Perform cross-correlation on `X1D` and `SX1` files to determine the wavelength offset between datasets\n", + "\n", + "- Perform wavelength shifts on COS and STIS data using each instrument's adjustment methods\n", + "\n", + "- Run the coadd code on new data and compare pre- and post-shift coadded spectra\n", + "\n", + "## Table of Contents\n", + "**0. [Introduction](#introduction)**\n", + "\n", + "**1. [Performing Wavelength Adjustments on COS Data](#maincos)**\n", + "\n", + "\\- 1.1 [Downloading and Organizing COS Data using `Astroquery`](#cosdownload)\n", + "\n", + "\\- 1.2 [Running the HASP Script on COS Data](#scriptrun1_cos)\n", + "\n", + "\\- 1.3 [Identifying a Wavelength Shift in COS Data](#id_shift)\n", + "\n", + "\\- 1.4 [Performing Cross-Correlation on `X1D` COS data and Calculating Lag](#cc_cos)\n", + "\n", + "\\- 1.5 [Applying Shifts to COS data and Re-creating Coadd](#apply_shifts_cos)\n", + "\n", + "\\- 1.6 [Comparing Non-Shifted and Shifted COS Coadds](#compare_shifted_cos)\n", + "\n", + "**2. [Performing Wavelength Adjustments on STIS data](#mainstis)**\n", + "\n", + "\\- 2.1 [Downloading and Organizing STIS Data using `Astroquery`](#downloadstis)\n", + "\n", + "\\- 2.2 [Running the HASP Script on STIS Data](#initial_coadd_stis)\n", + "\n", + "\\- 2.3 [Performing Cross-Correlation on STIS data and Calculating Lag](#cross_correlate_stis)\n", + "\n", + "\\- 2.4 [Applying Shifts to Raw STIS data and Re-creating Coadd](#run_calstis)\n", + "\n", + "\\- 2.5 [Comparing Non-Adjusted and Adjusted STIS Coadds](#compare_shifted_stis)\n", + "\n", + "\n", + "## 0. Introduction\n", + "\n", + "The [Hubble Advanced Spectral Products (HASP) `coadd` code](https://github.com/spacetelescope/hasp) is a script that coadds spectra of the same target within a program. This software can coadd data taken with the spectrographs onboard the [Hubble Space Telescope (HST)](https://www.stsci.edu/hst): the [Space Telescope Imaging Spectrograph (STIS)](https://www.stsci.edu/hst/instrumentation/stis) and the [Cosmic Origins Spectrograph (COS)](https://www.stsci.edu/hst/instrumentation/cos). The [Hubble Spectroscopic Legacy Archive (HSLA)](https://archive.stsci.edu/missions-and-data/hst/hasp) uses this script to coadd these instruments’ data from [The Mikulski Archive for Space Telescopes (MAST)](https://archive.stsci.edu/) to create high-quality spectra with a broad wavelength coverage that is publicly available for the scientific community. The script first coadds the observations for each grating for a given visit, then it combines all gratings for the observation set. Finally, it coadds the spectra of each observation set in the program to produce a fully coadded spectrum for each target in a program. [Check out the COS 2024-1 ISR for more information about the HASP script and products](https://www.stsci.edu/files/live/sites/www/files/home/hst/instrumentation/cos/documentation/instrument-science-reports-isrs/_documents/ISR2024-01.pdf).\n", + "\n", + "Wavelength alignment between exposures is crucial for performing accurate spectroscopic analysis on your data. Wavelength offsets can occur for a variety of reasons, such as poor centering during target acquisition or thermal fluctuations in between exposures. In addition to inaccurate measurements and the misidentification of spectral features, incorrect wavelength registration can result in the wrong flux sensitivity being applied at each wavelength. For COS in particular, the wavelength throughput decreases for low-resolution spectra and rapidly at shorter wavelengths, thus being more affected by incorrect flux sensitivity applications. It is sometimes necessary for users to manually apply shifts to their raw data to ensure wavelength alignment between files. More information about the wavelength accuracies for the instruments can be found in [Chapter 5](https://hst-docs.stsci.edu/cosihb/chapter-5-spectroscopy-with-cos/5-1-the-capabilities-of-cos#id-5.1TheCapabilitiesofCOS-Section5.1.115.1.11WavelengthAccuracy) of the COS Instrument Handbook and [Chapter 4](https://hst-docs.stsci.edu/stisdhb/chapter-4-stis-error-sources/4-2-summary-of-accuracies) in the STIS Instrument Handbook.\n", + "\n", + "This notebook will show users how to recognize if their data needs wavelength adjustments, and how to perform these corrections for both COS and STIS data. The example datasets that are used in this tutorial are abnormal -- most cases of wavelength shifts are small. It assumes that the user knows how to download MAST data and how to run the coadd code. Check out our [CoaddTutorial.ipynb](https://github.com/spacetelescope/hst_notebooks/tree/main/notebooks/HASP/CoaddTutorial) notebook for a detailed explanation of how to do both of these things.\n", + "\n", + "**Please check out our [Setup.ipynb](https://github.com/spacetelescope/hst_notebooks/tree/main/notebooks/HASP/Setup) notebook before running this tutorial to learn how to install and run the coadd code.**\n", + "\n", + "### Imports\n", + "\n", + "We will be using multiple libraries to retrieve and analyze data. We will use: \n", + "* `os`, `shutil`, `Path.pathlib`, and `glob` to work with our system's directories and filepaths\n", + "* `matplotlib.pyplot` to plot our data and shifts\n", + "* `numpy` to perform analysis on our data\n", + "* `astroquery.mast` to download our STIS and COS datasets for the exercises\n", + "* `astropy.io fits` to read the corresponding `FITS` files\n", + "* `calcos` and `stistools` to run the COS and STIS calibration pipelines, respectively\n", + "* `scipy.optimize` to fit a Gaussian when determining the radial velocity\n", + "* `astropy.modeling` and `scipy.signal` to help us fit a quadratic when calculating our shifts\n", + "\n", + "We recommend creating a HASP-specific `conda` environment when coadding spectra. You can check out our [Setup.ipynb](https://github.com/spacetelescope/hst_notebooks/tree/main/notebooks/HASP/Setup) notebook to create such an environment. Alternatively, you can also download the required dependencies to run this notebook with the terminal command:\n", + "```\n", + "pip install -r requirements.txt\n", + "```\n", + "This will download the dependencies that are necessary to run this current notebook. You will also need to install `hstcal` to run the STIS portion of this notebook. You can run the command below in your terminal after you've activated your `conda` environment:\n", + "```\n", + "conda install -y -c conda-forge hstcal==2.7.4\n", + "```\n", + "Or, similar to the `requirements.txt` file, you can run:\n", + "```\n", + "bash pre-requirements.sh\n", + "```\n", + "Let's import all of the packages that we will use in this notebook and print our `conda` environment by running the next cell:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "aa6a80d5-7b6c-49f0-8d68-94150cf9b08a", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "import os\n", + "import shutil\n", + "from pathlib import Path\n", + "\n", + "import glob\n", + "\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "\n", + "from astroquery.mast import Observations\n", + "from astropy.io import fits\n", + "\n", + "import calcos\n", + "import stistools\n", + "\n", + "from astropy.modeling import fitting\n", + "from astropy.modeling.models import Polynomial1D\n", + "\n", + "from scipy.signal import correlate\n", + "from scipy.signal import correlation_lags\n", + "\n", + "from scipy import optimize\n", + "\n", + "print(\"Currently active conda environment:\", os.environ.get(\"CONDA_PREFIX\"))" + ] + }, + { + "cell_type": "markdown", + "id": "500891d8", + "metadata": {}, + "source": [ + "Since we will be running both `CalCOS` and `CalSTIS`, we will need to set up environment variables for when we use the [Calibration Reference Data System](https://hst-crds.stsci.edu/). Uncomment these lines if you have not set up environment variables for both COS and STIS:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1a16bb3d", + "metadata": {}, + "outputs": [], + "source": [ + "# os.environ['CRDS_SERVER_URL'] = 'https://hst-crds.stsci.edu'\n", + "# os.environ['CRDS_SERVER'] = 'https://hst-crds.stsci.edu'\n", + "# os.environ['CRDS_PATH'] = './crds_cache'\n", + "# os.environ['lref'] = './crds_cache/references/hst/cos/'\n", + "# os.environ['oref'] = './crds_cache/references/hst/stis/'" + ] + }, + { + "cell_type": "markdown", + "id": "278b7b9a-8b7b-42a0-bfc8-fba3ab26b3a4", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "We will create a function to consolidate and organize our data since we will be utilizing it for two datasets in the tutorial. The coadd script requires all `X1D` and `SX1` files to be in a single directory to be coadded. When downloading data from MAST, it creates multiple directories: `./mastDownload/HST` and within this directory are sub-directories for each observation set. We will move all of the files in these sub-directories to the respective STIS/COS directory." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "243598dd-7d3b-4bde-ae16-dc6fa5f6bee6", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "def consolidate_files(data_path):\n", + " '''\n", + " Consolidate all files to single directory; necessary for HASP script run.\n", + " ---------------\n", + " Input:\n", + " str data_path : ./mastDownload/HST folders paths; files to be moved here\n", + " ---------------\n", + " Output:\n", + " None. Files moved to data_path. ./mastDownload/HST directory is deleted.\n", + " '''\n", + " # The path to all obs_id folders\n", + " mast_path = f\"{data_path}/mastDownload/HST/\"\n", + "\n", + " try:\n", + " # Check if mastDownload exists\n", + " if not os.path.exists(mast_path):\n", + " print(f\"Directory {mast_path} doesn't exist.\")\n", + " return\n", + "\n", + " # Get a list of the obs_id paths in mastDownload\n", + " obs_id_dirs = os.listdir(mast_path)\n", + "\n", + " # Iterate through each obs_id folder and move the files\n", + " for obs_id in obs_id_dirs:\n", + " obs_id_path = os.path.join(mast_path, obs_id)\n", + " files = glob.glob(obs_id_path + \"/*fits\")\n", + "\n", + " for file in files:\n", + " file_path = Path(file)\n", + " new_path = data_path / file_path.name\n", + " shutil.move(file, new_path)\n", + "\n", + " # Now we can remove the mastDownload directory\n", + " if os.path.exists(mast_path):\n", + " shutil.rmtree(f\"{data_path}/mastDownload\")\n", + "\n", + " except Exception as e:\n", + " print(f\"An error occurred: {e}\")" + ] + }, + { + "cell_type": "markdown", + "id": "fac0fdca-c47a-41d9-a774-855105d4e04d", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "**Note: In the cells using `astroquery`, `CalSTIS`, and `CalCOS`, you will see the Jupyter magic command `%%capture`, which keep the output of the cell in memory rather than printing it. This is to reduce the length of the `HTML` notebook. If you would like to see the output of these cells, simply delete the `%%capture` line for each relevant cell and re-run.**\n", + "\n", + "\n", + "## 1. Performing Wavelength Adjustments on COS Data\n", + "\n", + "During science exposures, the wavelength calibration aperture is illuminated by two Pt-Ne lamps, which are used as a template by `CalCOS` to assign wavelengths to positions on the COS detector. A zero-point shift is then applied to the science data to align the wavelength calibration exposure and the science exposure. These wavelength calibration exposures are taken simultaneously with science exposures (except LP6 exposures; more information about LP6 wavecals in [Chapter 5.7.6](https://hst-docs.stsci.edu/cosihb/chapter-5-spectroscopy-with-cos/5-7-internal-wavelength-calibration-exposures#id-5.7InternalWavelengthCalibrationExposures-Section5.7.65.7.6SPLITWavecals(defaultnon-concurrentwavelengthcalibrationatLP6)) of the COS Instrument Handbook). Information about COS wavelength accuracies can be found in [Chapter 5.1.11](https://hst-docs.stsci.edu/cosihb/chapter-5-spectroscopy-with-cos/5-1-the-capabilities-of-cos#id-5.1TheCapabilitiesofCOS-Section5.1.115.1.11WavelengthAccuracy) of the COS Instrument Handbook.\n", + "\n", + "\n", + "We will create a few different directories for our COS example." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4d2ffdc4-05c3-4a8d-8626-9778a36f8148", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "# To hold all of our COS example's output and data\n", + "cos_ex = Path(\"./cos_ex/\")\n", + "\n", + "# Contains the RAWTAG and ASN files\n", + "cos_data = Path(\"./cos_ex/cos_data/\")\n", + "\n", + "# Holds the unshifted X1Ds and the respective coadd\n", + "cos_unshifted = Path(\"./cos_ex/cos_unshifted/\")\n", + "cos_unshifted_coadd = Path(\"./cos_ex/cos_unshifted/coadd/\")\n", + "\n", + "# Will contain the X1Ds created with the shift applied, and respective coadd\n", + "cos_shifted = Path(\"./cos_ex/cos_shifted/\")\n", + "cos_shifted_coadd = Path(\"./cos_ex/cos_shifted/coadd/\")\n", + "\n", + "# If the directory doesn't exist, then create it\n", + "cos_ex.mkdir(exist_ok=True)\n", + "\n", + "cos_data.mkdir(exist_ok=True)\n", + "\n", + "cos_unshifted.mkdir(exist_ok=True)\n", + "cos_unshifted_coadd.mkdir(exist_ok=True)\n", + "\n", + "cos_shifted.mkdir(exist_ok=True)\n", + "cos_shifted_coadd.mkdir(exist_ok=True)" + ] + }, + { + "cell_type": "markdown", + "id": "00fc0317-6144-4c73-bd72-dd35fe71aed4", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "We will be shifting COS FUV spectra of the target [LMC079-1](http://simbad.cds.unistra.fr/simbad/sim-id?Ident=%409822397&Name=%5bM2002%5d%20LMC%20%2071747&submit=submit), a [Wolf-Rayet star](https://en.wikipedia.org/wiki/Wolf%E2%80%93Rayet_star) observed in the first visit of Program 15824. We will download the dataset from the first visit of this program, `LE2701010`, which observes the target at `LP4` using the grating and central wavelength `G160M/1577`.\n", + "\n", + "\n", + "### 1.1 Downloading and Organizing COS Data using `Astroquery`\n", + "\n", + "We will start by downloading our dataset using `astroquery.mast Observations` module; we will go through this *very* briefly, so if you are unfamiliar with downloading MAST data using `astroquery`, please check out our [CoaddTutorial.ipynb](https://github.com/spacetelescope/hst_notebooks/tree/main/notebooks/HASP/CoaddTutorial) notebook.\n", + "\n", + "To run the coadd script, we will need to download the `X1D` files. To perform our wavelength adjustments, we will need to download the `RAWTAG` and `ASN` files since we will be calculating the pixel offset between the `X1D` files and re-running the COS calibration pipeline (`CalCOS`). Let's do so in the cell below:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7c6f2a80", + "metadata": {}, + "outputs": [], + "source": [ + "%%capture\n", + "# Querying MAST for our dataset\n", + "cos_query = Observations.query_criteria(\n", + " obs_id=\"LE2701010\"\n", + ")\n", + "\n", + "# Getting a list of products for this dataset\n", + "cos_products = Observations.get_product_list(\n", + " cos_query\n", + ")\n", + "\n", + "# A list of files that we need to run the coadd script and CalCOS\n", + "necessary_filetypes = [\"X1D\", \"RAWTAG_A\", \"RAWTAG_B\", \"ASN\"]\n", + "\n", + "# Downloading the products\n", + "Observations.download_products(\n", + " cos_products,\n", + " download_dir=str(cos_data),\n", + " productSubGroupDescription=necessary_filetypes\n", + ")\n", + "\n", + "# Moving all files to single directory and deleting mastDownload\n", + "consolidate_files(cos_data)\n", + "\n", + "# Now moving all of the original X1Ds to cos_unshifted\n", + "for file in glob.glob(f\"{cos_data}/*x1d*\"):\n", + " shutil.move(file, f\"{cos_unshifted}/{os.path.basename(file)}\")" + ] + }, + { + "cell_type": "markdown", + "id": "bbc83a6f-2dd1-4052-9825-63a9eb5f8f38", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "\n", + "### 1.2 Running the HASP Script on COS Data\n", + "\n", + "Now that we've downloaded our data, we can run our coadd code. We can run the HASP coadd script, saving our coadded data products to the `./cos_ex/cos_unshifted/coadd/` directory:\n", + "\n", + "**Note: Again, this notebook will not go into specifics about running the script and analyzing the output files. Please check out our [CoaddTutorial.ipynb](https://github.com/spacetelescope/hst_notebooks/tree/main/notebooks/HASP/CoaddTutorial) to learn more about running the script.**" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fc44bbb5-2611-44d1-a60e-aa03e3bb86fe", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "!swrapper -i ./cos_ex/cos_unshifted -o ./cos_ex/cos_unshifted/coadd" + ] + }, + { + "cell_type": "markdown", + "id": "e04fd97a-5938-4c8f-96ef-3ff6aa9676f4", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "We have created two different coadded files by running our script. The file `hst_15824_cos_lmc079-1_g160m_le27_cspec.fits` is our fully coadded spectrum for this program. The other file, `hst_15824_cos_lmc079-1_g160m_le2701_cspec.fits` is the visit-level coadd; if we were to include other visits' `X1D` files then there would be a similar file created for each visit. Each `X1D` was taken at a different grating offset (`FPPOS`), more information about `FPPOS` is in [Section 5.8.2](https://hst-docs.stsci.edu/cosihb/chapter-5-spectroscopy-with-cos/5-8-fixed-pattern-noise) of the COS Instrument Handbook. Let's plot the coadd in the cell below with the `X1D` files:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f984c074-c4ce-4ef4-969c-7b6a7b631878", + "metadata": {}, + "outputs": [], + "source": [ + "# Coadd filename\n", + "coadd_cos = \"hst_15824_cos_lmc079-1_g160m_le27_cspec.fits\"\n", + "\n", + "# List of X1D files\n", + "x1ds = glob.glob(f\"{cos_unshifted}/*x1d*\")\n", + "\n", + "plt.figure(figsize=(15, 6))\n", + "\n", + "colors = [\"blue\", \"red\", \"cyan\", \"orange\"]\n", + "\n", + "# Plotting each X1D file\n", + "for i, file in enumerate(sorted(x1ds)):\n", + " with fits.open(file) as hdul:\n", + "\n", + " wl = hdul[1].data[\"WAVELENGTH\"].ravel()\n", + " flux = hdul[1].data[\"FLUX\"].ravel()\n", + "\n", + " plt.plot(wl, flux,\n", + " lw=2,\n", + " label=f'FPPOS {hdul[0].header[\"FPPOS\"]}',\n", + " color=colors[i],\n", + " alpha=0.6)\n", + "\n", + "# Plotting the coadd\n", + "with fits.open(f\"{cos_unshifted_coadd}/{coadd_cos}\") as hdul:\n", + "\n", + " wl_coadd = hdul[1].data[\"WAVELENGTH\"].ravel()\n", + " flux_coadd = hdul[1].data[\"FLUX\"].ravel()\n", + "\n", + " plt.plot(wl_coadd, flux_coadd,\n", + " lw=2,\n", + " label=\"Coadd\",\n", + " color=\"black\",\n", + " alpha=0.6)\n", + "\n", + "# Formatting and titles\n", + "plt.title(\"Coadd and X1D files - LMC079-1 G160M\")\n", + "plt.xlabel(r'Wavelength [$\\AA$]')\n", + "plt.ylabel(r'Flux [$erg\\ s^{-1}\\ cm^{-2}\\ \\AA^{-1}$]')\n", + "\n", + "plt.ylim(0, 1.5e-13)\n", + "\n", + "plt.legend()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "56883704", + "metadata": {}, + "source": [ + "\n", + "### 1.3 Identifying a Wavelength Shift in COS Data\n", + "\n", + "It is up to the user to determine if their data needs to be shifted, and there are multiple different ways to determine if you need to calculate and apply shifts. It can sometimes be visually apparent that there is a wavelength offset when looking at spectral features as you'll see misalignment, but a more absolute method of determining if your data needs to be shifted is if you compare the coadd's radial velocity to the literature value. In our COS example, we will compare the coadd's radial velocity to the [ULLYSES](https://ullyses.stsci.edu/) High-Level Science Product's (HLSP) radial velocity. We will do a rough, quick calculation to demonstrate that there's an actual shift between the datasets; there are more accurate methods of determining the radial velocity but it is outside the scope of this notebook. The ULLYSES team has created [tutorial notebooks](https://github.com/spacetelescope/ullyses/blob/main/notebooks/ullyses-getrawdata.ipynb) that detail how to download this data via `astroquery`. We will download the HLSP for our target, specifically the file `hlsp_ullyses_hst_cos_lmc079-1_g160m_dr7_cspec.fits`, which is the HLSP for the `G160M` grating." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c6b148f7", + "metadata": {}, + "outputs": [], + "source": [ + "%%capture\n", + "ullyses_query = Observations.query_criteria(\n", + " target_name=\"LMC079-1\",\n", + " provenance_name=\"ULLYSES\",\n", + ")\n", + "\n", + "ullyses_prodlist = Observations.get_product_list(\n", + " ullyses_query\n", + ")\n", + "\n", + "ullyses_prodlist = Observations.filter_products(\n", + " ullyses_prodlist,\n", + " productFilename=\"hlsp_ullyses_hst_cos_lmc079-1_g160m_dr7_cspec.fits\"\n", + ")\n", + "\n", + "Observations.download_products(\n", + " ullyses_prodlist,\n", + " download_dir=str(cos_ex)\n", + ")\n", + "\n", + "# Going to move file to ./cos_ex and get rid of mastDownload + sub-directories\n", + "f = \"hlsp_ullyses_hst_cos_lmc079-1_g160m_dr7_cspec.fits\"\n", + "\n", + "hlsp_path = f\"{cos_ex}/{f}\"\n", + "\n", + "shutil.move(f\"{cos_ex}/mastDownload/HLSP/hlsp_ullyses_hst_cos_lmc079-1_uv/{f}\",\n", + " hlsp_path)\n", + "\n", + "# Now we will remove mastDownload and all the subdirectories\n", + "shutil.rmtree(f\"{cos_ex}/mastDownload\")" + ] + }, + { + "cell_type": "markdown", + "id": "4ecc8ab5", + "metadata": {}, + "source": [ + "Now we can plot our coadd with the ULLYSES HLSP. We will zoom in on two different spectral features: `SiII` at `1526 Å` and `FeII` at `1608 Å`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e472dac5", + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure(figsize=(15, 6))\n", + "\n", + "wavelength_range = ([1607, 1611], [1525, 1530])\n", + "\n", + "# Plotting the coadd\n", + "plt.plot(wl_coadd, flux_coadd,\n", + " lw=2,\n", + " label=\"Coadd\",\n", + " color=\"black\",\n", + " alpha=0.6)\n", + "\n", + "# Plotting the ULLYSES HLSP\n", + "with fits.open(hlsp_path) as hdul:\n", + "\n", + " wl_hlsp = hdul[1].data[\"WAVELENGTH\"].ravel()\n", + " flux_hlsp = hdul[1].data[\"FLUX\"].ravel()\n", + "\n", + "plt.plot(wl_hlsp, flux_hlsp,\n", + " lw=2,\n", + " label=\"ULLYSES HLSP\",\n", + " color=\"blue\",\n", + " alpha=0.6)\n", + "\n", + "plt.title(\"Coadd and ULLYSES HLSP - LMC079-1 G160M\")\n", + "plt.xlabel(r'Wavelength [$\\AA$]')\n", + "plt.ylabel(r'Flux [$erg\\ s^{-1}\\ cm^{-2}\\ \\AA^{-1}$]')\n", + "\n", + "plt.legend()\n", + "\n", + "plt.show()\n", + "\n", + "# Zooming in on two spectral features\n", + "fig, ax = plt.subplots(1, 2, figsize=(15, 6))\n", + "\n", + "for subplot in ax.ravel():\n", + "\n", + " # Plotting the coadd\n", + " subplot.plot(wl_coadd, flux_coadd,\n", + " lw=2,\n", + " color=\"black\",\n", + " alpha=0.6)\n", + "\n", + " # Plotting the HLSP\n", + " subplot.plot(wl_hlsp, flux_hlsp,\n", + " lw=2,\n", + " color=\"blue\",\n", + " alpha=0.6)\n", + "\n", + " subplot.set_xlabel(r'Wavelength [$\\AA$]')\n", + " subplot.set_ylabel(r'Flux [$erg\\ s^{-1}\\ cm^{-2}\\ \\AA^{-1}$]')\n", + "\n", + "ax[0].set_xlim(wavelength_range[0])\n", + "ax[0].set_title(f\"FUVA - ZOOMED {wavelength_range[0]} Å\",\n", + " fontweight=\"bold\")\n", + "\n", + "ax[1].set_xlim(wavelength_range[1])\n", + "ax[1].set_title(f\"FUVB - ZOOMED {wavelength_range[1]} Å\",\n", + " fontweight=\"bold\")\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "c71edeb8", + "metadata": {}, + "source": [ + "We can see visually that there is a slight wavelength shift between the two data products. We will do a very rough estimate of the radial velocity using the `SiII` line at `1526.71 Å`.\n", + "\n", + "We will first fit a Gaussian to our spectral feature to determine the model fit parameters of the absorption feature. Let's define the function below:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ca956262", + "metadata": {}, + "outputs": [], + "source": [ + "def fit_gaussian(wavelength, amplitude, mean, stddev):\n", + " '''\n", + " Fit a Gaussian to the data.\n", + " Gaussian is defined as:\n", + " f(x) = c + amplitude * np.exp(-(((wavelength - mean) / stddev)**2)/2)\n", + " ----------\n", + " Input:\n", + " arr wavelength : Data of the independent variable (wl in our context)\n", + " float amplitude : Amplitude of Gaussian. Represents the peak of the data\n", + " float mean : Mean wl of Gaussian. Represents central wavelength of feature\n", + " float stddev : Standard deviation of Gaussian. Represents function's width\n", + " ----------\n", + " Output:\n", + " The value of flux at each wavelength using the fit.\n", + " '''\n", + " # An estimate of the continuum flux value\n", + " c = 6.5e-14\n", + " return c + amplitude * np.exp(-(((wavelength - mean) / stddev)**2)/2)" + ] + }, + { + "cell_type": "markdown", + "id": "f7fb0acf", + "metadata": {}, + "source": [ + "We will use `scipy.optimize` to help fit our model to the absorption feature by optimizing the parameters of the `fit_gaussian` function. We'll calculate the radial velocity of the HLSP first. This will allow us to determine the exact wavelength corresponding to the peak of the fitted Gaussian, which we'll use to calculate the radial velocity of the HLSP." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6a628010", + "metadata": {}, + "outputs": [], + "source": [ + "# Getting the indicies within the wavelength range\n", + "indicies = np.where((wl_hlsp >= 1525.1) & (wl_hlsp <= 1526.8))\n", + "\n", + "trimmed_flux = flux_hlsp[indicies]\n", + "trimmed_wl = wl_hlsp[indicies]\n", + "\n", + "# We need to make an initial guess for the parameter values\n", + "# The array is in form [PEAK OF GAUSSIAN, CENTRAL WAVELENGTH, WIDTH OF LINE]\n", + "p0 = [.8e-14, 1526.7, 0.0053]\n", + "\n", + "popt_hlsp, _ = optimize.curve_fit(fit_gaussian,\n", + " trimmed_wl,\n", + " trimmed_flux,\n", + " p0)\n", + "\n", + "print(f\"The SiII peak is fitted to be {np.round(popt_hlsp[1], 3)} Å\")\n", + "\n", + "# Plotting the absorption feature of the HLSP with our fit\n", + "plt.figure(figsize=(8, 6))\n", + "\n", + "# Plotting the actual data\n", + "plt.plot(\n", + " wl_hlsp,\n", + " flux_hlsp,\n", + " label=\"HLSP Data\"\n", + ")\n", + "\n", + "# Plotting our fit\n", + "plt.plot(\n", + " wl_hlsp,\n", + " fit_gaussian(wl_hlsp, *popt_hlsp),\n", + " label=\"Model Fit\",\n", + " color=\"red\",\n", + " alpha=0.6\n", + ")\n", + "\n", + "# Adding vertical lines to show fitted \"peak\"\n", + "plt.axvline(popt_hlsp[1],\n", + " label='Model\\'s Peak',\n", + " color=\"black\",\n", + " linestyle=\"dotted\")\n", + "\n", + "plt.xlim(wavelength_range[1])\n", + "plt.legend()\n", + "\n", + "plt.xlabel(r'Wavelength [$\\AA$]')\n", + "plt.ylabel(r'Flux [$erg\\ s^{-1}\\ cm^{-2}\\ \\AA^{-1}$]')\n", + "plt.title(\"HLSP of LMC079-1 with Modeled Fit\",\n", + " fontweight=\"bold\")\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "6673fbc8", + "metadata": {}, + "source": [ + "Now that we've determined the observed wavelength of the `SiII` line, we can calculate the radial velocity. The rest wavelength of this line is `1526.71 Å`. We will use the equation below to calculate the radial velocity:\n", + "\n", + "$ v_r = c \\left( \\frac{\\lambda - \\lambda_0}{\\lambda_0} \\right) $\n", + "\n", + "Where:\n", + "- $ v_r $ is the radial velocity,\n", + "- $ c $ is the speed of light ($ \\approx 2.99 * 10^5 $ km/s),\n", + "- $ \\lambda $ is the observed wavelength of the spectral line,\n", + "- $ \\lambda_0 $ is the rest wavelength of the spectral line.\n", + "\n", + "Let's calculate it for the HLSP below:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d68066ba", + "metadata": {}, + "outputs": [], + "source": [ + "# Speed of light (in km/s)\n", + "c = 2.99e5\n", + "\n", + "# Getting the observed wavelength of our HLSP\n", + "observed_wavelength = popt_hlsp[1]\n", + "\n", + "# Rest wavelength of SiII\n", + "rest_wavelength = 1526.71\n", + "\n", + "v_r = c * ((observed_wavelength - rest_wavelength) / rest_wavelength)\n", + "v_r = np.round(v_r, 2)\n", + "\n", + "print(f\"The radial velocity of LMC079-1 using the HLSP is roughly {v_r} km/s\")" + ] + }, + { + "cell_type": "markdown", + "id": "ea76f7b5", + "metadata": {}, + "source": [ + "Now let's calculate the radial velocity of our coadd using the same process." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6f012032", + "metadata": {}, + "outputs": [], + "source": [ + "# Getting the indicies within the wavelength range\n", + "indicies = np.where((wl_coadd >= 1525.) & (wl_coadd <= 1526.8))\n", + "\n", + "trimmed_flux = flux_coadd[indicies]\n", + "trimmed_wl = wl_coadd[indicies]\n", + "\n", + "# We need to make an initial guess for the parameter values\n", + "# The array is in form [PEAK OF GAUSSIAN, CENTRAL WAVELENGTH, WIDTH OF LINE]\n", + "p0 = [.75e-14, 1526.7, 0.0055]\n", + "\n", + "popt, _ = optimize.curve_fit(fit_gaussian,\n", + " trimmed_wl,\n", + " trimmed_flux,\n", + " p0)\n", + "\n", + "print(f\"The SiII peak is fitted to be {np.round(popt[1], 3)} Å\")\n", + "\n", + "# Plotting the absorption feature of the HLSP with our fit\n", + "plt.figure(figsize=(8, 6))\n", + "\n", + "# Plotting the actual data\n", + "plt.plot(\n", + " wl_coadd,\n", + " flux_coadd,\n", + " label=\"Coadd Data\"\n", + ")\n", + "\n", + "# Plotting our fit\n", + "plt.plot(\n", + " wl_coadd,\n", + " fit_gaussian(wl_coadd, *popt),\n", + " label=\"Model Fit\",\n", + " color=\"red\",\n", + " alpha=0.6\n", + ")\n", + "\n", + "# Adding vertical lines to show fitted \"peak\"\n", + "plt.axvline(popt[1],\n", + " label='Coadd Model\\'s Peak',\n", + " color=\"black\",\n", + " linestyle=\"dotted\")\n", + "\n", + "# Adding vertical line to show observed wl of HLSP\n", + "plt.axvline(popt_hlsp[1],\n", + " label='HLSP Model\\'s Peak',\n", + " color=\"red\",\n", + " linestyle=\"dashed\")\n", + "\n", + "plt.xlim(wavelength_range[1])\n", + "plt.legend()\n", + "\n", + "plt.xlabel(r'Wavelength [$\\AA$]')\n", + "plt.ylabel(r'Flux [$erg\\ s^{-1}\\ cm^{-2}\\ \\AA^{-1}$]')\n", + "\n", + "plt.title(\"Coadd of LMC079-1 with Modeled Fit\",\n", + " fontweight=\"bold\")\n", + "\n", + "plt.show()\n", + "\n", + "# Getting the observed wavelength of our HLSP\n", + "observed_wavelength = popt[1]\n", + "\n", + "v_r = c * ((observed_wavelength - rest_wavelength) / rest_wavelength)\n", + "v_r = np.round(v_r, 2)\n", + "\n", + "print(f\"The radial velocity of LMC079-1 using the coadd is roughly {v_r} km/s\")" + ] + }, + { + "cell_type": "markdown", + "id": "f12aaba8-4b5f-4ae3-a390-95bdebfbcc27", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "We can see that there is a very significant difference between the target's calculated radial velocity using the HLSP (roughly `-17.78 km/s`) and the coadd (roughly `-3.5 km/s`). This is an indicator that there is a wavelength offset between the coadd's constituent `X1D` files. We will need to determine this wavelength offset between the `X1D` files, apply this shift when running `CalCOS`, and redo the coadd. \n", + "\n", + "\n", + "### 1.4 Performing Cross-Correlation on `X1D` COS Data to Calculate Lag\n", + "\n", + "We will determine the shift of our spectrum by cross-correlating one of the `X1D` files to the other files, determining the lag with the greatest cross-correlation coefficient, and then inputting this lag into `CalCOS` to re-generate the `X1D` files. We will need to calculate the shift, in pixels, for both the `FUVA` and `FUVB` segments of each `X1D`. Users will want to cross-correlate along a wavelength range that captures absorption features unique to the target and is free of instrumental effects. Additionally, cross-correlating across a unique feature of the spectrum rather than the entire spectrum greatly improves the cross-correlation accuracy.\n", + "\n", + "We will cross-correlate along the ranges `1525 - 1530 Å` for Segment `FUVB` and `1607 - 1611 Å` for Segment `FUVA`, which capture `SiII` and `FeII` absorption features of this spectrum, respectively. Let's take a closer look at these spectra:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "516c66c0-54a8-490b-a5b8-ee700758ca81", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(3, 1, figsize=(14, 15))\n", + "\n", + "for i, file in enumerate(sorted(x1ds)):\n", + " # Getting the flux and wavelength for our dataset\n", + " data = fits.open(file)[1].data\n", + " wavelength = data[\"WAVELENGTH\"].ravel()\n", + " flux = data[\"FLUX\"].ravel()\n", + "\n", + " fppos = fits.open(file)[0].header[\"FPPOS\"]\n", + "\n", + " # ax[0] is plot with total spectrum, ax[1] is FUVB, ax[2] is FUVA\n", + " for subplot in ax:\n", + " # Plotting the spectrum\n", + " subplot.plot(wavelength, flux,\n", + " lw=2,\n", + " alpha=0.6,\n", + " color=colors[i],\n", + " label=f\"FP-POS {fppos}\")\n", + "\n", + " subplot.set_xlabel(r'Wavelength [$\\AA$]')\n", + " subplot.set_ylabel(r'Flux [$erg\\ s^{-1}\\ cm^{-2}\\ \\AA^{-1}$]')\n", + "\n", + " subplot.set_ylim(0, 1.5e-13)\n", + "\n", + " subplot.legend()\n", + "\n", + "# Adding formatting to the subplots\n", + "ax[0].set_title(\"Flux vs Wavelength, LMC079-1, both segments\",\n", + " fontweight=\"bold\")\n", + "ax[1].set_title(f'Flux vs Wavelength, {wavelength_range[1]} Å, Seg FUVB',\n", + " fontweight=\"bold\")\n", + "ax[2].set_title(f'Flux vs Wavelength, {wavelength_range[0]} Å, Seg FUVA',\n", + " fontweight=\"bold\")\n", + "\n", + "ax[1].set_xlim(wavelength_range[1])\n", + "ax[2].set_xlim(wavelength_range[0])\n", + "\n", + "ax[0].legend()\n", + "\n", + "fig.tight_layout()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "a444bda4-d89f-4ead-b3d6-f7bd568b503c", + "metadata": {}, + "source": [ + "We can see an apparent shift between the `X1D` files, especially between the `FPPOS 1+2` and `FPPOS 3+4` exposures.\n", + "\n", + "We will now determine which `X1D` file will be the reference file. Users have different methods for determining the reference file best suited for your dataset. For this example, we will choose the reference file to be the `X1D` with a model peak closest to the HLSP's peak of the `SiII` absorption feature.\n", + "\n", + "We'll define a function to help us fit the peaks:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b200b2ad", + "metadata": {}, + "outputs": [], + "source": [ + "def fit_model(i, p0, ax, wavelength_range, file):\n", + " '''\n", + " Fit the Gaussian to the feature, return the \"peak\" wavelength\n", + " ----------\n", + " arr i : The wavelength range to trim flux/wl to\n", + " arr p0 : The initial guess for the Gaussian fit\n", + " plt.subplot ax : A specific subplot, used for plotting each X1D\n", + " arr wavelength_range : Wavelength range for the plot\n", + " str file : Path to the X1D file\n", + " '''\n", + " with fits.open(file) as hdul:\n", + " wl = hdul[1].data[\"WAVELENGTH\"].ravel()\n", + " flux = hdul[1].data[\"FLUX\"].ravel()\n", + "\n", + " # Getting the indicies within the wavelength range\n", + " indicies = np.where((wl >= i[0]) & (wl <= i[1]))\n", + "\n", + " trimmed_flux = flux[indicies]\n", + " trimmed_wl = wl[indicies]\n", + "\n", + " popt, _ = optimize.curve_fit(fit_gaussian,\n", + " trimmed_wl,\n", + " trimmed_flux,\n", + " p0)\n", + "\n", + " # Plotting the X1D's data\n", + " ax.plot(\n", + " wl,\n", + " flux\n", + " )\n", + "\n", + " # Plotting the fit\n", + " ax.plot(\n", + " wl,\n", + " fit_gaussian(wl, *popt),\n", + " label=\"Modeled Fit\",\n", + " color=\"purple\"\n", + " )\n", + "\n", + " # Adding vertical line to show peak wl of HLSP\n", + " ax.axvline(popt_hlsp[1],\n", + " label='HLSP Model\\'s Peak',\n", + " color=\"red\",\n", + " linestyle=\"dashed\")\n", + "\n", + " # Adding vertical line to show peak wl of X1D\n", + " ax.axvline(popt[1],\n", + " label='X1D Model\\'s Peak',\n", + " color=\"black\",\n", + " linestyle=\"dotted\")\n", + "\n", + " ax.set_xlim(wavelength_range)\n", + " ax.set_ylim(0, 1.5e-13)\n", + "\n", + " ax.legend()\n", + "\n", + " fppos = fits.getval(file, \"FPPOS\")\n", + " ax.set_title(f'Flux vs Wavelength with Fit, FPPOS {fppos}',\n", + " fontweight=\"bold\")\n", + "\n", + " ax.set_xlabel(r'Wavelength [$\\AA$]')\n", + " ax.set_ylabel(r'Flux [$erg\\ s^{-1}\\ cm^{-2}\\ \\AA^{-1}$]')\n", + "\n", + " peak_diff = np.round(popt[1] - popt_hlsp[1], 3)\n", + "\n", + " ax.text(wavelength_range[1] - 1.35,\n", + " 1.2e-13,\n", + " f'Difference in peak wavelength = {peak_diff} Å',\n", + " ha='center', va='center')\n", + "\n", + " return peak_diff" + ] + }, + { + "cell_type": "markdown", + "id": "420fd779", + "metadata": {}, + "source": [ + "Now we will determine the file to be used as a reference, which will be the file that has the smallest difference in peak wavelength for the `SiII` feature compared to the HLSP peak." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "814d10b1", + "metadata": {}, + "outputs": [], + "source": [ + "fig, axes = plt.subplots(2, 2, figsize=(16, 16))\n", + "\n", + "# A placeholder reference file\n", + "ref_file = x1ds[0]\n", + "\n", + "# A placeholder value for the smallest peak\n", + "smallest_peak = 10000\n", + "\n", + "for file, ax in zip(sorted(x1ds), axes.ravel()):\n", + " fppos = fits.getval(file, \"FPPOS\")\n", + " i = [1525., 1525.8]\n", + "\n", + " if fppos == 1:\n", + " p0 = [-6e-14, 1526.75, 0.13]\n", + "\n", + " elif fppos == 2:\n", + " p0 = [-6e-14, 1526.75, 0.125]\n", + "\n", + " elif fppos == 3:\n", + " p0 = [-5.75e-14, 1526.65, 0.125]\n", + "\n", + " elif fppos == 4:\n", + " p0 = [-5.8e-14, 1526.62, 0.13]\n", + "\n", + " curr_peak = fit_model(i, p0, ax, wavelength_range[1], file)\n", + "\n", + " # Is the current iteration's file closer to the HLSP?\n", + " if np.abs(curr_peak) < smallest_peak:\n", + " smallest_peak = np.abs(curr_peak)\n", + " ref_file = file\n", + "\n", + "plt.show()\n", + "\n", + "fppos = fits.getval(ref_file, \"FPPOS\")\n", + "print(f'The file to be used as reference was taken at FPPOS {fppos}')" + ] + }, + { + "cell_type": "markdown", + "id": "01033e07-3fec-4b9f-a871-0fe63b611d7b", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "As we can see, the file that is closest to the HLSP is the observation taken at `FPPOS 4`.\n", + "\n", + "We will now define a function to determine the lag and cross-correlation coefficient between our reference and shifted spectra. We will use the `correlation_lags` function from `scipy.signal` to compute the cross-correlation coefficient for each lag. The lag is the displacement (in pixels) between the reference and current spectra. The cross-correlation coefficient is in the range `[-1, 1]`, a value of `1` being perfectly correlated, `0` being no correlation, and `-1` being inversely correlated. The shift that we will input to `CalCOS` is going to be the lag at the maximum cross-correlation coefficient value. We will define the function to calculate the lag and cross-correlation coefficient in the cell below (this function is from the [STIS Cross Correlation Notebook](https://github.com/spacetelescope/hst_notebooks/blob/main/notebooks/STIS/cross-correlation/cross-correlation.ipynb))." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f9c18c00-cab6-4553-9ef6-fecdb49903f2", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "# From the STIS cross correlation notebook\n", + "def cross_correlate(ref_flux, shifted_flux):\n", + " assert len(shifted_flux) == len(ref_flux), \"Arrays must be same size\"\n", + "\n", + " # Normalize inputs:\n", + " shifted_flux = shifted_flux - shifted_flux.mean()\n", + " shifted_flux /= shifted_flux.std()\n", + " ref_flux = ref_flux - ref_flux.mean()\n", + " ref_flux /= ref_flux.std()\n", + "\n", + " # centered at the median of len(a)\n", + " lag = correlation_lags(len(shifted_flux), len(ref_flux), mode=\"same\")\n", + " # find the cross-correlation coefficient\n", + " cc = correlate(shifted_flux, ref_flux, mode=\"same\") / float(len(ref_flux))\n", + "\n", + " return lag, cc" + ] + }, + { + "cell_type": "markdown", + "id": "50be84e3-fff3-42bd-a28c-8fadd919dabf", + "metadata": {}, + "source": [ + "Now we will define a function to calculate the shift between our reference file (`FPPOS 4`) and the `FPPOS 1+2+3` exposures. We will plot the correlation coefficient vs lag as well. We will have to get shift values for each segment; for `FUVA` we will correlate to the wavelength range `[1607, 1611]`, and for `FUVB` we will correlate to the wavelength range `[1525, 1530]`. Since our `cross_correlate` function returns only discrete value shifts, we will need to fit a curve to the correlation lag plot; the shift we will apply with `CalCOS` is at the maxima of this curve. \n", + "\n", + "**Note: You will need to experiment with different fits for your dataset.**" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "31223243", + "metadata": {}, + "outputs": [], + "source": [ + "def find_shift_normalized(ref_file, file, wl_range, segment, width, plot):\n", + " '''\n", + " Determine shift amount in pixels between the reference file and the file.\n", + " Can also plot wavelength range and the correlation coeff. vs lag as well.\n", + " ----------\n", + " Input:\n", + " Pathlib path ref_file : Path to reference file's spectrum to shift to.\n", + " Pathlib path file : Path to file of spectrum to apply shifts to\n", + " arr wl_range : Wavelength range to correlate to. Form: [START_WL, END_WL]\n", + " str segment : The segment that is being correlated. Either 'FUVA' or 'FUVB'\n", + " int width : Width of fit, needs to be greater than 2\n", + " bool plot : Do you want to show plots?\n", + " ----------\n", + " Output:\n", + " float shift : Calculated lag [px] between ref and current spectrum\n", + " '''\n", + " # Grabbing the reference spectrum's data\n", + " with fits.open(ref_file) as hdul:\n", + " seg_i = np.where(hdul[1].data[\"SEGMENT\"] == segment)\n", + "\n", + " ref_wl = hdul[1].data[\"WAVELENGTH\"][seg_i]\n", + " ref_flux = hdul[1].data[\"FLUX\"][seg_i]\n", + "\n", + " # Getting the flux data for the current FPPOS file\n", + " with fits.open(file) as hdul:\n", + " seg_i = np.where(hdul[1].data[\"SEGMENT\"] == segment)\n", + " flux = hdul[1].data[\"FLUX\"][seg_i]\n", + "\n", + " range_i = np.where((ref_wl >= wl_range[0]) & (ref_wl <= wl_range[1]))\n", + "\n", + " # Trimming the flux and wavelength ranges to match the user inputted range\n", + " flux = flux[range_i]\n", + " ref_flux = ref_flux[range_i]\n", + "\n", + " wavelength = ref_wl[range_i]\n", + "\n", + " # Calculating lag and correlation coeff. using STIS notebook's function\n", + " lag, cc = cross_correlate(ref_flux, flux)\n", + "\n", + " # Fitting a quadratic to the peak to find pixel shift\n", + " fitter = fitting.LinearLSQFitter()\n", + "\n", + " # Getting the points near the peak\n", + " low, hi = np.argmax(cc) - width//2, np.argmax(cc) + width//2 + 1\n", + " fit = fitter(Polynomial1D(degree=2), x=lag[low:hi], y=cc[low:hi])\n", + "\n", + " x_c = np.arange(np.min(lag), np.max(lag), 0.01)\n", + "\n", + " # Determining the shift, in pixels\n", + " shift1 = -fit.parameters[1] / (2. * fit.parameters[2])\n", + "\n", + " if segment == \"FUVA\":\n", + " shift = shift1 + fits.getval(file, \"SHIFT1A\", ext=1)\n", + "\n", + " else:\n", + " shift = shift1 + fits.getval(file, \"SHIFT1B\", ext=1)\n", + "\n", + " # If not plotting, then just return the shift\n", + " if plot is False:\n", + " return shift\n", + "\n", + " fppos = fits.getval(file, 'FPPOS')\n", + "\n", + " print(f\"Shift for {segment} between ref and FPPOS {fppos} is {shift1} px\")\n", + "\n", + " # Setting up the subplots\n", + " fig, ax = plt.subplots(2, 1, figsize=(16, 8))\n", + "\n", + " ax[0].plot(\n", + " wavelength,\n", + " ref_flux,\n", + " lw=2,\n", + " label=\"Reference File\"\n", + " )\n", + "\n", + " ax[0].plot(\n", + " wavelength,\n", + " flux,\n", + " lw=2,\n", + " label=f'FPPOS {fits.getval(file, \"FPPOS\")}'\n", + " )\n", + "\n", + " ax[1].plot(\n", + " lag,\n", + " cc,\n", + " lw=2,\n", + " color=\"black\"\n", + " )\n", + "\n", + " ax[1].plot(\n", + " x_c,\n", + " fit(x_c),\n", + " lw=2,\n", + " color=\"red\",\n", + " label=\"Fitted quadratic\",\n", + " alpha=0.6\n", + " )\n", + "\n", + " ax[1].axvline(\n", + " shift1,\n", + " linestyle=\"dashed\",\n", + " label=\"Maxima\",\n", + " )\n", + "\n", + " ax[1].legend()\n", + " ax[1].set_ylim(-0.3, 0.8)\n", + " ax[1].set_xlim(shift1-100, shift1+100)\n", + "\n", + " # Adjusting formatting + adding titles\n", + " ax[0].set_xlabel(\"Wavelength\")\n", + " ax[0].set_ylabel(r'Flux [$erg\\ s^{-1}\\ cm^{-2}\\ \\AA^{-1}$]')\n", + " ax[0].set_title(f'Flux vs Wavelength, FPPOS {fppos} & Ref, {segment}',\n", + " weight=\"bold\")\n", + "\n", + " ax[1].set_xlabel(\"Lags (px)\")\n", + " ax[1].set_ylabel(\"Correlation Coeff.\")\n", + " ax[1].set_title(f'CC vs Lags (px), FPPOS {fppos} & Ref, {segment}',\n", + " weight=\"bold\")\n", + "\n", + " ax[0].legend()\n", + "\n", + " plt.tight_layout()\n", + " plt.show()\n", + "\n", + " return" + ] + }, + { + "cell_type": "markdown", + "id": "58efe3da", + "metadata": {}, + "source": [ + "Let's calculate the shift of both segments for the `FPPOS 1` exposure:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b6d7f73a", + "metadata": {}, + "outputs": [], + "source": [ + "fppos1 = f\"{cos_unshifted}/le2701q0q_x1d.fits\"\n", + "\n", + "widths1 = [25, 25]\n", + "\n", + "find_shift_normalized(ref_file,\n", + " fppos1,\n", + " wavelength_range[0],\n", + " \"FUVA\",\n", + " width=widths1[0],\n", + " plot=True)\n", + "\n", + "find_shift_normalized(ref_file,\n", + " fppos1,\n", + " wavelength_range[1],\n", + " \"FUVB\",\n", + " width=widths1[1],\n", + " plot=True)" + ] + }, + { + "cell_type": "markdown", + "id": "51a3d9dd", + "metadata": {}, + "source": [ + "Now we can do the same for our `FPPOS 2` exposure:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2f99cbea", + "metadata": {}, + "outputs": [], + "source": [ + "fppos2 = f\"{cos_unshifted}/le2701q2q_x1d.fits\"\n", + "\n", + "widths2 = [25, 35]\n", + "\n", + "find_shift_normalized(ref_file,\n", + " fppos2,\n", + " wavelength_range[0],\n", + " \"FUVA\",\n", + " width=widths2[0],\n", + " plot=True)\n", + "\n", + "find_shift_normalized(ref_file,\n", + " fppos2,\n", + " wavelength_range[1],\n", + " \"FUVB\",\n", + " width=widths2[1],\n", + " plot=True)" + ] + }, + { + "cell_type": "markdown", + "id": "44623b93", + "metadata": {}, + "source": [ + "And finally, let's calculate the shift for our `FPPOS 3` exposure:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "221939a9", + "metadata": {}, + "outputs": [], + "source": [ + "fppos3 = f\"{cos_unshifted}/le2701q4q_x1d.fits\"\n", + "\n", + "widths3 = [25, 30]\n", + "\n", + "find_shift_normalized(ref_file,\n", + " fppos3,\n", + " wavelength_range[0],\n", + " \"FUVA\",\n", + " width=widths3[0],\n", + " plot=True)\n", + "\n", + "find_shift_normalized(ref_file,\n", + " fppos3,\n", + " wavelength_range[1],\n", + " \"FUVB\",\n", + " width=widths3[1],\n", + " plot=True)" + ] + }, + { + "cell_type": "markdown", + "id": "fcfe7c02-16e6-4498-9561-90d47f53761b", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "\n", + "### 1.5 Applying Shifts to COS data and Re-creating Coadd\n", + "\n", + "Now that we've determined the shift in pixels, we will re-calibrate our data using `CalCOS`. There are two ways that we can apply our shifts: we can either modify the `corrtag` files, or we can put them in a `.txt` file and feed that into `CalCOS`. We will do the latter option in this tutorial. The shift value that we will write in the `.txt` file is the sum of the shift we calculated using the above function and the value of the `SHIFT1A` header keyword of the `X1D` file (`SHIFT1A` for `FUVA`, `SHIFT1B` for `FUVB`).\n", + "\n", + "The `shifts.txt` file must have the following columns:\n", + "\n", + "\\- `ROOTNAME`\n", + "\n", + "\\- `FPOFFSET`\n", + "\n", + "\\- `FLASH/STRIPE/SEGMENT`\n", + "\n", + "\\- New `SHIFT1` value in pixels (dispersion direction shift)\n", + "\n", + "\\- New `SHIFT2` value in pixels (cross-dispersion direction shift)**\n", + "\n", + "_**Note: including the `SHIFT2` value is optional_\n", + "\n", + "If the `ROOTNAME` is not the association's rootname, then you can specify `FPOFFSET` as `any` in the file. This goes for `FLASH/STRIPE/SEGMENT` as well, but we will be specifying the segment in our example because we had to calculate the shifts for each segment. More information about running `CalCOS` when applying the shifts can be found in [Section 5.3.2 of the COS Data Handbook](https://hst-docs.stsci.edu/cosdhb/chapter-5-cos-data-analysis/5-3-working-with-extracted-spectra).\n", + "\n", + "Running the cell below will create a `shifts.txt` file, which will be inputted into `CalCOS` to apply the new shifts to the datasets." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "16e15622-b08e-4b59-9aa1-72ca74845c77", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "with open(f\"{cos_data}/shifts.txt\", \"w\") as shift_file:\n", + " widths = (widths1, widths2, widths3)\n", + " for file, width in zip([fppos1, fppos2, fppos3], widths):\n", + "\n", + " rootname = fits.getval(file, \"ROOTNAME\")\n", + "\n", + " shift_fuva = find_shift_normalized(ref_file,\n", + " file,\n", + " wavelength_range[0],\n", + " \"FUVA\",\n", + " width=width[0],\n", + " plot=False)\n", + "\n", + " shift_file.write(f'{rootname} any any FUVA {shift_fuva}\\n')\n", + "\n", + " shift_fuvb = find_shift_normalized(ref_file,\n", + " file,\n", + " wavelength_range[1],\n", + " \"FUVB\",\n", + " width=width[1],\n", + " plot=False)\n", + "\n", + " shift_file.write(f'{rootname} any any FUVB {shift_fuvb}\\n')" + ] + }, + { + "cell_type": "markdown", + "id": "3b3ecaa2-5284-4367-a249-5f157dbfd584", + "metadata": {}, + "source": [ + "Now we can run `CalCOS` on our data to perform the shift. Check out the COS team's [`CalCOS` notebook](https://github.com/spacetelescope/hst_notebooks/blob/main/notebooks/COS/CalCOS/CalCOS.ipynb) for a tutorial on how to run `CalCOS`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d77b0361-cd56-487c-add8-93509a9d5b6d", + "metadata": {}, + "outputs": [], + "source": [ + "%%capture\n", + "calcos.calcos(\n", + " f\"{cos_data}/le2701010_asn.fits\",\n", + " outdir=str(cos_shifted),\n", + " shift_file=f\"{cos_data}/shifts.txt\")" + ] + }, + { + "cell_type": "markdown", + "id": "202dda8b", + "metadata": {}, + "source": [ + "Let's compare the pre-shifted and shifted wavelength ranges of our observations for both segments:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "131a5c60", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(2, 2, figsize=(16, 8))\n", + "\n", + "# Plotting the unshifted data on the top row\n", + "for i, file in enumerate(sorted(glob.glob(f\"{cos_unshifted}/*x1d.fits*\"))):\n", + "\n", + " with fits.open(file) as hdul:\n", + " wavelength = hdul[1].data[\"WAVELENGTH\"].flatten()\n", + " flux = hdul[1].data[\"FLUX\"].flatten()\n", + "\n", + " for j in range(0, 2):\n", + " ax[0][j].plot(\n", + " wavelength,\n", + " flux,\n", + " lw=2,\n", + " alpha=0.6,\n", + " label=f'FPPOS {fits.getval(file, \"FPPOS\")}',\n", + " color=colors[i]\n", + " )\n", + "\n", + "# Now plotting the shifted data on the bottom row\n", + "for i, file in enumerate(sorted(glob.glob(f\"{cos_shifted}/*x1d.fits*\"))):\n", + " with fits.open(file) as hdul:\n", + "\n", + " wavelength = hdul[1].data[\"WAVELENGTH\"].flatten()\n", + " flux = hdul[1].data[\"FLUX\"].flatten()\n", + "\n", + " for j in range(0, 2):\n", + " ax[1][j].plot(\n", + " wavelength,\n", + " flux,\n", + " lw=2,\n", + " alpha=0.6,\n", + " label=f'FPPOS {fits.getval(file, \"FPPOS\")}',\n", + " color=colors[i]\n", + " )\n", + "\n", + "# Adding titles\n", + "fig.suptitle(\"Flux vs Wavelength, Pre-shifted (top) and Shifted (bottom)\",\n", + " weight=\"bold\",\n", + " fontsize=14)\n", + "\n", + "ax[0][0].set_title(f\"Pre-shifted, FUVA, {wavelength_range[0]}\")\n", + "ax[0][1].set_title(f\"Pre-shifted, FUVB, {wavelength_range[1]}\")\n", + "\n", + "ax[1][0].set_title(f\"Shifted, FUVA, {wavelength_range[0]}\")\n", + "ax[1][1].set_title(f\"Shifted, FUVB, {wavelength_range[1]}\")\n", + "\n", + "# Setting xlimits\n", + "ax[0][0].set_xlim(wavelength_range[0])\n", + "ax[0][1].set_xlim(wavelength_range[1][0], wavelength_range[1][1])\n", + "\n", + "ax[1][0].set_xlim(wavelength_range[0])\n", + "ax[1][1].set_xlim(wavelength_range[1][0], wavelength_range[1][1])\n", + "\n", + "# Adding a legend and setting ylim\n", + "for subplot in ax.flatten():\n", + " subplot.legend()\n", + " subplot.set_ylim(0, 1.5e-13)\n", + "\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "5e6a2733-b6ff-491f-bfea-554a44f13066", + "metadata": {}, + "source": [ + "We can visually see a clear improvement in the spectral features when we applied our shifts to `CalCOS` as the two absorption features have lined up for all `FPPOS`. \n", + "\n", + "\n", + "### 1.6 Comparing Non-Shifted and Shifted COS Coadds\n", + "\n", + "Now that we've successfully shifted our data, let's re-create our coadd. We will calculate the radial velocity of the new coadd as well and compare it with our values from the earlier section." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dfb56c66-0988-4d0b-b45c-b0b0081491c1", + "metadata": {}, + "outputs": [], + "source": [ + "!swrapper -i ./cos_ex/cos_shifted -o ./cos_ex/cos_shifted/coadd" + ] + }, + { + "cell_type": "markdown", + "id": "916a5850-e3d8-4128-988a-0bb572a5c66e", + "metadata": {}, + "source": [ + "Let's see the flux and SNR of our coadd now that we've recalibrated:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a5233a65", + "metadata": {}, + "outputs": [], + "source": [ + "coadd = \"hst_15824_cos_lmc079-1_g160m_le27_cspec.fits\"\n", + "\n", + "# The paths to both of our coadds\n", + "preshift = f\"{cos_unshifted_coadd}/{coadd}\"\n", + "shifted = f\"{cos_shifted_coadd}/{coadd}\"\n", + "\n", + "# Getting the data for our pre- and post-shifted coadd\n", + "with fits.open(preshift) as hdul:\n", + " pre_wl = hdul[1].data[\"WAVELENGTH\"].ravel()\n", + " pre_flux = hdul[1].data[\"FLUX\"].ravel()\n", + " pre_snr = hdul[1].data[\"SNR\"].ravel()\n", + "\n", + "with fits.open(shifted) as hdul:\n", + " shifted_wl = hdul[1].data[\"WAVELENGTH\"].ravel()\n", + " shifted_flux = hdul[1].data[\"FLUX\"].ravel()\n", + " shifted_snr = hdul[1].data[\"SNR\"].ravel()\n", + "\n", + "fig, ax = plt.subplots(2, 2, figsize=(14, 8))\n", + "\n", + "# Interpolating flux and SNR arrays since difference shapes from pre-shifted\n", + "shifted_flux = np.interp(pre_wl, shifted_wl, shifted_flux)\n", + "shifted_snr = np.interp(pre_wl, shifted_wl, shifted_snr)\n", + "\n", + "# Plotting flux vs wavelength\n", + "ax[0][0].plot(\n", + " pre_wl,\n", + " pre_flux,\n", + " label=\"Unshifted\",\n", + " lw=2,\n", + " alpha=0.6,\n", + " color=\"red\"\n", + ")\n", + "\n", + "ax[0][0].plot(\n", + " pre_wl,\n", + " shifted_flux,\n", + " label=\"Shifted\",\n", + " lw=2,\n", + " alpha=0.6,\n", + " color=\"blue\"\n", + ")\n", + "\n", + "# Plotting SNR vs wavelength\n", + "ax[0][1].plot(\n", + " pre_wl,\n", + " pre_snr,\n", + " label=\"Unshifted\",\n", + " lw=2,\n", + " alpha=0.6,\n", + " color=\"red\"\n", + ")\n", + "\n", + "ax[0][1].plot(\n", + " pre_wl,\n", + " shifted_snr,\n", + " label=\"Shifted\",\n", + " lw=2,\n", + " alpha=0.6,\n", + " color=\"blue\"\n", + ")\n", + "\n", + "# Plotting flux residuals\n", + "ax[1][0].plot(\n", + " pre_wl,\n", + " shifted_flux-pre_flux,\n", + " lw=2,\n", + " color=\"black\",\n", + " alpha=0.8\n", + ")\n", + "\n", + "# Plotting snr residuals\n", + "ax[1][1].plot(\n", + " pre_wl,\n", + " shifted_snr-pre_snr,\n", + " lw=2,\n", + " color=\"black\",\n", + " alpha=0.8\n", + ")\n", + "\n", + "for subplot in ax.ravel():\n", + " subplot.set_xlabel(r'Wavelength [$\\AA$]')\n", + " subplot.legend()\n", + "\n", + "ax[0][0].set_ylabel(r'Flux [$erg\\ s^{-1}\\ cm^{-2}\\ \\AA^{-1}$]')\n", + "ax[1][0].set_ylabel(r'Flux [$erg\\ s^{-1}\\ cm^{-2}\\ \\AA^{-1}$]')\n", + "\n", + "ax[0][1].set_ylabel(r'SNR')\n", + "ax[1][1].set_ylabel(r'SNR')\n", + "\n", + "ax[0][0].set_title(\"Flux vs Wavelength, Unshifted and Shifted\",\n", + " fontweight=\"bold\")\n", + "ax[1][0].set_title(\"Flux Residuals\",\n", + " fontweight=\"bold\")\n", + "\n", + "ax[0][1].set_title(\"SNR vs Wavelength, Unshifted and Shifted\",\n", + " fontweight=\"bold\")\n", + "ax[1][1].set_title(\"SNR Residuals\",\n", + " fontweight=\"bold\")\n", + "\n", + "plt.tight_layout()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "4f9ed39f", + "metadata": {}, + "source": [ + "While you won't necessarily see a big improvement in the SNR with a shift, we will see that our coadd has a better wavelength solution by calculating the radial velocity below:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8754ae22", + "metadata": {}, + "outputs": [], + "source": [ + "new_coadd = f\"{cos_shifted_coadd}/hst_15824_cos_lmc079-1_g160m_le27_cspec.fits\"\n", + "\n", + "with fits.open(new_coadd) as hdul:\n", + " data_new_coadd = hdul[1].data\n", + " wl_new_coadd = data_new_coadd[\"WAVELENGTH\"].ravel()\n", + " flux_new_coadd = data_new_coadd[\"FLUX\"].ravel()\n", + "\n", + "# Getting the indicies within the wavelength range\n", + "indicies = np.where((wl_new_coadd >= 1525.) & (wl_new_coadd <= 1525.75))\n", + "\n", + "trimmed_flux = flux_new_coadd[indicies]\n", + "trimmed_wl = wl_new_coadd[indicies]\n", + "\n", + "# We need to make an initial guess for the parameter values\n", + "# The array is in form [PEAK OF GAUSSIAN, CENTRAL WAVELENGTH, WIDTH OF LINE]\n", + "p0 = [-0.58e-13, 1526.63, 0.12]\n", + "\n", + "popt, _ = optimize.curve_fit(fit_gaussian,\n", + " trimmed_wl,\n", + " trimmed_flux,\n", + " p0)\n", + "\n", + "print(f\"The SiII peak is fitted to be {np.round(popt[1], 3)} Å\")\n", + "\n", + "# Plotting the absorption feature of the HLSP with our fit\n", + "plt.figure(figsize=(8, 6))\n", + "\n", + "# Plotting the actual data\n", + "plt.plot(\n", + " wl_new_coadd,\n", + " flux_new_coadd,\n", + " label=\"Coadd Data\"\n", + ")\n", + "\n", + "# Plotting our fit\n", + "plt.plot(\n", + " wl_new_coadd,\n", + " fit_gaussian(wl_new_coadd, *popt),\n", + " label=\"Model Fit\",\n", + " color=\"red\",\n", + " alpha=0.6\n", + ")\n", + "\n", + "# Adding vertical lines to show fitted \"peak\"\n", + "plt.axvline(popt[1],\n", + " label='Coadd Model\\'s Peak',\n", + " color=\"black\",\n", + " linestyle=\"dotted\")\n", + "\n", + "# Adding vertical line to show observed wl of HLSP\n", + "plt.axvline(popt_hlsp[1],\n", + " label='HLSP Model\\'s Peak',\n", + " color=\"red\",\n", + " linestyle=\"dashed\")\n", + "\n", + "plt.xlim(wavelength_range[1])\n", + "plt.legend()\n", + "\n", + "plt.xlabel(r'Wavelength [$\\AA$]')\n", + "plt.ylabel(r'Flux [$erg\\ s^{-1}\\ cm^{-2}\\ \\AA^{-1}$]')\n", + "\n", + "plt.title(\"Coadd of LMC079-1 with Modeled Fit\",\n", + " fontweight=\"bold\")\n", + "\n", + "plt.show()\n", + "\n", + "# Getting the observed wavelength of our HLSP\n", + "observed_wavelength = popt[1]\n", + "\n", + "v_r = c * ((observed_wavelength - rest_wavelength) / rest_wavelength)\n", + "v_r = np.round(v_r, 2)\n", + "\n", + "print(f\"The radial velocity of LMC079-1 using the coadd is roughly {v_r} km/s\")" + ] + }, + { + "cell_type": "markdown", + "id": "03f396c0-15b2-4f13-bb9b-80016f65c047", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "The radial velocity of the coadd created with the shifted `X1D` files (`-15.67 km/s`) is a lot closer to the HLSP's radial velocity (`-17.78 km/s`) compared to the radial velocity of the pre-shifted coadd (`-3.5 km/s`). We can now move on to the STIS example.\n", + "\n", + "\n", + "## 2. Performing Wavelength Adjustments on STIS Data\n", + "\n", + "At the beginning of a science exposure, two Pt-Cr/Ne line lamps are flashed to obtain wavelength calibration exposures and assign a zero-point for the wavelength scale. These wavecals are taken at the beginning of virtually every exposure and are automatically taken at each `2300sec` interval thereafter. Similar to COS, shifts between exposures can occur due to imperfect target acquisitions, repositioning of the grating, and thermal flexure. In these instances, users will need to manually shift their spectra to ensure the correct wavelength and flux calibrations are being performed. Echelle data and low-resolution spectra are suspectable to rapid sensitivity changes with wavelength, which makes accurate wavelength solutions even more important. Information about STIS wavelength accuracies can be found in [Chapter 4.2](https://hst-docs.stsci.edu/stisdhb/chapter-4-stis-error-sources/4-2-summary-of-accuracies) of the STIS Instrument Handbook.\n", + "\n", + "We will be shifting data taken with the `STIS/CCD`, using the `52X2` aperture. There is a small pixel shift between data taken at the nominal position and data taken at the `E1` pseudo-aperture position for this STIS mode (mentioned in [Joyce 2018](https://ui.adsabs.harvard.edu/abs/2018MNRAS.481.2361J/abstract)). We will download our data using `astroquery`, run the HASP script, cross-correlate the spectra in the dataset to the brightest observation to calculate the wavelength shift, and then apply this shift using `CalSTIS`.\n", + "\n", + "**Note: many of the cross-correlation steps outlined in this notebook are detailed in the [STIS Cross-Correlation Jupyter Notebook](https://github.com/spacetelescope/hst_notebooks/blob/main/notebooks/STIS/cross-correlation/cross-correlation.ipynb)**\n", + "\n", + "\n", + "### 2.1 Downloading and Organizing STIS Data using `Astroquery`\n", + "\n", + "We will be downloading `STIS/CCD G430L` data for the object, `AGK+81D266`, a spectrophotometric target taken in Program 14423. Let's create a few directories for our investigation: " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e1d547ee-5a6d-4d29-a6f7-0e158107dfea", + "metadata": {}, + "outputs": [], + "source": [ + "# To hold all of our STIS example's output and data\n", + "stis_ex = Path(\"./stis_ex\")\n", + "\n", + "# To contain the raw files necessary to run CalSTIS\n", + "stis_data = Path(\"./stis_ex/stis_data/\")\n", + "\n", + "# To hold the unshifted data products and coadds\n", + "stis_unshifted = Path(\"./stis_ex/stis_unshifted/\")\n", + "stis_unshifted_coadd = Path(\"./stis_ex/stis_unshifted/coadd/\")\n", + "\n", + "# To contain the shifted data products and coadds\n", + "stis_shifted = Path(\"./stis_ex/stis_shifted/\")\n", + "stis_shifted_coadd = Path(\"./stis_ex/stis_shifted/coadd/\")\n", + "\n", + "# If the directory doesn't exist then create it\n", + "stis_ex.mkdir(exist_ok=True)\n", + "\n", + "stis_data.mkdir(exist_ok=True)\n", + "\n", + "stis_unshifted.mkdir(exist_ok=True)\n", + "stis_unshifted_coadd.mkdir(exist_ok=True)\n", + "\n", + "stis_shifted.mkdir(exist_ok=True)\n", + "stis_shifted_coadd.mkdir(exist_ok=True)" + ] + }, + { + "cell_type": "markdown", + "id": "eecde454-8f9a-41e3-bad9-ea6d70756cd4", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "We can use `astroquery` to query our datasets, get a list of associated products, and download the data to our `./stis_data` directory. This is similar to how we downloaded the COS data in the previous example." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "651fa78e", + "metadata": {}, + "outputs": [], + "source": [ + "%%capture\n", + "# Querying our datasets\n", + "stis_query = Observations.query_criteria(\n", + " proposal_id=14423,\n", + " instrument_name=\"STIS/CCD\",\n", + " filters=\"G430L\",\n", + " target_name=\"AGK+81D266\"\n", + ")\n", + "\n", + "# Getting the list of products for our datasets\n", + "stis_prodlist = Observations.get_product_list(\n", + " stis_query\n", + ")\n", + "\n", + "# Some of the files aren't coadded by default because POSTARG != 0\n", + "# We won't download these to save memory space\n", + "bad_obsids = [\n", + " \"od1am1090\",\n", + " \"od1am10a0\",\n", + " \"od1am10c0\",\n", + " \"od1am10d0\",\n", + " \"od1am10f0\"\n", + "]\n", + "\n", + "stis_prodlist = stis_prodlist[~np.isin(stis_prodlist[\"obs_id\"], bad_obsids)]\n", + "\n", + "# Filtering the product list to only get the required files\n", + "stis_prodlist = Observations.filter_products(\n", + " stis_prodlist,\n", + " productSubGroupDescription=[\"SX1\", \"RAW\", \"ASN\", \"FLT\"]\n", + ")\n", + "\n", + "# Downloading them to ./stis_data\n", + "Observations.download_products(\n", + " stis_prodlist,\n", + " download_dir=str(stis_data),\n", + ")\n", + "\n", + "# Moving all sub-folders to stis_data and deleting mastDownload folder\n", + "consolidate_files(stis_data)\n", + "\n", + "# Moving all sx1 to stis_unshifted\n", + "for file in glob.glob(f\"{stis_data}/*sx1.fits*\"):\n", + " # New path for file\n", + " newpath = f\"{stis_unshifted}/{fits.getval(file, 'FILENAME')}\"\n", + "\n", + " # Moving the actual file\n", + " shutil.move(file, newpath)" + ] + }, + { + "cell_type": "markdown", + "id": "66ef8215-5ace-4fc8-9757-6cb276353975", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "\n", + "### 2.2 Running the HASP Script on STIS Data\n", + "\n", + "Now that we have our data downloaded, we can run the coadd script and create abutted products. Check out our [CoaddTutorial.ipynb](https://github.com/spacetelescope/hst_notebooks/tree/main/notebooks/HASP/CoaddTutorial) notebook for a more detailed look as to how to run the script and change the flag parameters." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f76f03e9-16db-4c9b-bcec-d317454a5402", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "!swrapper -i ./stis_ex/stis_unshifted -o ./stis_ex/stis_unshifted/coadd" + ] + }, + { + "cell_type": "markdown", + "id": "fe615f38-c269-48ab-bdb1-d3f89f0c7f32", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "We will plot the coadd over our current `SX1` files; the aperture used is also listed in the legend. Observations that were taken at the `E1` position have the value `52X2E1` for the `PROPAPER` header keyword (and are dotted on the plot), whereas observations taken at the nominal position have the value `52X2`. Notice the area around `4860Å` in the second plot shows a visible shift between datasets taken at the nominal and `E1` positions. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "650c8622-4d35-4e44-8e50-58336550796e", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "coadd_stis = \"hst_14423_stis_agkp81d266_g430l_od1a_cspec.fits\"\n", + "\n", + "# A list of the SX1 files for our data\n", + "sx1s = glob.glob(f\"{stis_unshifted}/*sx1*\")\n", + "\n", + "# Creating a figure with two subplots\n", + "fig, ax = plt.subplots(2, 1, figsize=(14, 12))\n", + "\n", + "for curr_plot in ax:\n", + " # Plotting each SX1 file\n", + " for file in sx1s:\n", + " with fits.open(file) as hdul:\n", + " wl = hdul[1].data[\"WAVELENGTH\"].ravel()\n", + " flux = hdul[1].data[\"FLUX\"].ravel()\n", + "\n", + " aperture = hdul[0].header['PROPAPER']\n", + "\n", + " # Differentiating between nominal and E1 position via linestyle\n", + " if aperture == \"52X2\":\n", + " linestyle = \"solid\"\n", + " label = f\"{hdul[0].header['ROOTNAME']} - nom.\"\n", + " else:\n", + " linestyle = \"dotted\"\n", + " label = f\"{hdul[0].header['ROOTNAME']} - E1\"\n", + "\n", + " # Plotting the actual data\n", + " curr_plot.plot(wl, flux,\n", + " label=label,\n", + " alpha=0.75,\n", + " lw=2,\n", + " linestyle=linestyle)\n", + "\n", + " # Plotting the coadd in black\n", + " with fits.open(f\"{stis_unshifted_coadd}/{coadd_stis}\") as hdul:\n", + " wl = hdul[1].data[\"WAVELENGTH\"].ravel()\n", + " flux = hdul[1].data[\"FLUX\"].ravel()\n", + "\n", + " curr_plot.plot(wl, flux,\n", + " label=\"Coadd\",\n", + " alpha=0.75,\n", + " color=\"black\",\n", + " lw=2)\n", + "\n", + "# Adding labels, formatting, xlim\n", + "plt.suptitle(\"SX1 Files against Coadd - AGK+81D266\")\n", + "plt.xlabel(r'Wavelength [$\\AA$]')\n", + "ax[0].set_ylabel(r'Flux [$erg\\ s^{-1}\\ cm^{-2}\\ \\AA^{-1}$]')\n", + "ax[1].set_ylabel(r'Flux [$erg\\ s^{-1}\\ cm^{-2}\\ \\AA^{-1}$]')\n", + "\n", + "ax[1].set_xlim(4825, 4900)\n", + "ax[1].set_ylim(0.4e-13, 1.5e-13)\n", + "\n", + "ax[0].legend()\n", + "\n", + "plt.tight_layout()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "0c7bef35-5a11-4724-8f2b-77fefdf80cb7", + "metadata": {}, + "source": [ + "Since we see a consistent shift in the data between the observations that used the `E1` aperture and the nominal, we need to determine the actual shift in pixels and re-run `CalSTIS` on our data. We will not be calculating the radial velocity and comparing it to the literature value to reduce notebook length, and because the shift is known to be relative to the nominal position.\n", + "\n", + "\n", + "### 2.3 Performing Cross-Correlation on STIS data and Calculating Lag\n", + "\n", + "There are multiple ways to determine the best observation to set as the reference, and that is up to the user to do so. Since we are shifting all the `E1` observations to the nominal positions, we will use the brightest nominal observation as our reference. We are using the brightest because we are making the assumption that it is the observation most centered in the aperture, but you will need to use other methods when determining the file to reference to. Let's determine the brightest nominal file below:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cbb94765-72b7-4469-bb93-9c663b46abdf", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "# A variable for the brightest filename, arbitrarily\n", + "ref_file = sx1s[0]\n", + "\n", + "for file in sx1s:\n", + " # Getting flux data for current interation's file + current brightest file\n", + " curr_flux = fits.open(file)[1].data[\"FLUX\"].ravel()\n", + " ref_flux = fits.open(ref_file)[1].data[\"FLUX\"].ravel()\n", + "\n", + " aperture = fits.getval(file, \"PROPAPER\")\n", + "\n", + " # If current file's mean flux > ref file, replace ref file\n", + " if ((np.mean(curr_flux) > np.mean(ref_flux)) & (aperture == \"52X2\")):\n", + " ref_file = file\n", + "\n", + " else:\n", + " continue\n", + "\n", + "# Printing out our file\n", + "ref_head = fits.open(ref_file)[0].header\n", + "print(f\"The brightest file is {ref_head['FILENAME']}.\")" + ] + }, + { + "cell_type": "markdown", + "id": "f410862f-7e5a-4457-b76f-4beec661040e", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "We've now determined that `od1al2030_sx1.fits` is the brightest spectra of our data. We will go through each `SX1` file taken at the `E1` position in our data and cross-correlate it to our `od1al2030_sx1.fits` file. We will calculate the shifts for a single file first (to better illustrate each step), then define a function to perform these same steps for the rest of the files.\n", + "\n", + "We will be cross-correlating similarly to how we did our COS example in the previous section, but this section will mainly follow the steps in the official [STIS Cross-Correlation Notebook](https://github.com/spacetelescope/hst_notebooks/blob/main/notebooks/STIS/cross-correlation/cross-correlation.ipynb).\n", + "\n", + "We will need to calculate the mean plate scale, i.e. the dispersion per pixel, for our reference dataset. This will be used to determine the shift in pixels later on in the tutorial." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9a3aad93-4f90-4f8f-af8a-8ba58f42268d", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "# Getting a list of files taken at E1 position\n", + "e1_pos = [file for file in sx1s if fits.getval(file, \"PROPAPER\") != \"52X2\"]\n", + "\n", + "# Example file\n", + "stis_file = e1_pos[0]\n", + "\n", + "# Getting the data for the reference file\n", + "with fits.open(ref_file) as hdul:\n", + " ref_wl = hdul[1].data[\"WAVELENGTH\"].ravel()\n", + " ref_flux = hdul[1].data[\"FLUX\"].ravel()\n", + "\n", + "with fits.open(stis_file) as hdul:\n", + " ex_wl = hdul[1].data[\"WAVELENGTH\"].ravel()\n", + " ex_flux = hdul[1].data[\"FLUX\"].ravel()\n", + "\n", + "# Interpolating the example spectra\n", + "ex_flux = np.interp(ref_wl, ex_wl, ex_flux)\n", + "\n", + "# Getting the dispersion per pixel\n", + "mean_plate_scale = np.mean(ref_wl[1:] - ref_wl[:-1])\n", + "print(f\"The dispersion per pixel is {np.round(mean_plate_scale, 2)} Å/pixel\")" + ] + }, + { + "cell_type": "markdown", + "id": "0fb025ec-ad26-4686-9520-65f9acaa9eed", + "metadata": {}, + "source": [ + "Now we will use our cross-correlation function defined in the COS section to determine the lag and cross-correlation coefficient. We can determine the shift by getting the lag value with the highest corresponding cross-correlation coefficient. Our lags are discrete pixel integer values, so we will need to fit a quadratic curve using `fitting` and `Polynomial1D` from the `astropy.modeling` module. The lag at the maxima of our fit is the pixel shift that we will use when re-calibrating with `CalSTIS`.\n", + "\n", + "**Note: you will need to experiment with different fits when determining your lag; there is no single function that will work in all scenarios**" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fa56442e", + "metadata": {}, + "outputs": [], + "source": [ + "# Calculating the lag and cc coeff using the earlier defined function\n", + "lag, cc = cross_correlate(ref_flux, ex_flux)\n", + "\n", + "# Fitting a quadratic to the peak to find pixel shift\n", + "fitter = fitting.LinearLSQFitter()\n", + "\n", + "# Getting the 8 points near the peak\n", + "width = 8\n", + "\n", + "low, hi = np.argmax(cc) - width//2, np.argmax(cc) + width//2 + 1\n", + "fit = fitter(Polynomial1D(degree=2), x=lag[low:hi], y=cc[low:hi])\n", + "\n", + "x_c = np.arange(np.min(lag), np.max(lag), 0.01)\n", + "\n", + "# Determining the shift, in pixels\n", + "shift1 = -fit.parameters[1] / (2. * fit.parameters[2])\n", + "\n", + "print(f\"The shift, in pixels, is {np.round(shift1, 2)} pixels\")\n", + "print(f\"The shift, in angstroms, is {np.round(shift1*mean_plate_scale, 2)}Å\")\n", + "\n", + "plt.plot(\n", + " lag,\n", + " cc,\n", + " lw=2,\n", + " color=\"black\"\n", + ")\n", + "\n", + "plt.plot(\n", + " x_c,\n", + " fit(x_c),\n", + " lw=2,\n", + " label=\"Fitted Quadratic\",\n", + " alpha=0.6\n", + ")\n", + "\n", + "plt.axvline(\n", + " shift1,\n", + " color=\"red\",\n", + " label=\"Maxima of fit\"\n", + ")\n", + "\n", + "# Adding title\n", + "plt.title(\"Lag (px) vs Cross Correlation Coeff, with fitted curve and ZOOMED\",\n", + " weight=\"bold\")\n", + "\n", + "plt.xlabel(\"Lag (px)\")\n", + "plt.ylabel(\"Cross Correlation Coeff.\")\n", + "\n", + "plt.ylim(0.9, 1.05)\n", + "plt.xlim(shift1 - 25, shift1 + 25)\n", + "\n", + "plt.legend()\n", + "\n", + "plt.tight_layout()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "ebbc8604-4a1c-48ca-84a7-16985e90a9c6", + "metadata": {}, + "source": [ + "\n", + "### 2.4 Applying Shifts to Raw STIS data and Re-creating Coadd\n", + "Now that we've calculated the shift, we need to alter the header keywords of our data and then re-calibrate the spectrum using `CalSTIS`. We will extract the `SHIFTA1` and `SHIFTA2` keywords from the corresponding `FLT` file, add our calculated shift value to the value of this keyword, then set the `WAVECORR` calibration switch in the `RAW` file to `OMIT` and re-run `CalSTIS`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bfa759b2", + "metadata": {}, + "outputs": [], + "source": [ + "# Path to the FLT for our example file\n", + "flt_file = f\"{stis_data}/{fits.getval(stis_file, 'ROOTNAME')}_flt.fits\"\n", + "\n", + "# Extracting the shift from the FLT file\n", + "SHIFTA1 = fits.getval(flt_file, \"SHIFTA1\", ext=1)\n", + "SHIFTA2 = fits.getval(flt_file, \"SHIFTA2\", ext=1)\n", + "\n", + "# Update SHIFTA1 (spectral direction) to take into account the calculated shift\n", + "SHIFTA1 += shift1\n", + "\n", + "# Path to the raw file\n", + "raw_file = f\"{stis_data}/{fits.getval(stis_file, 'ROOTNAME')}_raw.fits\"\n", + "\n", + "# Updating the header to reflect updated shift + FLT SHIFTA2 value\n", + "fits.setval(raw_file, keyword=\"SHIFTA1\", value=SHIFTA1, ext=1)\n", + "fits.setval(raw_file, keyword=\"SHIFTA2\", value=SHIFTA2, ext=1)\n", + "\n", + "# Set WAVECORR to OMIT since we're applying the shift\n", + "fits.setval(raw_file, \"WAVECORR\", value=\"OMIT\", ext=0)" + ] + }, + { + "cell_type": "markdown", + "id": "c81a587c-4091-4b20-8385-42e8761d8bb7", + "metadata": {}, + "source": [ + "Now we can run `CalSTIS`:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "333a1f70-9f42-41a4-ba97-3abe549253ec", + "metadata": {}, + "outputs": [], + "source": [ + "%%capture\n", + "stistools.calstis.calstis(raw_file,\n", + " outroot=f\"{str(stis_shifted)}/\")" + ] + }, + { + "cell_type": "markdown", + "id": "622ef089", + "metadata": {}, + "source": [ + "We will now define a function to calculate and perform shifts on the rest of the `E1` position datasets. We have three remaining files to shift, so we will define a function to calculate the fit, plot, and run `CalSTIS`. \n", + "\n", + "**Note that you will need to experiment with different fits/functions for your dataset.**" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3735400c", + "metadata": {}, + "outputs": [], + "source": [ + "def shift_stis(ref_file, stis_file, width, plot, run_calstis=False):\n", + " '''\n", + " Determing the shift, in pixels, between the file and the reference.\n", + " Update the RAW file with the new shift and run CalSTIS.\n", + " ----------\n", + " Input:\n", + " pathlib.Path or str ref_file : Path to the reference file\n", + " pathlib.Path or str stis_file : Path to the file that needs to be shifted\n", + " int width : Number points to fit at the peak of the function, must be > 2\n", + " bool plot : Do you want to plot the cc function and the fit?\n", + " bool run_calstis : True if you want to run CalSTIS on the shifted dataset\n", + " '''\n", + " # Getting the data for the reference file\n", + " with fits.open(ref_file) as hdul:\n", + " ref_wl = hdul[1].data[\"WAVELENGTH\"].ravel()\n", + " ref_flux = hdul[1].data[\"FLUX\"].ravel()\n", + "\n", + " with fits.open(stis_file) as hdul:\n", + " ex_wl = hdul[1].data[\"WAVELENGTH\"].ravel()\n", + " ex_flux = hdul[1].data[\"FLUX\"].ravel()\n", + "\n", + " # Interpolating the example spectra\n", + " ex_flux = np.interp(ref_wl, ex_wl, ex_flux)\n", + "\n", + " # Getting the dispersion per pixel\n", + " mean_plate_scale = np.mean(ref_wl[1:] - ref_wl[:-1])\n", + "\n", + " print(f\"Dispersion per pixel is {np.round(mean_plate_scale, 2)} Å/pixel\")\n", + "\n", + " # Calculating the lag and cc coeff using the earlier defined function\n", + " lag, cc = cross_correlate(ref_flux, ex_flux)\n", + "\n", + " # Fitting a quadratic to the peak to find pixel shift\n", + " fitter = fitting.LinearLSQFitter()\n", + "\n", + " low, hi = np.argmax(cc) - width//2, np.argmax(cc) + width//2 + 1\n", + " fit = fitter(Polynomial1D(degree=2), x=lag[low:hi], y=cc[low:hi])\n", + "\n", + " x_c = np.arange(np.min(lag), np.max(lag), 0.01)\n", + "\n", + " # Determining the shift, in pixels\n", + " shift1 = -fit.parameters[1] / (2. * fit.parameters[2])\n", + "\n", + " print(f\"The shift, in pixels, is {np.round(shift1, 2)} pixels\")\n", + " print(f\"The shift, in Å, is {np.round(shift1*mean_plate_scale, 2)}Å\")\n", + "\n", + " if plot:\n", + " plt.plot(\n", + " lag,\n", + " cc,\n", + " lw=2,\n", + " color=\"black\",\n", + " alpha=0.8\n", + " )\n", + "\n", + " plt.plot(\n", + " x_c,\n", + " fit(x_c),\n", + " lw=2,\n", + " label=\"Fitted Quadratic\",\n", + " alpha=0.6\n", + " )\n", + "\n", + " plt.axvline(\n", + " shift1,\n", + " color=\"red\",\n", + " label=\"Maxima of fit\"\n", + " )\n", + "\n", + " # Adding title\n", + " plt.title(\"Lag (px) vs Cross Correlation Coeff, with fitted curve\",\n", + " weight=\"bold\")\n", + "\n", + " plt.xlabel(\"Lag (px)\")\n", + " plt.ylabel(\"Cross Correlation Coeff.\")\n", + "\n", + " plt.ylim(0.9, 1.05)\n", + " plt.xlim(shift1 - 25, shift1 + 25)\n", + "\n", + " plt.legend()\n", + "\n", + " plt.tight_layout()\n", + "\n", + " plt.show()\n", + "\n", + " if run_calstis:\n", + " # Path to the FLT for our example file\n", + " flt_file = f\"{stis_data}/{fits.getval(stis_file, 'ROOTNAME')}_flt.fits\"\n", + "\n", + " # Extracting the shift from the FLT file\n", + " SHIFTA1 = fits.getval(flt_file, \"SHIFTA1\", ext=1)\n", + " SHIFTA2 = fits.getval(flt_file, \"SHIFTA2\", ext=1)\n", + "\n", + " # Update SHIFTA1 (spectral direction) to include the calculated shift\n", + " SHIFTA1 += shift1\n", + "\n", + " # Path to the raw file\n", + " raw_file = f\"{stis_data}/{fits.getval(stis_file, 'ROOTNAME')}_raw.fits\"\n", + "\n", + " # Updating the header to reflect updated shift + FLT SHIFTA2 value\n", + " fits.setval(raw_file, keyword=\"SHIFTA1\", value=SHIFTA1, ext=1)\n", + " fits.setval(raw_file, keyword=\"SHIFTA2\", value=SHIFTA2, ext=1)\n", + "\n", + " # Set WAVECORR to OMIT since we're applying the shift\n", + " fits.setval(raw_file, \"WAVECORR\", value=\"OMIT\", ext=0)\n", + "\n", + " stistools.calstis.calstis(raw_file,\n", + " verbose=False,\n", + " outroot=f\"{str(stis_shifted)}/\")" + ] + }, + { + "cell_type": "markdown", + "id": "f23d5163", + "metadata": {}, + "source": [ + "Let's shift `od1al4040_sx1.fits`:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f078dfc1", + "metadata": {}, + "outputs": [], + "source": [ + "%%capture\n", + "width = 7\n", + "shift_stis(ref_file, e1_pos[1], width, plot=False, run_calstis=True)" + ] + }, + { + "cell_type": "markdown", + "id": "2c4a08e0", + "metadata": {}, + "source": [ + "Now `od1am10e0_sx1.fits`:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "83a76f5b", + "metadata": {}, + "outputs": [], + "source": [ + "%%capture\n", + "width = 9\n", + "shift_stis(ref_file, e1_pos[2], width, plot=False, run_calstis=True)" + ] + }, + { + "cell_type": "markdown", + "id": "ac6018c9", + "metadata": {}, + "source": [ + "And finally, `od1al1040_sx1.fits`:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "147565e3", + "metadata": {}, + "outputs": [], + "source": [ + "%%capture\n", + "width = 7\n", + "shift_stis(ref_file, e1_pos[3], width, plot=False, run_calstis=True)" + ] + }, + { + "cell_type": "markdown", + "id": "df08a39c", + "metadata": {}, + "source": [ + "A file taken in the nominal position, `od1am10b0` is also shifted, so we will fix that too:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e2672bab", + "metadata": {}, + "outputs": [], + "source": [ + "%%capture\n", + "width = 15\n", + "filepath = f\"{stis_unshifted}/od1am10b0_sx1.fits\"\n", + "\n", + "shift_stis(ref_file, filepath, width, plot=False, run_calstis=True)" + ] + }, + { + "cell_type": "markdown", + "id": "ef3ca05b", + "metadata": {}, + "source": [ + "\n", + "## 2.5 Comparing Non-Adjusted and Adjusted STIS Coadds\n", + "\n", + "Now that we've applied the shifts to our data, we will re-run the coadd code. To get a proper comparison, we will need to copy over the files taken at the nominal position to the `stis_shifted` directory since the coadd code needs all `X1D`/`SX1` files in the same directory." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b467bd2c", + "metadata": {}, + "outputs": [], + "source": [ + "for file in sx1s:\n", + " if \"od1am10b0_sx1.fits\" in file:\n", + " continue\n", + "\n", + " if fits.getval(file, \"PROPAPER\") == \"52X2\":\n", + " new_path = f\"{stis_shifted}/{os.path.basename(file)}\"\n", + " shutil.copy(file, new_path)" + ] + }, + { + "cell_type": "markdown", + "id": "9d5dffef", + "metadata": {}, + "source": [ + "And now we run the coadd script:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "34a1f456", + "metadata": {}, + "outputs": [], + "source": [ + "!swrapper -i ./stis_ex/stis_shifted -o ./stis_ex/stis_shifted/coadd" + ] + }, + { + "cell_type": "markdown", + "id": "005e24cc", + "metadata": {}, + "source": [ + "And we can now compare our pre-shifted and shifted data and their coadds." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4da9d1e6", + "metadata": {}, + "outputs": [], + "source": [ + "coadd = \"hst_14423_stis_agkp81d266_g430l_od1a_cspec.fits\"\n", + "\n", + "fig, ax = plt.subplots(2, 2, figsize=(14, 10))\n", + "\n", + "types = [stis_unshifted, stis_shifted]\n", + "\n", + "for i, is_shifted in enumerate(types):\n", + " for file in glob.glob(f\"{is_shifted}/*sx1*\"):\n", + " with fits.open(file) as hdul:\n", + " wl = hdul[1].data[\"WAVELENGTH\"].ravel()\n", + " flux = hdul[1].data[\"FLUX\"].ravel()\n", + "\n", + " aperture = hdul[0].header['PROPAPER']\n", + "\n", + " # Differentiating between nominal and E1 position via linestyle\n", + " if aperture == \"52X2\":\n", + " linestyle = \"solid\"\n", + " label = f\"{hdul[0].header['ROOTNAME']} - nom.\"\n", + " else:\n", + " linestyle = \"dotted\"\n", + " label = f\"{hdul[0].header['ROOTNAME']} - E1\"\n", + "\n", + " ax[0][i].plot(\n", + " wl,\n", + " flux,\n", + " lw=2,\n", + " alpha=0.6,\n", + " label=label,\n", + " linestyle=linestyle\n", + " )\n", + "\n", + " ax[1][i].plot(\n", + " wl,\n", + " flux,\n", + " lw=2,\n", + " alpha=0.6,\n", + " label=label,\n", + " linestyle=linestyle\n", + " )\n", + "\n", + " # Plotting the coadd\n", + " with fits.open(f\"{is_shifted}/coadd/{coadd}\") as hdul:\n", + " wl = hdul[1].data[\"WAVELENGTH\"].ravel()\n", + " flux = hdul[1].data[\"FLUX\"].ravel()\n", + "\n", + " ax[0][i].plot(\n", + " wl,\n", + " flux,\n", + " label=\"Coadd\",\n", + " color=\"black\",\n", + " alpha=0.7\n", + " )\n", + "\n", + " ax[1][i].plot(\n", + " wl,\n", + " flux,\n", + " label=\"Coadd\",\n", + " color=\"black\",\n", + " alpha=0.7\n", + " )\n", + "\n", + " # All plots will have same labels\n", + " for subplot in ax[i].ravel():\n", + " subplot.set_ylabel(r'Flux [$erg\\ s^{-1}\\ cm^{-2}\\ \\AA^{-1}$]')\n", + " subplot.set_xlabel(r'Wavelength [$\\AA$]')\n", + "\n", + " shifted = ['Unshifted', 'Shifted'][i]\n", + "\n", + " ax[0][i].set_title(f\"Flux vs Wavelength, {shifted}\",\n", + " fontweight=\"bold\")\n", + " ax[1][i].set_title(f\"Flux vs Wavelength 4840-4900Å, {shifted}\",\n", + " fontweight=\"bold\")\n", + "\n", + " ax[1][i].set_xlim(4840, 4900)\n", + " ax[1][i].set_ylim(0.4e-13, 1.5e-13)\n", + "\n", + " ax[0][i].legend()\n", + "\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "67a0af3e", + "metadata": {}, + "source": [ + "We can see the offsets between the datasets taken at the nominal position and the `E1` position are reduced after we've applied our shifts and run `CalSTIS`. While we did a very rough estimate of the shift, users will want to verify that they have the correct wavelength solution through their own methods (e.g. checking expected wavelengths to observed wavelengths, comparing radial velocities, etc.)." + ] + }, + { + "cell_type": "markdown", + "id": "8a63bde6-d262-4867-850b-53a537be0465", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "# Congrats on completing the notebook!\n", + "\n", + "### There are more tutorial notebooks for custom coaddition cases in [this](https://github.com/spacetelescope/hst_notebooks/tree/main/notebooks/HASP) repo, check them out!\n", + "\n", + "## About this Notebook\n", + "**Author:** Sierra Gomez (sigomez@stsci.edu)\n", + "\n", + "**Updated on:** July 2, 2024\n", + "\n", + "*This tutorial was generated to be in compliance with the [STScI style guides](https://github.com/spacetelescope/style-guides) and would like to cite the [Jupyter guide](https://github.com/spacetelescope/style-guides/blob/master/templates/example_notebook.ipynb) in particular.*\n", + "\n", + "## Citations\n", + "\n", + "If you use `astropy`, `astroquery`, or `matplotlib`, `numpy` for published research, please cite the authors. Follow these links for more information about citations:\n", + "\n", + "* [Citing `astropy`](https://docs.astropy.org/en/stable/index.html)\n", + "\n", + "* [Citing `astroquery`](https://astroquery.readthedocs.io/en/latest/)\n", + "\n", + "* [Citing `matplotlib`](https://matplotlib.org/stable/users/project/citing.html)\n", + "\n", + "* [Citing `numpy`](https://numpy.org/citing-numpy/)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.8" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/notebooks/HASP/WavelengthAdjustment/pre-requirements.sh b/notebooks/HASP/WavelengthAdjustment/pre-requirements.sh new file mode 100755 index 000000000..4211058bc --- /dev/null +++ b/notebooks/HASP/WavelengthAdjustment/pre-requirements.sh @@ -0,0 +1 @@ +conda install -y -c conda-forge hstcal==2.7.4 \ No newline at end of file diff --git a/notebooks/HASP/WavelengthAdjustment/requirements.txt b/notebooks/HASP/WavelengthAdjustment/requirements.txt new file mode 100644 index 000000000..bd7ecee44 --- /dev/null +++ b/notebooks/HASP/WavelengthAdjustment/requirements.txt @@ -0,0 +1,11 @@ +astropy==6.0.0 +astroquery==0.4.6 +matplotlib==3.8.2 +numpy==1.26.2 +stistools==1.4.4 +crds>=11.17 +ullyses>=4.1.0 +hasp>=0.9.6 +scipy>=1.11.4 +calcos>=3.5.1 +pandas==2.1.4 \ No newline at end of file From 8ab2316d7311b670d6d1cd87d89769bf699e448a Mon Sep 17 00:00:00 2001 From: Sierra Gomez Date: Fri, 12 Jul 2024 13:18:11 -0400 Subject: [PATCH 02/26] Remove capture statements to see if fixes timeout --- .../WavelengthAdjustment/WavelengthAdjustment.ipynb | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb b/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb index 8e8d7eb51..25c86733f 100644 --- a/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb +++ b/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb @@ -233,8 +233,6 @@ "tags": [] }, "source": [ - "**Note: In the cells using `astroquery`, `CalSTIS`, and `CalCOS`, you will see the Jupyter magic command `%%capture`, which keep the output of the cell in memory rather than printing it. This is to reduce the length of the `HTML` notebook. If you would like to see the output of these cells, simply delete the `%%capture` line for each relevant cell and re-run.**\n", - "\n", "\n", "## 1. Performing Wavelength Adjustments on COS Data\n", "\n", @@ -311,7 +309,6 @@ "metadata": {}, "outputs": [], "source": [ - "%%capture\n", "# Querying MAST for our dataset\n", "cos_query = Observations.query_criteria(\n", " obs_id=\"LE2701010\"\n", @@ -461,7 +458,6 @@ "metadata": {}, "outputs": [], "source": [ - "%%capture\n", "ullyses_query = Observations.query_criteria(\n", " target_name=\"LMC079-1\",\n", " provenance_name=\"ULLYSES\",\n", @@ -1401,7 +1397,6 @@ "metadata": {}, "outputs": [], "source": [ - "%%capture\n", "calcos.calcos(\n", " f\"{cos_data}/le2701010_asn.fits\",\n", " outdir=str(cos_shifted),\n", @@ -1800,7 +1795,6 @@ "metadata": {}, "outputs": [], "source": [ - "%%capture\n", "# Querying our datasets\n", "stis_query = Observations.query_criteria(\n", " proposal_id=14423,\n", @@ -2199,7 +2193,6 @@ "metadata": {}, "outputs": [], "source": [ - "%%capture\n", "stistools.calstis.calstis(raw_file,\n", " outroot=f\"{str(stis_shifted)}/\")" ] @@ -2347,7 +2340,6 @@ "metadata": {}, "outputs": [], "source": [ - "%%capture\n", "width = 7\n", "shift_stis(ref_file, e1_pos[1], width, plot=False, run_calstis=True)" ] @@ -2367,7 +2359,6 @@ "metadata": {}, "outputs": [], "source": [ - "%%capture\n", "width = 9\n", "shift_stis(ref_file, e1_pos[2], width, plot=False, run_calstis=True)" ] @@ -2387,7 +2378,6 @@ "metadata": {}, "outputs": [], "source": [ - "%%capture\n", "width = 7\n", "shift_stis(ref_file, e1_pos[3], width, plot=False, run_calstis=True)" ] @@ -2407,7 +2397,6 @@ "metadata": {}, "outputs": [], "source": [ - "%%capture\n", "width = 15\n", "filepath = f\"{stis_unshifted}/od1am10b0_sx1.fits\"\n", "\n", From 91b3265a9fab4feb8defd69a4c83a4b0fc25cb19 Mon Sep 17 00:00:00 2001 From: Sierra Gomez Date: Fri, 12 Jul 2024 13:42:32 -0400 Subject: [PATCH 03/26] Uncomment env var cell --- .../WavelengthAdjustment/WavelengthAdjustment.ipynb | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb b/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb index 25c86733f..e1ffc4c60 100644 --- a/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb +++ b/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb @@ -148,11 +148,11 @@ "metadata": {}, "outputs": [], "source": [ - "# os.environ['CRDS_SERVER_URL'] = 'https://hst-crds.stsci.edu'\n", - "# os.environ['CRDS_SERVER'] = 'https://hst-crds.stsci.edu'\n", - "# os.environ['CRDS_PATH'] = './crds_cache'\n", - "# os.environ['lref'] = './crds_cache/references/hst/cos/'\n", - "# os.environ['oref'] = './crds_cache/references/hst/stis/'" + "os.environ['CRDS_SERVER_URL'] = 'https://hst-crds.stsci.edu'\n", + "os.environ['CRDS_SERVER'] = 'https://hst-crds.stsci.edu'\n", + "os.environ['CRDS_PATH'] = './crds_cache'\n", + "os.environ['lref'] = './crds_cache/references/hst/cos/'\n", + "os.environ['oref'] = './crds_cache/references/hst/stis/'" ] }, { From b29bfb74dd11a648fe8970a8997096ac3a50bbe8 Mon Sep 17 00:00:00 2001 From: Sierra Gomez Date: Fri, 12 Jul 2024 14:20:34 -0400 Subject: [PATCH 04/26] change hstcal ver --- notebooks/HASP/WavelengthAdjustment/pre-requirements.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/notebooks/HASP/WavelengthAdjustment/pre-requirements.sh b/notebooks/HASP/WavelengthAdjustment/pre-requirements.sh index 4211058bc..cf220e231 100755 --- a/notebooks/HASP/WavelengthAdjustment/pre-requirements.sh +++ b/notebooks/HASP/WavelengthAdjustment/pre-requirements.sh @@ -1 +1 @@ -conda install -y -c conda-forge hstcal==2.7.4 \ No newline at end of file +conda install -y -c conda-forge hstcal \ No newline at end of file From fb734d692bfde6120ef441852734ce58e2f247f8 Mon Sep 17 00:00:00 2001 From: Sierra Gomez Date: Fri, 12 Jul 2024 14:59:05 -0400 Subject: [PATCH 05/26] re-commenting env --- .../WavelengthAdjustment/WavelengthAdjustment.ipynb | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb b/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb index e1ffc4c60..fff8fb6e1 100644 --- a/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb +++ b/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb @@ -148,11 +148,11 @@ "metadata": {}, "outputs": [], "source": [ - "os.environ['CRDS_SERVER_URL'] = 'https://hst-crds.stsci.edu'\n", - "os.environ['CRDS_SERVER'] = 'https://hst-crds.stsci.edu'\n", - "os.environ['CRDS_PATH'] = './crds_cache'\n", - "os.environ['lref'] = './crds_cache/references/hst/cos/'\n", - "os.environ['oref'] = './crds_cache/references/hst/stis/'" + "# os.environ['CRDS_SERVER_URL'] = 'https://hst-crds.stsci.edu'\n", + "# os.environ['CRDS_SERVER'] = 'https://hst-crds.stsci.edu'\n", + "# os.environ['CRDS_PATH'] = './crds_cache'\n", + "# os.environ['lref'] = './crds_cache/references/hst/cos/'\n", + "# os.environ['oref'] = './crds_cache/references/hst/stis/'" ] }, { @@ -2605,7 +2605,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.8" + "version": "3.11.9" } }, "nbformat": 4, From d6ac00abaf10e5f319e2d3b6a01fbb7cb7516d52 Mon Sep 17 00:00:00 2001 From: Sierra Gomez Date: Fri, 12 Jul 2024 15:27:22 -0400 Subject: [PATCH 06/26] add hstcal dependencies --- notebooks/HASP/WavelengthAdjustment/pre-requirements.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/notebooks/HASP/WavelengthAdjustment/pre-requirements.sh b/notebooks/HASP/WavelengthAdjustment/pre-requirements.sh index cf220e231..e6a0a65d8 100755 --- a/notebooks/HASP/WavelengthAdjustment/pre-requirements.sh +++ b/notebooks/HASP/WavelengthAdjustment/pre-requirements.sh @@ -1 +1 @@ -conda install -y -c conda-forge hstcal \ No newline at end of file +conda create -n hstcal==2.7.4 -c conda-forge cmake compilers cfitsio pkgconfig \ No newline at end of file From a6e871c7b622428f896300ee268c417865c9ae18 Mon Sep 17 00:00:00 2001 From: Sierra Gomez Date: Fri, 12 Jul 2024 15:30:56 -0400 Subject: [PATCH 07/26] edited conda command --- notebooks/HASP/WavelengthAdjustment/pre-requirements.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/notebooks/HASP/WavelengthAdjustment/pre-requirements.sh b/notebooks/HASP/WavelengthAdjustment/pre-requirements.sh index e6a0a65d8..fc2ae9296 100755 --- a/notebooks/HASP/WavelengthAdjustment/pre-requirements.sh +++ b/notebooks/HASP/WavelengthAdjustment/pre-requirements.sh @@ -1 +1 @@ -conda create -n hstcal==2.7.4 -c conda-forge cmake compilers cfitsio pkgconfig \ No newline at end of file +conda install hstcal==2.7.4 -c conda-forge cmake compilers cfitsio pkgconfig \ No newline at end of file From e604ad5e38b817f84dfb757c42adffdd0ce8397b Mon Sep 17 00:00:00 2001 From: Sierra Gomez Date: Fri, 12 Jul 2024 15:51:36 -0400 Subject: [PATCH 08/26] trying new install method for hstcal --- notebooks/HASP/WavelengthAdjustment/pre-requirements.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/notebooks/HASP/WavelengthAdjustment/pre-requirements.sh b/notebooks/HASP/WavelengthAdjustment/pre-requirements.sh index fc2ae9296..d7457350c 100755 --- a/notebooks/HASP/WavelengthAdjustment/pre-requirements.sh +++ b/notebooks/HASP/WavelengthAdjustment/pre-requirements.sh @@ -1 +1 @@ -conda install hstcal==2.7.4 -c conda-forge cmake compilers cfitsio pkgconfig \ No newline at end of file +conda install -c conda-forge cmake compilers cfitsio pkgconfig \ No newline at end of file From 442dbd51696a1424012813808ab2af92b1951464 Mon Sep 17 00:00:00 2001 From: Sierra Gomez Date: Fri, 12 Jul 2024 15:51:55 -0400 Subject: [PATCH 09/26] trying new install method for hstcal, with version num --- notebooks/HASP/WavelengthAdjustment/pre-requirements.sh | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/notebooks/HASP/WavelengthAdjustment/pre-requirements.sh b/notebooks/HASP/WavelengthAdjustment/pre-requirements.sh index d7457350c..d598c0c47 100755 --- a/notebooks/HASP/WavelengthAdjustment/pre-requirements.sh +++ b/notebooks/HASP/WavelengthAdjustment/pre-requirements.sh @@ -1 +1,2 @@ -conda install -c conda-forge cmake compilers cfitsio pkgconfig \ No newline at end of file +conda install -c conda-forge cmake compilers cfitsio pkgconfig +conda install -y -c conda-forge hstcal==2.7.4 \ No newline at end of file From ccd1f9e9345b8be41c7b8915cb6deec2321bb4f0 Mon Sep 17 00:00:00 2001 From: Sierra Gomez Date: Fri, 12 Jul 2024 15:56:06 -0400 Subject: [PATCH 10/26] add -y to line --- notebooks/HASP/WavelengthAdjustment/pre-requirements.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/notebooks/HASP/WavelengthAdjustment/pre-requirements.sh b/notebooks/HASP/WavelengthAdjustment/pre-requirements.sh index d598c0c47..8757d6d44 100755 --- a/notebooks/HASP/WavelengthAdjustment/pre-requirements.sh +++ b/notebooks/HASP/WavelengthAdjustment/pre-requirements.sh @@ -1,2 +1,2 @@ -conda install -c conda-forge cmake compilers cfitsio pkgconfig +conda install -y -c conda-forge cmake compilers cfitsio pkgconfig conda install -y -c conda-forge hstcal==2.7.4 \ No newline at end of file From 7822e2e101aba1b487d9a4213be660fee7008e84 Mon Sep 17 00:00:00 2001 From: Sierra Gomez Date: Fri, 12 Jul 2024 16:03:54 -0400 Subject: [PATCH 11/26] trying original method --- notebooks/HASP/WavelengthAdjustment/pre-requirements.sh | 1 - 1 file changed, 1 deletion(-) diff --git a/notebooks/HASP/WavelengthAdjustment/pre-requirements.sh b/notebooks/HASP/WavelengthAdjustment/pre-requirements.sh index 8757d6d44..4211058bc 100755 --- a/notebooks/HASP/WavelengthAdjustment/pre-requirements.sh +++ b/notebooks/HASP/WavelengthAdjustment/pre-requirements.sh @@ -1,2 +1 @@ -conda install -y -c conda-forge cmake compilers cfitsio pkgconfig conda install -y -c conda-forge hstcal==2.7.4 \ No newline at end of file From c16964e9ce7e153ee410718b7eab5a43966d1c0a Mon Sep 17 00:00:00 2001 From: Sierra Gomez Date: Fri, 12 Jul 2024 16:29:15 -0400 Subject: [PATCH 12/26] trying newest hstcal version --- notebooks/HASP/WavelengthAdjustment/pre-requirements.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/notebooks/HASP/WavelengthAdjustment/pre-requirements.sh b/notebooks/HASP/WavelengthAdjustment/pre-requirements.sh index 4211058bc..1409f75d6 100755 --- a/notebooks/HASP/WavelengthAdjustment/pre-requirements.sh +++ b/notebooks/HASP/WavelengthAdjustment/pre-requirements.sh @@ -1 +1 @@ -conda install -y -c conda-forge hstcal==2.7.4 \ No newline at end of file +conda install -y -c conda-forge hstcal==3.0.1 \ No newline at end of file From 34b674737665335ea8c4634f3349c568bd17437a Mon Sep 17 00:00:00 2001 From: Sierra Gomez Date: Fri, 12 Jul 2024 16:43:30 -0400 Subject: [PATCH 13/26] using most up-to-date packages --- .../HASP/WavelengthAdjustment/requirements.txt | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/notebooks/HASP/WavelengthAdjustment/requirements.txt b/notebooks/HASP/WavelengthAdjustment/requirements.txt index bd7ecee44..43d1279f1 100644 --- a/notebooks/HASP/WavelengthAdjustment/requirements.txt +++ b/notebooks/HASP/WavelengthAdjustment/requirements.txt @@ -1,11 +1,11 @@ -astropy==6.0.0 -astroquery==0.4.6 -matplotlib==3.8.2 -numpy==1.26.2 +astropy==6.1.1 +astroquery==0.4.8 +matplotlib==3.9.0 +numpy==2.0 stistools==1.4.4 crds>=11.17 ullyses>=4.1.0 -hasp>=0.9.6 -scipy>=1.11.4 +hasp>=1.0.0 +scipy>=1.14.0 calcos>=3.5.1 -pandas==2.1.4 \ No newline at end of file +pandas==2.2.2 \ No newline at end of file From 74ef9f5dec1d6dda87fec1c358d61c5116520cf1 Mon Sep 17 00:00:00 2001 From: Sierra Gomez Date: Fri, 12 Jul 2024 16:47:25 -0400 Subject: [PATCH 14/26] change astroquery v --- notebooks/HASP/WavelengthAdjustment/requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/notebooks/HASP/WavelengthAdjustment/requirements.txt b/notebooks/HASP/WavelengthAdjustment/requirements.txt index 43d1279f1..b6619820d 100644 --- a/notebooks/HASP/WavelengthAdjustment/requirements.txt +++ b/notebooks/HASP/WavelengthAdjustment/requirements.txt @@ -1,5 +1,5 @@ astropy==6.1.1 -astroquery==0.4.8 +astroquery==0.4.7 matplotlib==3.9.0 numpy==2.0 stistools==1.4.4 From ae1f8cec243877d891e1a92a490daa13c936df7c Mon Sep 17 00:00:00 2001 From: Sierra Gomez Date: Fri, 12 Jul 2024 16:56:40 -0400 Subject: [PATCH 15/26] use diff numpy version --- notebooks/HASP/WavelengthAdjustment/requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/notebooks/HASP/WavelengthAdjustment/requirements.txt b/notebooks/HASP/WavelengthAdjustment/requirements.txt index b6619820d..3fb992b73 100644 --- a/notebooks/HASP/WavelengthAdjustment/requirements.txt +++ b/notebooks/HASP/WavelengthAdjustment/requirements.txt @@ -1,7 +1,7 @@ astropy==6.1.1 astroquery==0.4.7 matplotlib==3.9.0 -numpy==2.0 +numpy==1.25.2 stistools==1.4.4 crds>=11.17 ullyses>=4.1.0 From 98388bb77f53581e0a5b8dea4615f72e1f1b5fe7 Mon Sep 17 00:00:00 2001 From: Sierra Gomez Date: Mon, 15 Jul 2024 12:26:39 -0400 Subject: [PATCH 16/26] Updating CRDS cells --- .../WavelengthAdjustment.ipynb | 22 +++++++++++++------ 1 file changed, 15 insertions(+), 7 deletions(-) diff --git a/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb b/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb index fff8fb6e1..68e5bf6f8 100644 --- a/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb +++ b/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb @@ -148,11 +148,11 @@ "metadata": {}, "outputs": [], "source": [ - "# os.environ['CRDS_SERVER_URL'] = 'https://hst-crds.stsci.edu'\n", - "# os.environ['CRDS_SERVER'] = 'https://hst-crds.stsci.edu'\n", - "# os.environ['CRDS_PATH'] = './crds_cache'\n", - "# os.environ['lref'] = './crds_cache/references/hst/cos/'\n", - "# os.environ['oref'] = './crds_cache/references/hst/stis/'" + "os.environ['CRDS_SERVER_URL'] = 'https://hst-crds.stsci.edu'\n", + "os.environ['CRDS_SERVER'] = 'https://hst-crds.stsci.edu'\n", + "os.environ['CRDS_PATH'] = './crds_cache'\n", + "os.environ['lref'] = './crds_cache/references/hst/cos/'\n", + "os.environ['oref'] = './crds_cache/references/hst/stis/'" ] }, { @@ -334,7 +334,11 @@ "\n", "# Now moving all of the original X1Ds to cos_unshifted\n", "for file in glob.glob(f\"{cos_data}/*x1d*\"):\n", - " shutil.move(file, f\"{cos_unshifted}/{os.path.basename(file)}\")" + " shutil.move(file, f\"{cos_unshifted}/{os.path.basename(file)}\")\n", + "\n", + "# Getting the references files for when we run CalCOS\n", + "%%capture\n", + "!crds bestrefs --files ./cos_ex/cos_data/*raw*.fits --sync-references=1 --update-bestrefs" ] }, { @@ -1841,7 +1845,11 @@ " newpath = f\"{stis_unshifted}/{fits.getval(file, 'FILENAME')}\"\n", "\n", " # Moving the actual file\n", - " shutil.move(file, newpath)" + " shutil.move(file, newpath)\n", + "\n", + "# Getting the references files for when we run CalSTIS\n", + "%%capture\n", + "!crds bestrefs --files ./stis_ex/stis_data/*raw*.fits --sync-references=1 --update-bestrefs" ] }, { From 04105f5fe7a1707ae12fe3e5060d599a77d8009b Mon Sep 17 00:00:00 2001 From: Sierra Gomez Date: Mon, 15 Jul 2024 12:37:10 -0400 Subject: [PATCH 17/26] CRDS changes, similar to Setup.ipynb notebook --- .../WavelengthAdjustment.ipynb | 43 ++++++++++++------- 1 file changed, 27 insertions(+), 16 deletions(-) diff --git a/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb b/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb index 68e5bf6f8..eb11f7033 100644 --- a/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb +++ b/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb @@ -138,7 +138,7 @@ "id": "500891d8", "metadata": {}, "source": [ - "Since we will be running both `CalCOS` and `CalSTIS`, we will need to set up environment variables for when we use the [Calibration Reference Data System](https://hst-crds.stsci.edu/). Uncomment these lines if you have not set up environment variables for both COS and STIS:" + "Since we will be running both `CalCOS` and `CalSTIS`, we will need to set up environment variables for when we use the [Calibration Reference Data System](https://hst-crds.stsci.edu/). Much of this is gone in detail in the [COS Setup.ipynb](https://github.com/spacetelescope/hst_notebooks/blob/main/notebooks/COS/Setup/Setup.ipynb) notebook." ] }, { @@ -148,11 +148,30 @@ "metadata": {}, "outputs": [], "source": [ - "os.environ['CRDS_SERVER_URL'] = 'https://hst-crds.stsci.edu'\n", - "os.environ['CRDS_SERVER'] = 'https://hst-crds.stsci.edu'\n", - "os.environ['CRDS_PATH'] = './crds_cache'\n", - "os.environ['lref'] = './crds_cache/references/hst/cos/'\n", - "os.environ['oref'] = './crds_cache/references/hst/stis/'" + "# We will set the CRDS_PATH environment_variable\n", + "# Creating path to where the files are saved\n", + "crds_path = os.path.join(os.environ[\"HOME\"], \"crds_cache\")\n", + "\n", + "# Setting the environment variable CRDS_PATH to our CRDS path\n", + "os.environ[\"CRDS_PATH\"] = crds_path\n", + "\n", + "# URL for the STScI CRDS page\n", + "crds_server_url = \"https://hst-crds.stsci.edu\"\n", + "\n", + "# Setting env variable to URL\n", + "os.environ[\"CRDS_SERVER_URL\"] = crds_server_url\n", + "\n", + "# Set the lref environment variable (COS)\n", + "lref = os.path.join(crds_path, \"references/hst/cos\")\n", + "os.environ[\"lref\"] = lref\n", + "\n", + "!crds sync --contexts hst_cos_0355.imap --fetch-references\n", + "\n", + "# Set the oref environment variable (STIS)\n", + "oref = os.path.join(crds_path, \"references/hst/stis\")\n", + "os.environ[\"oref\"] = oref\n", + "\n", + "!crds sync --contexts hst_cos_0368.imap --fetch-references" ] }, { @@ -334,11 +353,7 @@ "\n", "# Now moving all of the original X1Ds to cos_unshifted\n", "for file in glob.glob(f\"{cos_data}/*x1d*\"):\n", - " shutil.move(file, f\"{cos_unshifted}/{os.path.basename(file)}\")\n", - "\n", - "# Getting the references files for when we run CalCOS\n", - "%%capture\n", - "!crds bestrefs --files ./cos_ex/cos_data/*raw*.fits --sync-references=1 --update-bestrefs" + " shutil.move(file, f\"{cos_unshifted}/{os.path.basename(file)}\")" ] }, { @@ -1845,11 +1860,7 @@ " newpath = f\"{stis_unshifted}/{fits.getval(file, 'FILENAME')}\"\n", "\n", " # Moving the actual file\n", - " shutil.move(file, newpath)\n", - "\n", - "# Getting the references files for when we run CalSTIS\n", - "%%capture\n", - "!crds bestrefs --files ./stis_ex/stis_data/*raw*.fits --sync-references=1 --update-bestrefs" + " shutil.move(file, newpath)" ] }, { From ceac731fe6d69749b12b1fc308d8913b25ecd22a Mon Sep 17 00:00:00 2001 From: Hatice Karatay Date: Wed, 24 Jul 2024 15:31:19 -0400 Subject: [PATCH 18/26] Change the crds files to test the CI run --- notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb b/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb index eb11f7033..5b168cc23 100644 --- a/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb +++ b/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb @@ -171,7 +171,7 @@ "oref = os.path.join(crds_path, \"references/hst/stis\")\n", "os.environ[\"oref\"] = oref\n", "\n", - "!crds sync --contexts hst_cos_0368.imap --fetch-references" + "!crds sync --contexts hst_stis_0368.imap --fetch-references" ] }, { From f2476b178b54a4baa0817abbe3d3d411a5a8f3ba Mon Sep 17 00:00:00 2001 From: Hatice Karatay <66814693+haticekaratay@users.noreply.github.com> Date: Wed, 24 Jul 2024 15:52:45 -0400 Subject: [PATCH 19/26] Update numpy pinning --- notebooks/HASP/WavelengthAdjustment/requirements.txt | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/notebooks/HASP/WavelengthAdjustment/requirements.txt b/notebooks/HASP/WavelengthAdjustment/requirements.txt index 3fb992b73..8d935b67f 100644 --- a/notebooks/HASP/WavelengthAdjustment/requirements.txt +++ b/notebooks/HASP/WavelengthAdjustment/requirements.txt @@ -1,11 +1,11 @@ astropy==6.1.1 astroquery==0.4.7 matplotlib==3.9.0 -numpy==1.25.2 +numpy>=1.25.2 stistools==1.4.4 crds>=11.17 ullyses>=4.1.0 hasp>=1.0.0 scipy>=1.14.0 calcos>=3.5.1 -pandas==2.2.2 \ No newline at end of file +pandas==2.2.2 From 11b72729b0d4dc2710c40bd2ffdaf5a1d0781a2a Mon Sep 17 00:00:00 2001 From: Hatice Karatay Date: Wed, 24 Jul 2024 16:01:41 -0400 Subject: [PATCH 20/26] Add a minor update to trigger CI --- notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb b/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb index 5b168cc23..7fdd0e84a 100644 --- a/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb +++ b/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb @@ -17,7 +17,7 @@ "\n", "## Learning Goals:\n", "\n", - "By the end of this tutorial, you will learn how to:\n", + "By the end of this tutorial, you will learn how to: \n", "\n", "- Determine if you must make wavelength adjustments based on your coadded product\n", "\n", From 007d80b3043beafbf309c82241ed40fc68b7a31c Mon Sep 17 00:00:00 2001 From: Sierra Rosa Gomez Date: Fri, 26 Jul 2024 15:20:22 -0400 Subject: [PATCH 21/26] Incorporate Hatice's feedback into notebook --- .../WavelengthAdjustment.ipynb | 159 ++++++++++-------- 1 file changed, 90 insertions(+), 69 deletions(-) diff --git a/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb b/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb index 7fdd0e84a..040eef98b 100644 --- a/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb +++ b/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb @@ -85,7 +85,7 @@ "```\n", "This will download the dependencies that are necessary to run this current notebook. You will also need to install `hstcal` to run the STIS portion of this notebook. You can run the command below in your terminal after you've activated your `conda` environment:\n", "```\n", - "conda install -y -c conda-forge hstcal==2.7.4\n", + "conda install -y -c conda-forge hstcal==3.0.1\n", "```\n", "Or, similar to the `requirements.txt` file, you can run:\n", "```\n", @@ -138,7 +138,7 @@ "id": "500891d8", "metadata": {}, "source": [ - "Since we will be running both `CalCOS` and `CalSTIS`, we will need to set up environment variables for when we use the [Calibration Reference Data System](https://hst-crds.stsci.edu/). Much of this is gone in detail in the [COS Setup.ipynb](https://github.com/spacetelescope/hst_notebooks/blob/main/notebooks/COS/Setup/Setup.ipynb) notebook." + "Since we will be running both `CalCOS` and `CalSTIS`, we will need to set up environment variables for when we use the [Calibration Reference Data System](https://hst-crds.stsci.edu/). Much of this is gone in detail in the [COS Setup.ipynb](https://github.com/spacetelescope/hst_notebooks/blob/main/notebooks/COS/Setup/Setup.ipynb) notebook. We will download the STIS reference files in a later cell since we'll need the data downloaded." ] }, { @@ -165,13 +165,7 @@ "lref = os.path.join(crds_path, \"references/hst/cos\")\n", "os.environ[\"lref\"] = lref\n", "\n", - "!crds sync --contexts hst_cos_0355.imap --fetch-references\n", - "\n", - "# Set the oref environment variable (STIS)\n", - "oref = os.path.join(crds_path, \"references/hst/stis\")\n", - "os.environ[\"oref\"] = oref\n", - "\n", - "!crds sync --contexts hst_stis_0368.imap --fetch-references" + "!crds sync --contexts hst_cos_0355.imap --fetch-references" ] }, { @@ -206,36 +200,20 @@ " Consolidate all files to single directory; necessary for HASP script run.\n", " ---------------\n", " Input:\n", - " str data_path : ./mastDownload/HST folders paths; files to be moved here\n", + " Path data_path: Path object to the directory where the files should be moved.\n", " ---------------\n", " Output:\n", " None. Files moved to data_path. ./mastDownload/HST directory is deleted.\n", " '''\n", " # The path to all obs_id folders\n", - " mast_path = f\"{data_path}/mastDownload/HST/\"\n", + " mast_path = data_path / Path(\"mastDownload/HST/\")\n", "\n", " try:\n", - " # Check if mastDownload exists\n", - " if not os.path.exists(mast_path):\n", - " print(f\"Directory {mast_path} doesn't exist.\")\n", - " return\n", - "\n", - " # Get a list of the obs_id paths in mastDownload\n", - " obs_id_dirs = os.listdir(mast_path)\n", - "\n", - " # Iterate through each obs_id folder and move the files\n", - " for obs_id in obs_id_dirs:\n", - " obs_id_path = os.path.join(mast_path, obs_id)\n", - " files = glob.glob(obs_id_path + \"/*fits\")\n", + " for obs_id_path in mast_path.iterdir():\n", + " for file in obs_id_path.glob(\"*fits\"):\n", + " shutil.move(str(file), data_path / file.name)\n", "\n", - " for file in files:\n", - " file_path = Path(file)\n", - " new_path = data_path / file_path.name\n", - " shutil.move(file, new_path)\n", - "\n", - " # Now we can remove the mastDownload directory\n", - " if os.path.exists(mast_path):\n", - " shutil.rmtree(f\"{data_path}/mastDownload\")\n", + " shutil.rmtree(mast_path.parent)\n", "\n", " except Exception as e:\n", " print(f\"An error occurred: {e}\")" @@ -352,8 +330,8 @@ "consolidate_files(cos_data)\n", "\n", "# Now moving all of the original X1Ds to cos_unshifted\n", - "for file in glob.glob(f\"{cos_data}/*x1d*\"):\n", - " shutil.move(file, f\"{cos_unshifted}/{os.path.basename(file)}\")" + "for file in cos_data.glob(\"*x1d*\"):\n", + " shutil.move(str(file), cos_unshifted / file.name)" ] }, { @@ -416,7 +394,7 @@ "coadd_cos = \"hst_15824_cos_lmc079-1_g160m_le27_cspec.fits\"\n", "\n", "# List of X1D files\n", - "x1ds = glob.glob(f\"{cos_unshifted}/*x1d*\")\n", + "x1ds = cos_unshifted.glob(\"*x1d*\")\n", "\n", "plt.figure(figsize=(15, 6))\n", "\n", @@ -436,7 +414,7 @@ " alpha=0.6)\n", "\n", "# Plotting the coadd\n", - "with fits.open(f\"{cos_unshifted_coadd}/{coadd_cos}\") as hdul:\n", + "with fits.open(cos_unshifted_coadd / coadd_cos) as hdul:\n", "\n", " wl_coadd = hdul[1].data[\"WAVELENGTH\"].ravel()\n", " flux_coadd = hdul[1].data[\"FLUX\"].ravel()\n", @@ -496,16 +474,16 @@ " download_dir=str(cos_ex)\n", ")\n", "\n", - "# Going to move file to ./cos_ex and get rid of mastDownload + sub-directories\n", - "f = \"hlsp_ullyses_hst_cos_lmc079-1_g160m_dr7_cspec.fits\"\n", + "download_dir = cos_ex / \"mastDownload\"\n", + "hlsp_dir = download_dir / \"HLSP/hlsp_ullyses_hst_cos_lmc079-1_uv\"\n", + "target_filename = \"hlsp_ullyses_hst_cos_lmc079-1_g160m_dr7_cspec.fits\"\n", "\n", - "hlsp_path = f\"{cos_ex}/{f}\"\n", + "source_path = hlsp_dir / target_filename\n", + "hlsp_path = cos_ex / target_filename\n", "\n", - "shutil.move(f\"{cos_ex}/mastDownload/HLSP/hlsp_ullyses_hst_cos_lmc079-1_uv/{f}\",\n", - " hlsp_path)\n", + "shutil.move(str(source_path), str(hlsp_path))\n", "\n", - "# Now we will remove mastDownload and all the subdirectories\n", - "shutil.rmtree(f\"{cos_ex}/mastDownload\")" + "shutil.rmtree(cos_ex / \"mastDownload\")" ] }, { @@ -849,7 +827,10 @@ "source": [ "fig, ax = plt.subplots(3, 1, figsize=(14, 15))\n", "\n", + "x1ds = cos_unshifted.glob(\"*x1d*\")\n", + "\n", "for i, file in enumerate(sorted(x1ds)):\n", + "\n", " # Getting the flux and wavelength for our dataset\n", " data = fits.open(file)[1].data\n", " wavelength = data[\"WAVELENGTH\"].ravel()\n", @@ -1000,6 +981,8 @@ "source": [ "fig, axes = plt.subplots(2, 2, figsize=(16, 16))\n", "\n", + "x1ds = list(cos_unshifted.glob(\"*x1d*\"))\n", + "\n", "# A placeholder reference file\n", "ref_file = x1ds[0]\n", "\n", @@ -1241,7 +1224,7 @@ "metadata": {}, "outputs": [], "source": [ - "fppos1 = f\"{cos_unshifted}/le2701q0q_x1d.fits\"\n", + "fppos1 = cos_unshifted / \"le2701q0q_x1d.fits\"\n", "\n", "widths1 = [25, 25]\n", "\n", @@ -1275,7 +1258,7 @@ "metadata": {}, "outputs": [], "source": [ - "fppos2 = f\"{cos_unshifted}/le2701q2q_x1d.fits\"\n", + "fppos2 = cos_unshifted / \"le2701q2q_x1d.fits\"\n", "\n", "widths2 = [25, 35]\n", "\n", @@ -1309,7 +1292,7 @@ "metadata": {}, "outputs": [], "source": [ - "fppos3 = f\"{cos_unshifted}/le2701q4q_x1d.fits\"\n", + "fppos3 = cos_unshifted / \"le2701q4q_x1d.fits\"\n", "\n", "widths3 = [25, 30]\n", "\n", @@ -1376,7 +1359,7 @@ }, "outputs": [], "source": [ - "with open(f\"{cos_data}/shifts.txt\", \"w\") as shift_file:\n", + "with open(cos_data / \"shifts.txt\", \"w\") as shift_file:\n", " widths = (widths1, widths2, widths3)\n", " for file, width in zip([fppos1, fppos2, fppos3], widths):\n", "\n", @@ -1440,7 +1423,7 @@ "fig, ax = plt.subplots(2, 2, figsize=(16, 8))\n", "\n", "# Plotting the unshifted data on the top row\n", - "for i, file in enumerate(sorted(glob.glob(f\"{cos_unshifted}/*x1d.fits*\"))):\n", + "for i, file in enumerate(sorted(cos_unshifted.glob(\"*x1d.fits*\"))):\n", "\n", " with fits.open(file) as hdul:\n", " wavelength = hdul[1].data[\"WAVELENGTH\"].flatten()\n", @@ -1457,7 +1440,7 @@ " )\n", "\n", "# Now plotting the shifted data on the bottom row\n", - "for i, file in enumerate(sorted(glob.glob(f\"{cos_shifted}/*x1d.fits*\"))):\n", + "for i, file in enumerate(sorted(cos_shifted.glob(\"*x1d.fits*\"))):\n", " with fits.open(file) as hdul:\n", "\n", " wavelength = hdul[1].data[\"WAVELENGTH\"].flatten()\n", @@ -1541,8 +1524,8 @@ "coadd = \"hst_15824_cos_lmc079-1_g160m_le27_cspec.fits\"\n", "\n", "# The paths to both of our coadds\n", - "preshift = f\"{cos_unshifted_coadd}/{coadd}\"\n", - "shifted = f\"{cos_shifted_coadd}/{coadd}\"\n", + "preshift = cos_unshifted_coadd / coadd\n", + "shifted = cos_shifted_coadd / coadd\n", "\n", "# Getting the data for our pre- and post-shifted coadd\n", "with fits.open(preshift) as hdul:\n", @@ -1657,7 +1640,7 @@ "metadata": {}, "outputs": [], "source": [ - "new_coadd = f\"{cos_shifted_coadd}/hst_15824_cos_lmc079-1_g160m_le27_cspec.fits\"\n", + "new_coadd = cos_shifted_coadd / \"hst_15824_cos_lmc079-1_g160m_le27_cspec.fits\"\n", "\n", "with fits.open(new_coadd) as hdul:\n", " data_new_coadd = hdul[1].data\n", @@ -1855,14 +1838,49 @@ "consolidate_files(stis_data)\n", "\n", "# Moving all sx1 to stis_unshifted\n", - "for file in glob.glob(f\"{stis_data}/*sx1.fits*\"):\n", + "for file in stis_data.glob(\"*sx1.fits*\"):\n", " # New path for file\n", - " newpath = f\"{stis_unshifted}/{fits.getval(file, 'FILENAME')}\"\n", + " newpath = stis_unshifted / fits.getval(file, 'FILENAME')\n", "\n", " # Moving the actual file\n", " shutil.move(file, newpath)" ] }, + { + "cell_type": "markdown", + "id": "6dd0db27", + "metadata": {}, + "source": [ + "Now we will go and download the reference files needed for running `CalSTIS`:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "43331308", + "metadata": {}, + "outputs": [], + "source": [ + "# We will set the CRDS_PATH environment_variable\n", + "# Creating path to where the files are saved\n", + "crds_path = os.path.join(os.environ[\"HOME\"], \"crds_cache\")\n", + "\n", + "# Setting the environment variable CRDS_PATH to our CRDS path\n", + "os.environ[\"CRDS_PATH\"] = crds_path\n", + "\n", + "# URL for the STScI CRDS page\n", + "crds_server_url = \"https://hst-crds.stsci.edu\"\n", + "\n", + "# Setting env variable to URL\n", + "os.environ[\"CRDS_SERVER_URL\"] = crds_server_url\n", + "\n", + "# Set the oref environment variable (STIS)\n", + "oref = os.path.join(crds_path, \"references/hst/stis/\")\n", + "os.environ[\"oref\"] = oref\n", + "\n", + "!crds bestrefs --update-bestrefs --sync-references=1 --files ./stis_ex/stis_data/*_raw.fits" + ] + }, { "cell_type": "markdown", "id": "66ef8215-5ace-4fc8-9757-6cb276353975", @@ -1926,7 +1944,7 @@ "coadd_stis = \"hst_14423_stis_agkp81d266_g430l_od1a_cspec.fits\"\n", "\n", "# A list of the SX1 files for our data\n", - "sx1s = glob.glob(f\"{stis_unshifted}/*sx1*\")\n", + "sx1s = list(stis_unshifted.glob(\"*sx1*\"))\n", "\n", "# Creating a figure with two subplots\n", "fig, ax = plt.subplots(2, 1, figsize=(14, 12))\n", @@ -1956,7 +1974,7 @@ " linestyle=linestyle)\n", "\n", " # Plotting the coadd in black\n", - " with fits.open(f\"{stis_unshifted_coadd}/{coadd_stis}\") as hdul:\n", + " with fits.open(stis_unshifted_coadd / coadd_stis) as hdul:\n", " wl = hdul[1].data[\"WAVELENGTH\"].ravel()\n", " flux = hdul[1].data[\"FLUX\"].ravel()\n", "\n", @@ -2013,10 +2031,14 @@ "\n", "for file in sx1s:\n", " # Getting flux data for current interation's file + current brightest file\n", - " curr_flux = fits.open(file)[1].data[\"FLUX\"].ravel()\n", - " ref_flux = fits.open(ref_file)[1].data[\"FLUX\"].ravel()\n", + " with fits.open(file) as hdul:\n", + " curr_flux = hdul[1].data[\"FLUX\"].ravel()\n", + " aperture = hdul[0].header[\"PROPAPER\"]\n", + "\n", + " with fits.open(ref_file) as ref_hdul:\n", + " ref_flux = ref_hdul[1].data[\"FLUX\"].ravel()\n", "\n", - " aperture = fits.getval(file, \"PROPAPER\")\n", + " \n", "\n", " # If current file's mean flux > ref file, replace ref file\n", " if ((np.mean(curr_flux) > np.mean(ref_flux)) & (aperture == \"52X2\")):\n", @@ -2026,8 +2048,7 @@ " continue\n", "\n", "# Printing out our file\n", - "ref_head = fits.open(ref_file)[0].header\n", - "print(f\"The brightest file is {ref_head['FILENAME']}.\")" + "print(f\"The brightest file is {ref_file.name}.\")" ] }, { @@ -2177,7 +2198,7 @@ "outputs": [], "source": [ "# Path to the FLT for our example file\n", - "flt_file = f\"{stis_data}/{fits.getval(stis_file, 'ROOTNAME')}_flt.fits\"\n", + "flt_file = stis_data / f\"{fits.getval(stis_file, 'ROOTNAME')}_flt.fits\"\n", "\n", "# Extracting the shift from the FLT file\n", "SHIFTA1 = fits.getval(flt_file, \"SHIFTA1\", ext=1)\n", @@ -2187,7 +2208,7 @@ "SHIFTA1 += shift1\n", "\n", "# Path to the raw file\n", - "raw_file = f\"{stis_data}/{fits.getval(stis_file, 'ROOTNAME')}_raw.fits\"\n", + "raw_file = stis_data / f\"{fits.getval(stis_file, 'ROOTNAME')}_raw.fits\"\n", "\n", "# Updating the header to reflect updated shift + FLT SHIFTA2 value\n", "fits.setval(raw_file, keyword=\"SHIFTA1\", value=SHIFTA1, ext=1)\n", @@ -2212,8 +2233,8 @@ "metadata": {}, "outputs": [], "source": [ - "stistools.calstis.calstis(raw_file,\n", - " outroot=f\"{str(stis_shifted)}/\")" + "stistools.calstis.calstis(str(raw_file),\n", + " outroot=str(stis_shifted))" ] }, { @@ -2417,7 +2438,7 @@ "outputs": [], "source": [ "width = 15\n", - "filepath = f\"{stis_unshifted}/od1am10b0_sx1.fits\"\n", + "filepath = stis_unshifted / \"od1am10b0_sx1.fits\"\n", "\n", "shift_stis(ref_file, filepath, width, plot=False, run_calstis=True)" ] @@ -2441,11 +2462,11 @@ "outputs": [], "source": [ "for file in sx1s:\n", - " if \"od1am10b0_sx1.fits\" in file:\n", + " if file.name == \"od1am10b0_sx1.fits\":\n", " continue\n", "\n", " if fits.getval(file, \"PROPAPER\") == \"52X2\":\n", - " new_path = f\"{stis_shifted}/{os.path.basename(file)}\"\n", + " new_path = stis_shifted / os.path.basename(file)\n", " shutil.copy(file, new_path)" ] }, @@ -2489,7 +2510,7 @@ "types = [stis_unshifted, stis_shifted]\n", "\n", "for i, is_shifted in enumerate(types):\n", - " for file in glob.glob(f\"{is_shifted}/*sx1*\"):\n", + " for file in is_shifted.glob(\"*sx1*\"):\n", " with fits.open(file) as hdul:\n", " wl = hdul[1].data[\"WAVELENGTH\"].ravel()\n", " flux = hdul[1].data[\"FLUX\"].ravel()\n", @@ -2590,7 +2611,7 @@ "## About this Notebook\n", "**Author:** Sierra Gomez (sigomez@stsci.edu)\n", "\n", - "**Updated on:** July 2, 2024\n", + "**Updated on:** July 26, 2024\n", "\n", "*This tutorial was generated to be in compliance with the [STScI style guides](https://github.com/spacetelescope/style-guides) and would like to cite the [Jupyter guide](https://github.com/spacetelescope/style-guides/blob/master/templates/example_notebook.ipynb) in particular.*\n", "\n", From 0ef9ff4de146338a9f216eb04de4d69afa7b19bf Mon Sep 17 00:00:00 2001 From: Sierra Rosa Gomez Date: Fri, 26 Jul 2024 15:22:59 -0400 Subject: [PATCH 22/26] fix pep8 --- .../HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb | 4 ---- 1 file changed, 4 deletions(-) diff --git a/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb b/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb index 040eef98b..e4c9e5247 100644 --- a/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb +++ b/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb @@ -111,8 +111,6 @@ "import shutil\n", "from pathlib import Path\n", "\n", - "import glob\n", - "\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", @@ -2038,8 +2036,6 @@ " with fits.open(ref_file) as ref_hdul:\n", " ref_flux = ref_hdul[1].data[\"FLUX\"].ravel()\n", "\n", - " \n", - "\n", " # If current file's mean flux > ref file, replace ref file\n", " if ((np.mean(curr_flux) > np.mean(ref_flux)) & (aperture == \"52X2\")):\n", " ref_file = file\n", From fa559b7c27ec9fc161ff8dea3c995cb3a91ab52f Mon Sep 17 00:00:00 2001 From: Sierra Rosa Gomez Date: Fri, 26 Jul 2024 15:50:16 -0400 Subject: [PATCH 23/26] Add try except statements for when CI does second run of cal pipelines --- .../WavelengthAdjustment.ipynb | 39 ++++++++++++++----- 1 file changed, 29 insertions(+), 10 deletions(-) diff --git a/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb b/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb index e4c9e5247..86540baa6 100644 --- a/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb +++ b/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb @@ -1397,10 +1397,13 @@ "metadata": {}, "outputs": [], "source": [ - "calcos.calcos(\n", - " f\"{cos_data}/le2701010_asn.fits\",\n", - " outdir=str(cos_shifted),\n", - " shift_file=f\"{cos_data}/shifts.txt\")" + "try:\n", + " calcos.calcos(f\"{cos_data}/le2701010_asn.fits\",\n", + " outdir=str(cos_shifted),\n", + " shift_file=f\"{cos_data}/shifts.txt\")\n", + " \n", + "except RuntimeError as e:\n", + " print('An error occured', e)" ] }, { @@ -2229,8 +2232,12 @@ "metadata": {}, "outputs": [], "source": [ - "stistools.calstis.calstis(str(raw_file),\n", - " outroot=str(stis_shifted))" + "try:\n", + " stistools.calstis.calstis(str(raw_file),\n", + " outroot=str(stis_shifted))\n", + " \n", + "except RuntimeError as e:\n", + " print('An error occured', e)" ] }, { @@ -2377,7 +2384,10 @@ "outputs": [], "source": [ "width = 7\n", - "shift_stis(ref_file, e1_pos[1], width, plot=False, run_calstis=True)" + "try:\n", + " shift_stis(ref_file, e1_pos[1], width, plot=False, run_calstis=True)\n", + "except RuntimeError as e:\n", + " print('An error occured', e)" ] }, { @@ -2396,7 +2406,10 @@ "outputs": [], "source": [ "width = 9\n", - "shift_stis(ref_file, e1_pos[2], width, plot=False, run_calstis=True)" + "try:\n", + " shift_stis(ref_file, e1_pos[2], width, plot=False, run_calstis=True)\n", + "except RuntimeError as e:\n", + " print('An error occured', e)" ] }, { @@ -2415,7 +2428,10 @@ "outputs": [], "source": [ "width = 7\n", - "shift_stis(ref_file, e1_pos[3], width, plot=False, run_calstis=True)" + "try:\n", + " shift_stis(ref_file, e1_pos[3], width, plot=False, run_calstis=True)\n", + "except RuntimeError as e:\n", + " print('An error occured', e)" ] }, { @@ -2436,7 +2452,10 @@ "width = 15\n", "filepath = stis_unshifted / \"od1am10b0_sx1.fits\"\n", "\n", - "shift_stis(ref_file, filepath, width, plot=False, run_calstis=True)" + "try:\n", + " shift_stis(ref_file, filepath, width, plot=False, run_calstis=True)\n", + "except RuntimeError as e:\n", + " print('An error occured', e)" ] }, { From 85523330bb08924b4a6d830ddaeb057670932c70 Mon Sep 17 00:00:00 2001 From: Sierra Rosa Gomez Date: Tue, 30 Jul 2024 14:03:39 -0400 Subject: [PATCH 24/26] test adddition --- notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb b/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb index 86540baa6..fc00e1ed8 100644 --- a/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb +++ b/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb @@ -136,7 +136,7 @@ "id": "500891d8", "metadata": {}, "source": [ - "Since we will be running both `CalCOS` and `CalSTIS`, we will need to set up environment variables for when we use the [Calibration Reference Data System](https://hst-crds.stsci.edu/). Much of this is gone in detail in the [COS Setup.ipynb](https://github.com/spacetelescope/hst_notebooks/blob/main/notebooks/COS/Setup/Setup.ipynb) notebook. We will download the STIS reference files in a later cell since we'll need the data downloaded." + "Since we will be running both `CalCOS` and `CalSTIS`, we will need to set up environment variables for when we use the [Calibration Reference Data System](https://hst-crds.stsci.edu/). Much of this is gone in detail in the COS notebook, [Setup.ipynb](https://github.com/spacetelescope/hst_notebooks/blob/main/notebooks/COS/Setup/Setup.ipynb) notebook. We will download the STIS reference files in a later cell since we'll need the data downloaded." ] }, { From fda5d02d46d695137e12ebfccbe65f435b117f56 Mon Sep 17 00:00:00 2001 From: Sierra Rosa Gomez Date: Thu, 29 Aug 2024 11:49:33 -0400 Subject: [PATCH 25/26] Update dependencies and context files --- .../WavelengthAdjustment/WavelengthAdjustment.ipynb | 8 ++++---- .../HASP/WavelengthAdjustment/pre-requirements.sh | 2 +- notebooks/HASP/WavelengthAdjustment/requirements.txt | 10 +++++----- 3 files changed, 10 insertions(+), 10 deletions(-) diff --git a/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb b/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb index fc00e1ed8..41ca7c80b 100644 --- a/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb +++ b/notebooks/HASP/WavelengthAdjustment/WavelengthAdjustment.ipynb @@ -85,7 +85,7 @@ "```\n", "This will download the dependencies that are necessary to run this current notebook. You will also need to install `hstcal` to run the STIS portion of this notebook. You can run the command below in your terminal after you've activated your `conda` environment:\n", "```\n", - "conda install -y -c conda-forge hstcal==3.0.1\n", + "conda install -y -c conda-forge hstcal==3.0.3\n", "```\n", "Or, similar to the `requirements.txt` file, you can run:\n", "```\n", @@ -163,7 +163,7 @@ "lref = os.path.join(crds_path, \"references/hst/cos\")\n", "os.environ[\"lref\"] = lref\n", "\n", - "!crds sync --contexts hst_cos_0355.imap --fetch-references" + "!crds sync --contexts hst_cos_0357.imap --fetch-references" ] }, { @@ -2626,7 +2626,7 @@ "## About this Notebook\n", "**Author:** Sierra Gomez (sigomez@stsci.edu)\n", "\n", - "**Updated on:** July 26, 2024\n", + "**Updated on:** August 29, 2024\n", "\n", "*This tutorial was generated to be in compliance with the [STScI style guides](https://github.com/spacetelescope/style-guides) and would like to cite the [Jupyter guide](https://github.com/spacetelescope/style-guides/blob/master/templates/example_notebook.ipynb) in particular.*\n", "\n", @@ -2660,7 +2660,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.9" + "version": "3.11.8" } }, "nbformat": 4, diff --git a/notebooks/HASP/WavelengthAdjustment/pre-requirements.sh b/notebooks/HASP/WavelengthAdjustment/pre-requirements.sh index 1409f75d6..66ef220ef 100755 --- a/notebooks/HASP/WavelengthAdjustment/pre-requirements.sh +++ b/notebooks/HASP/WavelengthAdjustment/pre-requirements.sh @@ -1 +1 @@ -conda install -y -c conda-forge hstcal==3.0.1 \ No newline at end of file +conda install -y -c conda-forge hstcal==3.0.3 \ No newline at end of file diff --git a/notebooks/HASP/WavelengthAdjustment/requirements.txt b/notebooks/HASP/WavelengthAdjustment/requirements.txt index 8d935b67f..16c904e0e 100644 --- a/notebooks/HASP/WavelengthAdjustment/requirements.txt +++ b/notebooks/HASP/WavelengthAdjustment/requirements.txt @@ -1,11 +1,11 @@ -astropy==6.1.1 -astroquery==0.4.7 -matplotlib==3.9.0 +astropy>=6.1.1 +astroquery>=0.4.7 +matplotlib>=3.9.0 numpy>=1.25.2 -stistools==1.4.4 +stistools>=1.4.4 crds>=11.17 ullyses>=4.1.0 hasp>=1.0.0 scipy>=1.14.0 calcos>=3.5.1 -pandas==2.2.2 +pandas>=2.2.2 From 88dcb025644e1060023a57138d5780f3ab288ba2 Mon Sep 17 00:00:00 2001 From: Sierra Rosa Gomez Date: Fri, 30 Aug 2024 10:14:09 -0400 Subject: [PATCH 26/26] remove pinning of hstcal --- notebooks/HASP/WavelengthAdjustment/pre-requirements.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/notebooks/HASP/WavelengthAdjustment/pre-requirements.sh b/notebooks/HASP/WavelengthAdjustment/pre-requirements.sh index 66ef220ef..cf220e231 100755 --- a/notebooks/HASP/WavelengthAdjustment/pre-requirements.sh +++ b/notebooks/HASP/WavelengthAdjustment/pre-requirements.sh @@ -1 +1 @@ -conda install -y -c conda-forge hstcal==3.0.3 \ No newline at end of file +conda install -y -c conda-forge hstcal \ No newline at end of file