From c83a43e40f9ebb271cde5eef155d45d2e9bc3f22 Mon Sep 17 00:00:00 2001 From: Martin Kilbinger Date: Mon, 16 Oct 2023 18:06:55 +0200 Subject: [PATCH] wcs notebook to compare different WCS conventions --- scripts/jupyter/wcs.ipynb | 542 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 542 insertions(+) create mode 100644 scripts/jupyter/wcs.ipynb diff --git a/scripts/jupyter/wcs.ipynb b/scripts/jupyter/wcs.ipynb new file mode 100644 index 00000000..7a194e1b --- /dev/null +++ b/scripts/jupyter/wcs.ipynb @@ -0,0 +1,542 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "871d80f3-6570-4303-aeb3-f1ded8fb9d09", + "metadata": {}, + "source": [ + "# WCS issues" + ] + }, + { + "cell_type": "markdown", + "id": "26bf7743-4ad1-4275-884f-d9b68a0beb98", + "metadata": {}, + "source": [ + "## PV polynomials" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "a76aedde-b4b5-448e-902b-256952e99671", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "\n", + "from copy import copy\n", + "\n", + "from astropy.io import fits\n", + "from astropy.wcs import WCS\n", + "import sip_tpv as stp\n", + "from esutil import wcsutil\n", + "\n", + "from cs_util import canfar" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "43982ec7-8172-4951-96dc-1601f6759cf1", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Target file 2105209p.fits.fz exists, skipping download.\n" + ] + } + ], + "source": [ + "# Get image\n", + "\n", + "id = \"2105209\"\n", + "image_name = f\"{id}p.fits.fz\"\n", + "vos_dir = f\"vos:cfis/pitcairn/\"\n", + "canfar.download(f\"{vos_dir}/{image_name}\", image_name, verbose=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "cf4c69a6-60f2-4e0f-adb6-7791bd6a839d", + "metadata": {}, + "outputs": [], + "source": [ + "# Open file\n", + "dat = fits.open(image_name)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "3a7b5e25-3c79-4771-80f1-fc4f9bdab252", + "metadata": {}, + "outputs": [], + "source": [ + "# Get main image header\n", + "header_0 = dat[0].header" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "71197a6c-2da6-4fec-a195-2f59db937f27", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "WARNING: FITSFixedWarning: RADECSYS= 'FK5 ' / Coordinate system for equinox (FK4/FK5/GAPPT) \n", + "the RADECSYS keyword is deprecated, use RADESYSa. [astropy.wcs.wcs]\n", + "WARNING: FITSFixedWarning: The WCS transformation has more axes (2) than the image it is associated with (0) [astropy.wcs.wcs]\n", + "WARNING: FITSFixedWarning: 'datfix' made the change 'Invalid parameter values: MJD-OBS and DATE-OBS are inconsistent'. [astropy.wcs.wcs]\n" + ] + } + ], + "source": [ + "# Get WCS, warning might be raised\n", + "wcs_0 = WCS(header_0)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "5a62ffcf-3e9d-4a46-ad46-313d75f057e3", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[]" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Show PV polynomials: empty with standard \"TAN\" convention\n", + "wcs_0.wcs.get_pv()" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "73ce8d14-a3b2-4676-9717-6a932a5dc826", + "metadata": {}, + "outputs": [], + "source": [ + "# Get header of CCD extension\n", + "header_1 = fits.getheader(image_name, 1)\n", + "header_2 = fits.getheader(image_name, 2)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "f14e151e-5a8f-4eaa-87e0-0ef3187c4719", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ERROR 6 in wcsset() at line 2594 of file cextern/wcslib/C/wcs.c:\n", + "PV1_5 : Unrecognized coordinate transformation parameter.\n", + "\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "WARNING: FITSFixedWarning: 'celfix' made the change 'PV1_5 : Unrecognized coordinate transformation parameter'. [astropy.wcs.wcs]\n" + ] + } + ], + "source": [ + "# Try to get WCS: will raise error\n", + "try:\n", + " wcs_1 = WCS(header_1)\n", + "except Exception as e:\n", + " print(e)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "1da0a4bf-098b-4ab5-871f-869aaaf398b8", + "metadata": {}, + "outputs": [], + "source": [ + "# Change convention\n", + "header_0[\"CTYPE1\"] = \"RA-TPV\"\n", + "header_0[\"CTYPE2\"] = \"DEC-TPV\"\n", + "\n", + "header_1[\"CTYPE1\"] = \"RA-TPV\"\n", + "header_1[\"CTYPE2\"] = \"DEC-TPV\"" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "11ccbb22-a684-480e-9f54-392338fa5b26", + "metadata": {}, + "outputs": [], + "source": [ + "wcs_0_tpv = WCS(header_0)\n", + "wcs_1_tpv = WCS(header_1)" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "511bf40f-0e96-472c-bd92-327d1401c55e", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[]" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Show PV polynomials: from main header still empty\n", + "wcs_0_tpv.wcs.get_pv()" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "df253e8f-4e46-4f42-bc62-aae431deb7b6", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[(1, 0, -0.006879527402639),\n", + " (1, 1, 1.01650122237),\n", + " (1, 2, 0.007779393963485),\n", + " (1, 3, 0.0),\n", + " (1, 4, 0.001658429411009),\n", + " (1, 5, 0.001232870290071),\n", + " (1, 6, 0.0006421933266491),\n", + " (1, 7, -0.02501061777044),\n", + " (1, 8, -0.001885796980293),\n", + " (1, 9, -0.02476149608525),\n", + " (1, 10, -0.0005079512200027),\n", + " (2, 0, -0.006116094003451),\n", + " (2, 1, 1.01486402431),\n", + " (2, 2, 0.007783886945903),\n", + " (2, 3, 0.0),\n", + " (2, 4, 0.001436380491668),\n", + " (2, 5, 0.001297612914507),\n", + " (2, 6, 0.0006046712116443),\n", + " (2, 7, -0.02474775482412),\n", + " (2, 8, -0.001916401228446),\n", + " (2, 9, -0.02480256398885),\n", + " (2, 10, -0.0004942902025637)]" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Show PV polynomials: now not empty\n", + "wcs_1_tpv.wcs.get_pv()" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "12636f27-ebd0-4fd1-a373-fb2a1c316d00", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/mkilbing/.conda/envs/shapepipe/lib/python3.9/site-packages/sip_tpv/reverse.py:77: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n", + "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n", + " apcoeffs, r, rank, s = np.linalg.lstsq(a, b)\n", + "/home/mkilbing/.conda/envs/shapepipe/lib/python3.9/site-packages/sip_tpv/reverse.py:82: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n", + "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n", + " bpcoeffs, r, rank, s = np.linalg.lstsq(a, b)\n" + ] + } + ], + "source": [ + "header_1_sip = fits.getheader(image_name, 1)\n", + "\n", + "# Transform header from PV to SIP\n", + "stp.pv_to_sip(header_1_sip)" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "15897199-88d8-4c94-a504-406419431b36", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "DEC--TAN-SIP\n", + "DEC-TPV\n" + ] + } + ], + "source": [ + "# Check key words\n", + "print(header_1_sip[\"CTYPE2\"])\n", + "print(header_1[\"CTYPE2\"])" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "id": "140fcdcc-eba6-48e3-aaca-20aad00ed8e9", + "metadata": {}, + "outputs": [], + "source": [ + "wcs_1_sip = WCS(header_1_sip)" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "4cc2169b-ae94-40dc-967d-bb7df9e316c1", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[]" + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "wcs_1_sip.wcs.get_pv()" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "id": "5a1daffb-ec3f-40b7-a8b2-6071e8132065", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "4644 2112\n" + ] + } + ], + "source": [ + "n_pix_x, n_pix_y = dat[1].data.shape\n", + "print(n_pix_x, n_pix_y)" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "id": "163c68e5-8dc4-4b6b-a893-10b01b8cbb7d", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[[2000. 1500.]\n", + " [2322. 1500.]\n", + " [2000. 1056.]\n", + " [2322. 1056.]]\n", + "[2000. 2322. 2000. 2322.]\n", + "[1500. 1500. 1056. 1056.]\n" + ] + } + ], + "source": [ + "# Define some coordinate pairs\n", + "x = np.array([2000, n_pix_x/2])\n", + "y = np.array([1500, n_pix_y/2])\n", + "xx, yy = np.meshgrid(x, y)\n", + "pixargs = np.vstack([xx.reshape(-1), yy.reshape(-1)]).T\n", + "print(pixargs)\n", + "print(pixargs[:,0])\n", + "print(pixargs[:,1])" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "id": "01ac1dda-5e82-4202-b4bb-65ace2054bce", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[205.70844802, 49.69450748],\n", + " [205.72490131, 49.69446538],\n", + " [205.70815285, 49.71724671],\n", + " [205.72460613, 49.71720461]])" + ] + }, + "execution_count": 19, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Transform to world coordinates using PV\n", + "wcs_1_tpv.all_pix2world(pixargs, 1)" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "id": "0f306759-05f1-497f-beb3-f572edfec9a5", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[205.96952336, 49.69219154],\n", + " [205.99483982, 49.69195117],\n", + " [205.96933441, 49.7148398 ],\n", + " [205.9946501 , 49.7145905 ]])" + ] + }, + "execution_count": 20, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Transform to world coordinates using SIP: different results\n", + "wcs_1_sip.all_pix2world(pixargs, 1)" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "id": "fb696e11-c066-4297-96c6-acc5245b2716", + "metadata": {}, + "outputs": [], + "source": [ + "# Transform sip back to pv\n", + "header_1_tpv_back = copy(header_1_sip)\n", + "stp.sip_to_pv(header_1_tpv_back)\n", + "wcs_1_tpv_back = WCS(header_1_tpv_back)" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "id": "209e5bfa-f4f5-4070-9a3b-9057ee69c697", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[205.96952336, 49.69219154],\n", + " [205.99483982, 49.69195117],\n", + " [205.96933441, 49.7148398 ],\n", + " [205.9946501 , 49.7145905 ]])" + ] + }, + "execution_count": 22, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Transform to world coordinates: same as SIP\n", + "wcs_1_tpv_back.all_pix2world(pixargs, 1)" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "id": "49f45f87-dc7d-4187-ba7d-597370a3d20b", + "metadata": {}, + "outputs": [], + "source": [ + "# Compare with Erin Shepdon's utils\n", + "header_1[\"CTYPE1\"] = \"RA- -TPV\"\n", + "header_1[\"CTYPE2\"] = \"DEC- -TPV\"\n", + "wcs1_tpv_es = wcsutil.WCS(header_1)" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "id": "dcee6113-fd9f-4c06-a5ae-e8703ff12b15", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(array([205.96952336, 205.99483982, 205.96933441, 205.9946501 ]),\n", + " array([49.69219154, 49.69195117, 49.7148398 , 49.7145905 ]))" + ] + }, + "execution_count": 25, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Transform to world coordinates: same as SIP\n", + "wcs1_tpv_es.image2sky(pixargs[:,0], pixargs[:,1])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3f75ed37-6ac1-48b4-b2fd-590c2e209224", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "shapepipe", + "language": "python", + "name": "shapepipe" + }, + "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.9.10" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}