docs/notebook/arcwidget_tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ArcWidget Tutorial"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib widget\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This tutorial provides an introduction to how to use the ArcWidget package to wavelength calibrate a spectroscopic observatoin. The tutorial starts with a 1D spectrum of an observations of an AGN and an Argon lamp.  The spectra have already been extracted from the observed data and are currently in units of pixels and electrons.  The tutorial demonstrates how the ArcWidget class in the [Specidentify](https://github.com/crawfordsm/specidentify) package can be used to wavelength calibrate the dispersion axis.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: leap-second auto-update failed due to the following exception: ValueError(\"Malformed URL: '//anaconda3/envs/specreduce/lib/python3.8/site-packages/astropy/utils/iers/data/Leap_Second.dat'\") [astropy.time.core]\n"
     ]
    }
   ],
   "source": [
    "# imports needed from \n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "from astropy import units as u\n",
    "from astropy import modeling as mod\n",
    "from astropy.io import fits\n",
    "from astropy.table import Table\n",
    "\n",
    "from specutils import Spectrum1D\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# imports needed from specidentify \n",
    "from specidentify import WavelengthSolution\n",
    "from specidentify import ArcIdentify"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read in the data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The two data sets made availabe are an Argon arc lamp and an observations of an AGN.  These are already 1 dimensional representations of the spectra as they have already been extracted from the original observations.   The files are written as a FITS data table where the first column is the position in pixels and the second column is the flux in electrons. These observations have been modified from their original observations with the Robert Stobie Spectrograph on the Southern African Large Telescope. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "hdu = fits.open('agn.fits')\n",
    "agn = Spectrum1D(spectral_axis=hdu[1].data['unknown'] * u.pixel, flux=hdu[1].data['flux']*u.electron)\n",
    "hdu.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "907aaa7aa1874debadfaeffb97b00c1d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Flux [electron]')"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# a plot of the uncalibrated AGN spectra\n",
    "ax = plt.subplots()[1]  \n",
    "ax.plot(agn.spectral_axis, agn.flux)  \n",
    "ax.set_xlabel(\"Dispersion [pixels]\")  \n",
    "ax.set_ylabel(\"Flux [electron]\")  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "hdu = fits.open('arc_ar_lamp.fits')\n",
    "arc = Spectrum1D(spectral_axis=hdu[1].data['unknown'] * u.pixel, flux=hdu[1].data['flux']*u.electron)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "af3b0b72892048f8ba56e8112cceb086",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Flux [electron]')"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# a plot of the uncalibrated arc spectra\n",
    "ax = plt.subplots()[1]  \n",
    "ax.plot(arc.spectral_axis, arc.flux)  \n",
    "ax.set_xlabel(\"Dispersion [pixels]\")  \n",
    "ax.set_ylabel(\"Flux [electron]\")  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In addition to the observed spectra, a table of known wavelengths for the arc lamp is also required. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "line_table = np.loadtxt('Ar.salt', unpack=True, usecols=(0,1))\n",
    "line_table = Table(line_table.T, names=('wavelength', 'flux'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create an initial wavelength solution "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our goal is to produce a wavelength calibrated spectrum of the AGN.  In order to do so, we must determine the transformation between pixel position and wavelength.  We use the `specidentify.WavelengthSolution` object to represent that transformation with an `Astropy.modeling.Models` object as the transformation.  In this case, we will use a third-order Chebyshev polynomial as the model for the transformation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "xarr = np.arange(len(arc.data))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "ws_init = mod.models.Chebyshev1D(3)\n",
    "ws_init.domain = [xarr.min(), xarr.max()]\n",
    "ws_init = mod.fitting.LinearLSQFitter()(ws_init, xarr, xarr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "ws = WavelengthSolution(x=None, wavelength=None, model=ws_init)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If a set of points that are already match an x-position to a wavelength, these can be supplied in the x and wavelength values and the Wavelength solution can be fit to those values to provide an initial guess for the wavelength solution. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = [358,390,2300,1211,2489,694, 1421]\n",
    "w = [4510.733, 4522.323,5162.285,4806.019,5221.271,4628.441, 4876.261]\n",
    "w = [4515.733, 4527.323,5167.285,4811.019,5226.271,4633.441, 4881.261]\n",
    "ws = WavelengthSolution(x=x, wavelength=w, model=ws_init)\n",
    "ws.fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Likewise, if there exists a model for the transfomration, then that can also be provided to provide a better initial guess. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interactively wavelength calibrate the spectra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `ArcIdentify` task provides an interactive way to identify wavelengths that match features in the observed spectrum.  The task is built on ipython widgets and meant to have similar capabilities to the IRAF identify task. \n",
    "\n",
    "Once started, it will display the observed arc.  If there is a wavelength solution, then the version of the line list simulated to look like the observed arc will also be plotted in gray.  \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5a7c8381911a46eab0d002e7886bb9bc",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "7e3b9bc8abe3443cad10a0a580129260",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(ColorPicker(value='#2d34ff', description='Color for Arc'), BoundedFloatText(value=0.0, descript…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "aw = ArcIdentify(arc, ws, line_table, wavelength_range=[4300, 5500], flatten_order=9)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Via the widgets or key presses, the user has a number of different options for how to interact with the tools.  The user can change the color of the arc, switch the dispersion axis between pixels and wavelength if the wavelength solution is set, turn off the dispaly of the simulated spectra of the line list, or display features that have been matched.\n",
    "\n",
    "By clicking on an area of the plot, the pixel and wavelength position will be updated to where the user has clicked.  The user can enter in their own values for pixel and wavelength, press 'add point', and then that feature will be added to the list of matches between pixel position and wavelength. \n",
    "\n",
    "There are also a number of key press events as well, and these can be see by running the `show_commands` function.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      " x - centroid on line    a - Display spectrum\n",
      " b - identify features   f - fit solution\n",
      " r - redraw spectrum     e - add closest line\n",
      " d - delete feature\n",
      "       \n"
     ]
    }
   ],
   "source": [
    "aw.show_commands()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The 'x' command will centroid on the nearest spectral line. \n",
    "\n",
    "The 'a' command will toggle the spectrum simulated from the line list to display.\n",
    "\n",
    "The 'e' command will add match the closest feature to the closest line and add that to the list of matched features.\n",
    "\n",
    "The 'b' command will automatically match all features to their closest line. \n",
    "\n",
    "The 'd' command will delete a feature.\n",
    "\n",
    "The 'f' command will fit the solution to the existing list of matched features. \n",
    "\n",
    "The 'r' command will redraw the plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The user can also display the fit to look for any points that need to be manual deleted. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a6a45bbf5c58483e9a11f7db308e969d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "7d1da1b0cd9c41fe9949cd7511929d43",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(ColorPicker(value='#FF0011', description='Color for Arc'), Checkbox(value=True, description='Di…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "aw.draw_error()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.10583618385498994"
      ]
     },
     "execution_count": 77,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Displays the rms of the fit\n",
    "aw.ws.sigma(aw.xp, aw.wp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Wavelength calibrate the spectra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once the users has found a satisfactory solution, the solution can be applied to the observed data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "hdu = fits.open('agn.fits')\n",
    "agn = Spectrum1D(spectral_axis=aw.ws(hdu[1].data['unknown']) * u.angstrom, flux=hdu[1].data['flux']*u.electron)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "35a7f9bec34a4e39ad7ff4f23bb46bf0",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Flux [electron]')"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# a plot of the uncalibrated AGN spectra\n",
    "ax = plt.subplots()[1]  \n",
    "ax.plot(agn.spectral_axis, agn.flux)  \n",
    "ax.set_xlabel(\"Dispersion [Angstrom]\")  \n",
    "ax.set_ylabel(\"Flux [electron]\")  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### References"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "iraf \n",
    "\n",
    "salt\n",
    "\n",
    "pysalt\n",
    "\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}