WCS Comparison.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "from astropy.nddata import NDDataArray\n",
    "import astropy.units as u\n",
    "from astropy.wcs import WCS\n",
    "\n",
    "from gwcs import WCS as GWCS\n",
    "from gwcs import coordinate_frames as cf\n",
    "from astropy.modeling.tabular import Tabular1D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "spectral_axis = np.arange(10) * u.AA\n",
    "flux = np.random.sample(10) * u.Jy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a gwcs\n",
    "coord_frame = cf.CoordinateFrame(naxes=1,\n",
    "                                 axes_type=('SPECTRAL',),\n",
    "                                 axes_order=(0,))\n",
    "spec_frame = cf.SpectralFrame(unit=spectral_axis.unit, \n",
    "                              axes_order=(0,))\n",
    "\n",
    "forward_transform = Tabular1D(np.arange(len(spectral_axis)),\n",
    "                              lookup_table=spectral_axis)\n",
    "forward_transform.inverse = Tabular1D(spectral_axis,\n",
    "                                      lookup_table=np.arange(len(spectral_axis)))\n",
    "\n",
    "gwcs = GWCS(forward_transform=forward_transform,\n",
    "            input_frame=coord_frame,\n",
    "            output_frame=spec_frame)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a fitswcs\n",
    "fitswcs = WCS(header={'CDELT1': 1, 'CRVAL1': 0, 'CUNIT1': 'Angstrom',\n",
    "                      'CTYPE1': 'WAVE', 'RESTFRQ': 1400000000, 'CRPIX1': 1},\n",
    "              naxis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using GWCS\n",
    "\n",
    "Let's create a regular `NDData` object that with the random data and the gwcs and do a simple `pixel_to_world` transformation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1. 2.] Angstrom\n"
     ]
    }
   ],
   "source": [
    "arr = NDDataArray(data=flux, wcs=gwcs)\n",
    "print(arr.wcs.pixel_to_world([1, 2]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Awesome! But while slicing seems to work..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "arr_sliced = arr[0:8]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see an issue when we try and do a `pixel_to_world` transformation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1. 2.]] Angstrom\n"
     ]
    }
   ],
   "source": [
    "print(arr_sliced.wcs.pixel_to_world([1, 2]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note the extra dimesion added."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using FITS WCS\n",
    "\n",
    "We also get weird inconsistencies when using a simple fits wcs. We can create the `NDData` object just the same:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1. 2.] Angstrom\n"
     ]
    }
   ],
   "source": [
    "arr = NDDataArray(data=flux, wcs=fitswcs)\n",
    "print(arr.wcs.pixel_to_world([1, 2]).to(u.AA)) # Note, for some reason, the unit is meters for the returned transformation..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.] Angstrom\n"
     ]
    }
   ],
   "source": [
    "arr_sliced = arr[0:8]\n",
    "print(arr_sliced.wcs.pixel_to_world([1, 2]).to(u.AA))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, slicing seems to work, but we notice that we only get the first element in the world values array... Looking at the code (mainly `astropy/wcs/wcsapi/sliced_low_level_wcs.py`):\n",
    "\n",
    "1. Not shown in this notebook is a fix on line 162 in `pixel_to_world_values` in `SlicedLowLevelWCS` which tries to concatanate a list and int when the pixel values are not passed as an ndarray. The list are the pixel values to be converted to world values; the int is (presumably) the offset due to the start of the slicing operation. Converting the list to an ndarray and adding the start value seems to fix this problem.\n",
    "2. The `pixel_to_world_values` on a gwcs and fits wcs object return different things (line 167). For gwcs, it returns a list with a single ndarray containing the world values (e.g. `[array([1, 2])]`), while in fits wcs, it returns just the ndarray (e.g. `array([1, 2])`). This seems to be important because when `SlicedLowLevelWCS` is initialized the pixel dimensions (and world dimensions) to be kept are calculated based on the correlation matrices. In the case of a 1d array wcs, it's just the 0th dimesion. This `world_keep` value is applied to the return value of the `pixel_to_world_values` call; in the gwcs case, it correctly gets `array([1, 2])`, but for fits wcs, it returns `1`.\n",
    "3. Related to 2, applying the index for the kept pixel dimension is done inside a list comprehension, meaning that for gwcs the final returns world values is `[[1, 2]]`. For fits, this means `[1]`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}