Coordinate_Conversion.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "using WCS\n",
    "using FITSIO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "FITS_pix_to_wcs (generic function with 1 method)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function FITS_pix_to_wcs(path, pix_x, pix_y)\n",
    "    file = FITS(path)\n",
    "\n",
    "    header_image = read_header(file[2])\n",
    "    header_keys  = keys(header_image)\n",
    "\n",
    "    # Keys used during WCS transform:\n",
    "    # cdelt -> [TCDLT'X' TCDLT'Y']\n",
    "    # ctype -> [TCTYP'X' TCTYP'Y']\n",
    "    # crpix -> [TCRPX'X' TCRPX'Y']\n",
    "    # crval -> [TCRVL'X' TCRVL'Y']\n",
    "\n",
    "    keys_useful = []\n",
    "\n",
    "    # To select the useful key types\n",
    "    # for ALL COORDINATES (DET1, DET2, etc... and desired declination/right ascension)\n",
    "    for key in header_keys\n",
    "        if length(key) >= 5\n",
    "            if key[1:5] in [\"TCDLT\" \"TCTYP\" \"TCRPX\" \"TCRVL\" \"TCNULL\"]\n",
    "                append!(keys_useful, [key])\n",
    "            end\n",
    "        end\n",
    "    end\n",
    "\n",
    "    # Now that ALL keys have been found, select the ones corresponding to dec/ra\n",
    "    # to find the key numbers\n",
    "    keys_radec = []\n",
    "    for key in keys_useful\n",
    "        if header_image[key] == \"DEC--TAN\"\n",
    "            append!(keys_radec, [key])\n",
    "        elseif header_image[key] == \"RA---TAN\"\n",
    "            append!(keys_radec, [key])\n",
    "        end\n",
    "    end\n",
    "\n",
    "    @assert length(keys_radec) == 2 \"Incorrect number of keys found\"\n",
    "\n",
    "    # Keys ending in these numbers correspond to the desired degree units\n",
    "    keys_radec_id = replace.(keys_radec, \"TCTYP\", \"\")\n",
    "\n",
    "    # In the next part we assume that the key id has TWO DIGITS, which seems to always be true\n",
    "    # check here anyway\n",
    "    @assert !(false in [parse.(Int, keys_radec_id) .>= 10]) \"Key ID under two digits\"\n",
    "\n",
    "    keys_selected = []\n",
    "    for key in keys_useful\n",
    "        if key[end-1:end] in keys_radec_id # Select just keys ending in the key nunber\n",
    "            append!(keys_selected, [key])\n",
    "        end\n",
    "    end\n",
    "\n",
    "    # Set up the transform\n",
    "    cdelt = [header_image[string(\"TCDLT\", keys_radec_id[1])], header_image[string(\"TCDLT\", keys_radec_id[2])]]\n",
    "    ctype = [header_image[string(\"TCTYP\", keys_radec_id[1])], header_image[string(\"TCTYP\", keys_radec_id[2])]]\n",
    "    crpix = [header_image[string(\"TCRPX\", keys_radec_id[1])], header_image[string(\"TCRPX\", keys_radec_id[2])]]\n",
    "    crval = [header_image[string(\"TCRVL\", keys_radec_id[1])], header_image[string(\"TCRVL\", keys_radec_id[2])]]\n",
    "\n",
    "    wcs = WCSTransform(2;\n",
    "                        cdelt = cdelt,\n",
    "                        ctype = ctype,\n",
    "                        crpix = crpix,\n",
    "                        crval = crval)\n",
    "\n",
    "    # *1.0 to convert to float\n",
    "    pixcoords = [pix_x*1.0;  # x coordinates\n",
    "                 pix_y*1.0]  # y coordinates\n",
    "\n",
    "    # DS9 only saves coordinates to EIGHT significan figures in .reg files\n",
    "    # stick to that for consistency\n",
    "    # Maybe it doesn't...\n",
    "    # signif.(pix_to_world(wcs, pixcoords), 8)\n",
    "\n",
    "    pix_to_world(wcs, pixcoords)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "radec_to_sxgm"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "    radec_to_sxgm(ra_in, dec_in)\n",
    "\n",
    "Convert right-ascension and declination from degrees to sexagesimal (h:m:s)\n",
    "\n",
    "Taken from: http://www.bdnyc.org/2012/10/decimal-deg-to-hms/\n",
    "\"\"\"\n",
    "function radec_to_sxgm(ra_in, dec_in)\n",
    "    ra_sign = sign(ra_in)\n",
    "    ra      = abs(ra_in)\n",
    "    \n",
    "    raH = floor(ra/15)\n",
    "    raM = floor(((ra/15) - raH) * 60)\n",
    "    \n",
    "    raS = ((((ra/15) - raH)*60) - raM)*60\n",
    "    \n",
    "    ra_sxgm = (ra_sign*raH, raM, raS)\n",
    "    \n",
    "    dec_sign = sign(dec_in)\n",
    "    dec      = abs(dec_in)\n",
    "\n",
    "    deg = floor(dec)\n",
    "\n",
    "    decM = abs(floor((dec - deg)*60))\n",
    "\n",
    "    decS = (abs((dec - deg)*60) - decM)*60\n",
    "\n",
    "    dec_sxgm = (dec_sign*deg, decM, decS)\n",
    "    \n",
    "    return ra_sxgm, dec_sxgm\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Selecting pixel coordinates 500, 500 in DS9 gives:\n",
    "\n",
    "FK5 degrees: 187.27823, 2.0521586\n",
    "\n",
    "FK5 sexagesimal: 12:29:06.776, +2:03:07.77\n",
    "\n",
    "Using my code for this file "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Array{Float64,1}:\n",
       " 187.278  \n",
       "   2.05216"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fits_path = path = \"/mnt/hgfs/Ubuntu-MPhys-Shared/data_ftp/nustar/data/obs/00/1/10012001002/event_cl/nu10012001002A01_cl.evt.gz\"\n",
    "\n",
    "coords_degrees = FITS_pix_to_wcs(fits_path, 500, 500)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "187.27824162290239"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coords_degrees[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.0521585961635624"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coords_degrees[2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Within rounding errors these are the same as the decimal values given by DS9\n",
    "\n",
    "Now for sexagesimal ones:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((12.0, 29.0, 6.777989496573724), (2.0, 3.0, 7.770946188824794))"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "radec_to_sxgm(coords_degrees[1], coords_degrees[2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So the right ascension is within rounding, as is the declination\n",
    "\n",
    "For another example:\n",
    "\n",
    "Degrees: 299.59612, 35.225458\n",
    "\n",
    "Sexagesimal: 19:58:23.069, +35:13:31.65"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Array{Float64,1}:\n",
       " 299.596 \n",
       "  35.2255"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fits_path = path = \"/home/robert/Documents/10002003001/pipeline_out/nu10002003001A01_cl.evt\"\n",
    "\n",
    "coords_degrees = FITS_pix_to_wcs(fits_path, 500, 500)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "299.59611793184934"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coords_degrees[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "35.225458595481825"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coords_degrees[2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((19.0, 58.0, 23.068303643836146), (35.0, 13.0, 31.65094373457123))"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "radec_to_sxgm(coords_degrees[1], coords_degrees[2])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.0",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}