{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Simple WCS with astropy modeling and gwcs\n", "- categories: [astropy,wcs]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Some notes on how to convert a FITS WCS to a WCS from `gwcs` (Generalized WCS) to be used within the JWST pipeline" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import astropy\n", "from astropy.io import fits\n", "import numpy as np\n", "from astropy import units as u" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "from astropy import wcs" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "hdulist = fits.open(\"example_field/iris_sim_gc_filterKN3.fits\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The conventional FITS WCS is defined by keywords in the FITS file\n", "and is automatically parsed by `astropy.wcs.WCS`" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "SIMPLE = T / Written by IDL: Mon Dec 16 16:40:46 2019 \n", "BITPIX = -64 / IEEE double precision floating point \n", "NAXIS = 2 / \n", "NAXIS1 = 4096 / \n", "NAXIS2 = 4096 / \n", "INSTR = 'IRIS ' / \n", "SCALE = 0.00400000 / pixel scale (arcsec) \n", "UNITS = 'electrons' / \n", "COORD_SY= 'C ' / \n", "RADECSYS= 'FK5 ' / \n", "CTYPE1 = 'RA--LINEAR' / \n", "CTYPE2 = 'DEC-LINEAR' / \n", "CUNIT1 = 'deg ' / \n", "CUNIT2 = 'deg ' / \n", "CRPIX1 = 2048.12 / \n", "CRPIX2 = 2048.12 / \n", "CRVAL1 = 265.197723389 / \n", "CRVAL2 = -28.9921894073 / \n", "CDELT1 = 1.11111013731E-06 /0.00400000 pixel scale \n", "CDELT2 = 1.11111013731E-06 /0.00400000 pixel scale \n", "CROTA1 = 0 / \n", "CROTA2 = 0 / " ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "hdulist[0].header[:22]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Cannot make `LINEAR` to work, so let's instead replace it with Gnomonic" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "hdulist[0].header[\"CTYPE1\"] = \"RA---TAN\"\n", "hdulist[0].header[\"CTYPE2\"] = \"DEC--TAN\"" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "WARNING: FITSFixedWarning: RADECSYS= 'FK5 ' / \n", "the RADECSYS keyword is deprecated, use RADESYSa. [astropy.wcs.wcs]\n", "WARNING: FITSFixedWarning: EPOCH = '2019-12-17T00:40:46.48107737302794Z' / \n", "a floating-point value was expected. [astropy.wcs.wcs]\n" ] } ], "source": [ "w = wcs.WCS(hdulist[0])" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "WCSAXES = 2 / Number of coordinate axes \n", "CRPIX1 = 2048.12 / Pixel coordinate of reference point \n", "CRPIX2 = 2048.12 / Pixel coordinate of reference point \n", "CDELT1 = 1.11111013731E-06 / [deg] Coordinate increment at reference point \n", "CDELT2 = 1.11111013731E-06 / [deg] Coordinate increment at reference point \n", "CUNIT1 = 'deg' / Units of coordinate increment and value \n", "CUNIT2 = 'deg' / Units of coordinate increment and value \n", "CTYPE1 = 'RA---TAN' / Right ascension, gnomonic projection \n", "CTYPE2 = 'DEC--TAN' / Declination, gnomonic projection \n", "CRVAL1 = 265.197723389 / [deg] Coordinate value at reference point \n", "CRVAL2 = -28.9921894073 / [deg] Coordinate value at reference point \n", "LONPOLE = 180.0 / [deg] Native longitude of celestial pole \n", "LATPOLE = -28.9921894073 / [deg] Native latitude of celestial pole \n", "MJDREFI = 0.0 / [d] MJD of fiducial time, integer part \n", "MJDREFF = 0.0 / [d] MJD of fiducial time, fractional part \n", "RADESYS = 'FK5' / Equatorial coordinate system \n", "EQUINOX = 2000.0 / [yr] Equinox of equatorial coordinates " ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "w.to_header()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can then convert between the pixel indices and the coordinates in the sky" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "pixcrd = np.array([[0, 0], [0, 4096],[4096, 4096], [4096,0]], dtype=np.float64)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[265.19512288 -28.99446396]\n", " [265.195123 -28.98991285]\n", " [265.20032602 -28.98991285]\n", " [265.20032613 -28.99446396]]\n" ] } ], "source": [ "world = w.wcs_pix2world(pixcrd, 0)\n", "print(world)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can calculate the size of the instrument field of view" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$18.731693 \\; \\mathrm{{}^{\\prime\\prime}}$" ], "text/plain": [ "" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "((-world[0][0]+ world[-1][0]) * u.deg).to(u.arcsec)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$16.383986 \\; \\mathrm{{}^{\\prime\\prime}}$" ], "text/plain": [ "" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "((-world[0][1]+ world[1][1]) * u.deg).to(u.arcsec)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "WCS Keywords\n", "\n", "Number of WCS axes: 2\n", "CTYPE : 'RA---TAN' 'DEC--TAN' \n", "CRVAL : 265.197723389 -28.9921894073 \n", "CRPIX : 2048.12 2048.12 \n", "NAXIS : 4096 4096" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "w" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Create a `gwcs` WCS object\n", "\n", "We want now to use `astropy.modeling` to build a transformation that is equivalent to the FITS WCS transformation defined above" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "from gwcs import WCS" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "from astropy.modeling import models\n", "from astropy import coordinates as coord\n", "from astropy import units as u\n", "from gwcs import coordinate_frames as cf" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "shift_by_crpix = models.Shift(-(hdulist[0].header[\"CRPIX1\"] - 1)*u.pix) & models.Shift(-(hdulist[0].header[\"CRPIX2\"] - 1)*u.pix)\n" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "tan = models.Pix2Sky_TAN()\n", "celestial_rotation = models.RotateNative2Celestial(\n", " hdulist[0].header[\"CRVAL1\"]*u.deg, hdulist[0].header[\"CRVAL2\"]*u.deg, 180*u.deg)\n" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "tan.input_units_equivalencies = {\"x\": u.pixel_scale(hdulist[0].header[\"CDELT1\"] *u.deg/u.pix),\n", " \"y\": u.pixel_scale(hdulist[0].header[\"CDELT2\"] *u.deg/u.pix)}" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "det2sky = shift_by_crpix | tan | celestial_rotation\n", "det2sky.name = \"linear_transform\"" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": [ "detector_frame = cf.Frame2D(name=\"detector\", axes_names=(\"x\", \"y\"),\n", " unit=(u.pix, u.pix))\n", "sky_frame = cf.CelestialFrame(reference_frame=coord.FK5(), name='fk5',\n", " unit=(u.deg, u.deg))" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " From Transform \n", "-------- ----------------\n", "detector linear_transform\n", " fk5 None\n" ] } ], "source": [ "pipeline = [(detector_frame, det2sky),\n", " (sky_frame, None)\n", " ]\n", "wcsobj = WCS(pipeline)\n", "print(wcsobj)" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(,\n", " )" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ " wcsobj(pixcrd[:,0]*u.pix, pixcrd[:,1]*u.pix, with_units=False)" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$18.731693 \\; \\mathrm{{}^{\\prime\\prime}}$" ], "text/plain": [ "" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "((-_[0][0]+ _[0][-1])).to(u.arcsec)" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$16.383986 \\; \\mathrm{{}^{\\prime\\prime}}$" ], "text/plain": [ "" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "((-__[1][0]+ __[1][1])).to(u.arcsec)" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\n", "\n", " [1]: \n", "\n", " [2]: \n", "\n", " [3]: \n", "Parameters:\n", " offset_0 offset_1 lon_3 lat_3 lon_pole_3\n", " pix pix deg deg deg \n", " -------- -------- ------------- ------------------- ----------\n", " -2047.12 -2047.12 265.197723389 -28.992189407300003 180.0)>" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "wcsobj" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "IRIS", "language": "python", "name": "iris" }, "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 }