figure-out-which-wcs-astropy-uses.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.io import fits\n",
    "from astropy.wcs import WCS\n",
    "from astropy.nddata import CCDData"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## This notebook is checking what astropy does when `PC` and `CD` are included in a FITS header\n",
    "\n",
    "The file below has both."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "bad_file = 'TIC_83210867.01-S001-R089-C001-rp-as_pc.new'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "with fits.open(bad_file) as f:\n",
    "    raw_bad = f[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def print_cd_pc(header):\n",
    "    \"\"\"\n",
    "    Print which of PC and CD are in the header.\n",
    "    \"\"\"\n",
    "    cd_pc = [k for k in header.keys() if k.startswith('CD') or k.startswith('PC')]\n",
    "    for k in sorted(cd_pc):\n",
    "        print(f'{k:6} = {header[k]:e}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CD1_1  = 1.522566e-04\n",
      "CD1_2  = -3.573089e-05\n",
      "CD2_1  = 3.575750e-05\n",
      "CD2_2  = 1.522348e-04\n",
      "PC1_1  = 1.522621e-04\n",
      "PC1_2  = -3.564429e-05\n",
      "PC2_1  = 3.566437e-05\n",
      "PC2_2  = 1.521763e-04\n"
     ]
    }
   ],
   "source": [
    "print_cd_pc(raw_bad.header)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Yep, both are there and they are not equal to each other. If they were, there would be no isue."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read this file as a `CCDData`, which will convert the header to a WCS."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: FITSFixedWarning: RADECSYS= 'FK5 ' / Equatorial coordinate system \n",
      "the RADECSYS keyword is deprecated, use RADESYSa. [astropy.wcs.wcs]\n"
     ]
    }
   ],
   "source": [
    "ccdd = CCDData.read(bad_file, format='fits')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is nothing at all special about this test point except that I know it is not `CRPIX` in this header."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_point = [4000, 4000]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This will be handy for deleting either CD or PC in a moment\n",
    "CDs = [k for k in raw_bad.header if k.startswith('CD')]\n",
    "PCs = [k for k in raw_bad.header if k.startswith('PC')]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Make a wcs without PC and confirm that it has no PC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CD1_1  = 1.522566e-04\n",
      "CD1_2  = -3.573089e-05\n",
      "CD2_1  = 3.575750e-05\n",
      "CD2_2  = 1.522348e-04\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: FITSFixedWarning: RADECSYS= 'FK5 ' / Equatorial coordinate system \n",
      "the RADECSYS keyword is deprecated, use RADESYSa. [astropy.wcs.wcs]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "WCS Keywords\n",
       "\n",
       "Number of WCS axes: 2\n",
       "CTYPE : 'RA---TAN-SIP'  'DEC--TAN-SIP'  \n",
       "CRVAL : 297.459006197  34.8301203697  \n",
       "CRPIX : 2055.0  2048.5  \n",
       "CD1_1 CD1_2  : 0.000152256648923  -3.57308903932e-05  \n",
       "CD2_1 CD2_2  : 3.57575046574e-05  0.000152234782375  \n",
       "NAXIS : 4109  4096"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "head_no_pcs = raw_bad.header.copy()\n",
    "for pc in PCs:\n",
    "    del head_no_pcs[pc]\n",
    "\n",
    "print_cd_pc(head_no_pcs)\n",
    "wcs_no_pcs = WCS(head_no_pcs)\n",
    "wcs_no_pcs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Make a wcs without CD and confirm that it has no CD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PC1_1  = 1.522621e-04\n",
      "PC1_2  = -3.564429e-05\n",
      "PC2_1  = 3.566437e-05\n",
      "PC2_2  = 1.521763e-04\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: FITSFixedWarning: RADECSYS= 'FK5 ' / Equatorial coordinate system \n",
      "the RADECSYS keyword is deprecated, use RADESYSa. [astropy.wcs.wcs]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "WCS Keywords\n",
       "\n",
       "Number of WCS axes: 2\n",
       "CTYPE : 'RA---TAN-SIP'  'DEC--TAN-SIP'  \n",
       "CRVAL : 297.459006197  34.8301203697  \n",
       "CRPIX : 2055.0  2048.5  \n",
       "PC1_1 PC1_2  : 0.000152262101123  -3.56442880876e-05  \n",
       "PC2_1 PC2_2  : 3.56643743388e-05  0.000152176346785  \n",
       "CDELT : 1.0  1.0  \n",
       "NAXIS : 4109  4096"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "head_no_cds = raw_bad.header.copy()\n",
    "for cd in CDs:\n",
    "    del head_no_cds[cd]\n",
    "\n",
    "print_cd_pc(head_no_cds)\n",
    "wcs_no_cds = WCS(head_no_cds)\n",
    "wcs_no_cds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare the CCDData-generated coordinates with these two headers\n",
    "\n",
    "This comparison ignores the SIP distortions because I want to look at just the transform, not the distortion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WCS from   CCDData [array(297.73642431), array(35.19632622)]\n",
      "WCS no CD, just PC [array(297.73642431), array(35.19632622)]\n",
      "WCS no PC, just CD [array(297.73620541), array(35.19662204)]\n"
     ]
    }
   ],
   "source": [
    "print('WCS from   CCDData', ccdd.wcs.wcs_pix2world(*test_point, 0))\n",
    "print('WCS no CD, just PC', wcs_no_cds.wcs_pix2world(*test_point, 0))\n",
    "print('WCS no PC, just CD', wcs_no_pcs.wcs_pix2world(*test_point, 0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Conclusion: astropy uses PC, not CD, if both are present"
   ]
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}