ASDF_and_World_Coordinate_Systems.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ASDF and World Coordinate Systems\n",
    "\n",
    "The handling of World Coordinate Systems (WCS) is inextricably bound to \n",
    "digital observational data. For example, given a pixel in the data array \n",
    "from a detector, what RA, Dec does it map to? Or what pixel does an \n",
    "RA, Dec map to? Likewise for spectrographic data, what wavelength, \n",
    "RA, and Dec does a pixel correspond to? Astronomers have had various \n",
    "solutions to this problem in the past. Most of them have been awkward\n",
    "to deal with.\n",
    "\n",
    "Dealing with WCS issues flexibly was a major driver for ASDF and it \n",
    "deserves its own tutorial.\n",
    "\n",
    "## FITS and WCS\n",
    "\n",
    "FITS has had a series of papers and resulting standards to deal with\n",
    "WCS issues. Generally speaking, these have presumed mostly issues of \n",
    "idealized detectors and optics, where distortion is not a significant\n",
    "factor. There was much focus on handling many kinds of map projections\n",
    "for images, but little provision for distortions beyond what projections\n",
    "introduce. Likewise, most of the spectrosopic WCS models either presume\n",
    "specific dispersion models, or the old fallback of wavelength tables,\n",
    "which to this day many astronomers consider the only solution to \n",
    "spectroscopic WCS.\n",
    "\n",
    "Handling distortions in FITS has met with a failure to settle on \n",
    "any general approach. Currently there are only two practical imaging\n",
    "solutions, namely using the SIP convention, or a pixel translation\n",
    "table. The first is not suitable for high spatial position precision\n",
    "space telescopes \n",
    "such as HST, JWST or NGRST, the second is wasteful of disk space\n",
    "while being inflexible to adjustments (more on that later). Since\n",
    "most spectrographs use effectively the pixel translation approach,\n",
    "they are also wasteful of space for IFUs and long slit data, while\n",
    "also being inflexible to adjustments.\n",
    "\n",
    "## A more flexible approach\n",
    "\n",
    "Ideally, one would want to be able to construct coordinate transforms\n",
    "out of simpler models to form a pipeline of transformations from \n",
    "one coordinate frame to another. The combination of the Python packages - astropy modeling, GWCS (Generalized World Coordinate Systems) and asdf\n",
    "provide a very flexible and powerful means of dealing with WCS.\n",
    "\n",
    "This tutorial will not address how one develops a WCS model as that\n",
    "is almost always a complex issue requiring very deep understanding\n",
    "of the instrument involved and careful treatment of coordinates, and\n",
    "usually is best left to experts (though being able to tweak the \n",
    "parameters of such WCS models is something than can and often should\n",
    "be exposed to users).\n",
    "\n",
    "This tutorial will focus more on using existing WCS models and seeing\n",
    "how they can be modified to improve results.\n",
    "\n",
    "## Using GWCS with ASDF\n",
    "\n",
    "This initial imaging example will use the same wcs that was in the\n",
    "JWST data file used in the tutorial: [Reading a JWST ASDF file]\n",
    "(Reading_a_JWST_ASDF_file.ipynb) except only the WCS has been saved\n",
    "to its own file (ASDF makes saving components of ASDF files to their\n",
    "own files trivial).\n",
    "\n",
    "## Required Software\n",
    "\n",
    "This notebook requires that asdf 2.7.1 (or later), astropy 4.0.1 (or later) and gwcs 0.14 (or later) are installed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import asdf\n",
    "from astropy.utils.data import download_file\n",
    "import matplotlib\n",
    "matplotlib.rcParams['figure.figsize'] = (10, 10)\n",
    "from matplotlib import pyplot as plt\n",
    "file1_url_1 = 'https://data.science.stsci.edu/redirect/ASDF/asdf_tut_gwcs.asdf'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filepath = download_file(file1_url_1)\n",
    "# Wait until download completes (when 'In [*]' changes the  '* to a number')\n",
    "# Ignore warnings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "af = asdf.open(filepath)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the GWCS for these data\n",
    "wcs = af['wcs']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Illustrate use in converting pixel coordinate to sky\n",
    "wcs(100,200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# But what coordinate system for the outputs? What units?\n",
    "wcs.output_frame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The WCS may have intermediate frames defined as well\n",
    "wcs.pipeline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From this we see that there are three frames defined, the detector frame\n",
    "(e.g., pixel coordinates), the v2v3 frame, and the wolrd coordinate frame,\n",
    "which is the output frame.\n",
    "\n",
    "It is possible to access the corresponding coordinates for any intermediate frame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v2v3coords = wcs.transform('detector','v2v3', 100, 200)\n",
    "v2v3coords"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Units?\n",
    "wcs.pipeline[1].frame.unit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Complete computation to world\n",
    "wcs.transform('v2v3', 'world', *v2v3coords)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Relevance to ASDF?\n",
    "\n",
    "OK, GWCS has some extra features. But why is this part of an ASDF tutorial? Mainly to illustrate that ASDF preserves all the intrinsic information about the WCS and does not restrict transforms to a finite set of tools as FITS WCS does. We will illustrate this by examining this WCS in more detail. First note that ASDF allows an arbitrary number of intermediate frames (especially useful for spectral data), as well as an arbitrary number of different WCS models within a file. ASDF provides a way of saving most of astropy's models, of which there is a large number. More important, it can save astropy compound models. Arbitrarily complex combinations of models can be used and saved in ASDF, including all the projection models that FITS provides.\n",
    "\n",
    "So while most of this tutorial is focused on GWCS, it is ASDF that enables this to be used with data files without forcing the user to obtain reference files and special software to determine the correct WCS as had to be done during most of HST lifetime.\n",
    "\n",
    "Let's look at the details of the first step's transform."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "transform1 = wcs.pipeline[0].transform\n",
    "type(transform1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(transform1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This compound model is sufficiently complicated to require some detailed explanation (but not nearly as complicated as some JWST spectral WCS models!). First a brief explanation of the notation. Inputs and outputs indicate that this model takes two input values (or two arrays of values and generates two output values, i.e., two input coordinates and two output coordinates. Expression is the expression used to combine astropy models into a more complex compound model. To keep the expression concise, each term is replaced by a number and the model associated with that number is listed below. The two operators in particular that need special explanation are & and |. The & operator means that two functions are joined by adding dimensionality. If a Shift function shifts the value of the input by the given value, the joining two Shift functions shift the first input by the first value, and the second input by the second value. It is a way of making two 1D functions a 2D function, albeit with each input coordinate treated independently of the other. The | operator feeds the output of the preceding function into the next. The number of outputs of the preceding function must match the number of inputs of the following function. Finally, the Mapping function acts as a switchboard between inputs and output. For example, Mapping(0, 1, 0, 1) takes two inputs and generates 4 outputs, where outputs 0 and 2 are both equal to the first input, and outputs 1 and 3 are both equal to the second input.\n",
    "\n",
    "Now to summarize in words what this transform does: shifts both input values first by 1 and then by -1024.5. This to make the coordinates relative to the center of the detector and center of that pixel. The next step splits the two input values so that both go to two different 2D polynomials (each of which produces only one value, so the net result is two output values. Then these values go through a similar split and to a second polynomial that is effectively applying a small rotation, followed by a shift in the V2,V3 coordinate system.\n",
    "\n",
    "Now to look at the second step transform."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "transform2 = wcs.pipeline[1].transform\n",
    "print(transform2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This transform is much simpler and consists of a scale change then conversion to cartesian coordinates to perform Euler rotations to a RA, Dec system and then conversion back to spherical coordinates.\n",
    "\n",
    "Note that there is no use of a tangent projection, and such distortions that would result from that have been folded into the polynomial distortion model.\n",
    "\n",
    "## Modifying the WCS\n",
    "\n",
    "This is a somewhat contrived example, but it is to illustrate that one can make customizations to a WCS fairly easily. In this particular case one could change any of the parameters in the transforms and save the new WCS to an ASDF file without any concern for how to do that. Or one can modify the existing transform by adding new elements to the transform. So supposing the existing plate scale used does not include velocity aberration, then one would want to include that as well as a plate scale change, and perhaps not change  the existing plate scale application to be clear there are two different scales applied. \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.modeling.models import Scale\n",
    "va_plate_scale_correction = 1.00008\n",
    "transform3 = (Scale(va_plate_scale_correction) & Scale(va_plate_scale_correction)) | transform2\n",
    "wcs.pipeline[1].transform = transform3\n",
    "af.write_to('modified_wcs.asdf', all_array_storage='inline')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "af2 = asdf.open('modified_wcs.asdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wcs2 = af2['wcs']\n",
    "mtransform2 = wcs2.pipeline[1].transform\n",
    "print(mtransform2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For those interested, here are the actual contents of the asdf file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('modified_wcs.asdf') as asdftext: print(asdftext.read()) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result is fairly extensive and not nearly as easy to interpret as the python printed version. Nevertheless, it demonstrates ASDF is capable of saving arbitrarily complex expressions without the user needing to understand the details of saving or reading those expressions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "As befits a complex topic, we have more challenging exercises. These exercises are more GWCS exercises but GWCS is intimately connected to ASDF.\n",
    "\n",
    "## Exercise 1\n",
    "You decide you would rather obtain galactic coordinates from the WCS. How? Hint: the output of the WCS in this case can be treated as an astropy SkyCoord object, so looking at that documentation is recommended. The answer is quite simple."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 2 (needs some Astropy modeling expertise)\n",
    "\n",
    "Change the pixel shifts in the transform by 1 pixel (direction at your pleasure) for both x and y, and save the new WCS to a new file."
   ]
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}