From 0c78f98858ff007ef4c3b90fda0b83e4b0ecd293 Mon Sep 17 00:00:00 2001 From: Paul Goudfrooij Date: Tue, 17 Sep 2024 19:55:03 -0400 Subject: [PATCH] adding update_pure_parallel_wcs (#241) * adding update_pure_parallel_wcs * adding one sentence to NIRISS_correct_pure_parallel_WCS.ipynb * adding update_parallel_wcs.py * adding missing whitespace in code * some indents were not multiples of 4 * correcting one more indent * adding one black line (gosh...) * few pep8 errors in script addressed * few more pep8 errors in script addressed * implementing suggestions of reviewer * adding mastquery to requirements * adding mastquery to requirements * clearing output cells * correcting typo * trying to clear one import * Add new nb to the toc --------- Co-authored-by: Hatice Karatay --- _toc.yml | 2 + .../NIRISS_correct_pure_parallel_WCS.ipynb | 374 ++++++++++++++++++ .../update_pure_parallel_wcs/README.md | 30 ++ .../update_pure_parallel_wcs/requirements.txt | 7 + .../update_parallel_wcs.py | 163 ++++++++ 5 files changed, 576 insertions(+) create mode 100644 notebooks/cross_instrument/update_pure_parallel_wcs/NIRISS_correct_pure_parallel_WCS.ipynb create mode 100644 notebooks/cross_instrument/update_pure_parallel_wcs/README.md create mode 100644 notebooks/cross_instrument/update_pure_parallel_wcs/requirements.txt create mode 100644 notebooks/cross_instrument/update_pure_parallel_wcs/update_parallel_wcs.py diff --git a/_toc.yml b/_toc.yml index d0a6aab14..9df60f39c 100644 --- a/_toc.yml +++ b/_toc.yml @@ -49,6 +49,8 @@ parts: title: RGB images with Imviz - file: notebooks/cross_instrument/specviz_notebookGUI_interaction/specviz_notebook_gui_interaction_redshift.ipynb title: Specviz Simple Demo + - file: notebooks/cross_instrument/update_pure_parallel_wcs/NIRISS_correct_pure_parallel_WCS.ipynb + title: Improving Accuracy of WCS of Pure Parallel Dataset - caption: MIRI chapters: - file: notebooks/MIRI/MRS_Mstar_analysis/JWST_Mstar_dataAnalysis_runpipeline.ipynb diff --git a/notebooks/cross_instrument/update_pure_parallel_wcs/NIRISS_correct_pure_parallel_WCS.ipynb b/notebooks/cross_instrument/update_pure_parallel_wcs/NIRISS_correct_pure_parallel_WCS.ipynb new file mode 100644 index 000000000..052773471 --- /dev/null +++ b/notebooks/cross_instrument/update_pure_parallel_wcs/NIRISS_correct_pure_parallel_WCS.ipynb @@ -0,0 +1,374 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "18f3654c-92ae-4074-a19a-158ab5bc8ada", + "metadata": {}, + "source": [ + "# Improving Accuracy of World Coordinate System of Pure Parallel Dataset\n", + "\n", + "## Introduction\n", + "\n", + "In this notebook, we will go through the steps needed to improve the accuracy of the World Coordinate System (WCS) in the headers of pure parallel datasets that were observed with JWST prior to the installation of DMS 10.2 (April 2024). The example dataset is jw01571078001_03201_00001_nis_rate.fits, which is a direct image from pure parallel Program ID 1571 (PI: Malkan). Note that the `update_parallel_wcs.py` script has to be run on all individual datasets taken during a pure parallel visit in order to get the improved WCS for that visit. \n", + "\n", + "In this notebook, we assume that all relevant files are located in the current working directory.\n", + "\n", + "### Install pipeline and other required packages\n", + "\n", + "The required packages are in the provided file `requirements.txt`. We generally recommend to create a fresh conda environment followed by the installation of those required packages:\n", + "\n", + "```\n", + "conda create -n improve_pure_parallel_wcs\n", + "conda activate improve_pure_parallel_wcs\n", + "pip install -r requirements.txt\n", + "```\n", + "\n", + "Date last published: September 17, 2024\n", + "\n", + "### Imports\n", + "\n", + "The imports in the next cell are needed to run this notebook." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "fc6bcd19", + "metadata": {}, + "outputs": [], + "source": [ + "import copy\n", + "import jwst.datamodels as dm\n", + "import numpy as np\n", + "from astropy.io import ascii\n", + "from astropy import units as u\n", + "from astropy.coordinates import SkyCoord\n", + "from astropy.table import Table\n", + "from astropy.time import Time\n", + "from astropy.convolution import interpolate_replace_nans, Gaussian2DKernel\n", + "import astropy.io.fits as pyfits\n", + "import matplotlib.pyplot as plt\n", + "import matplotlib.patches as patches\n", + "import pysiaf\n", + "import warnings\n", + "import urllib.request\n", + "from astroquery.gaia import Gaia\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "markdown", + "id": "78ec2247-e0c8-4c96-9a28-39b090337576", + "metadata": {}, + "source": [ + "### Define a few functions" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "b221f2b6", + "metadata": {}, + "outputs": [], + "source": [ + "# This example is to show the (mis-)alignment of the World Coordinate System of a Pure Parallel imaging dataset \n", + "# taken before the installation of DMS 10.2, relative to GAIA DR3, and how the application of the \n", + "# \"update_parallel_wcs.py\" script improves that alignment.\n", + "#\n", + "# First some functions to obtain a catalog of GAIA stars near a certain sky position:\n", + "#\n", + "def query_gaia(ra, dec, radius, verbose=False, epoch=None, filename=None):\n", + " \"\"\"\n", + " Execute a Gaia DR3 query using astroquery, return the table of sources.\n", + "\n", + " Parameters\n", + " ----------\n", + "\n", + " ra: a float value, the right ascension in decimal degrees\n", + "\n", + " dec: a float value, the declination in decimal degrees\n", + "\n", + " radius: an optional float value, the search radius in arc-seconds, \n", + " default 1.0\n", + "\n", + " epoch: an optional float value, the epoch for the positions in decimal \n", + " years (2023.197588611 for example, for 2021-05-24T17:46:12.814) \n", + " if given, the proper motion is allowed for in the positions that \n", + " are returned\n", + "\n", + " verbose: an optional boolean value, if True print a list of the sources\n", + "\n", + " filename: an optional string value, used as the output file name if \n", + " the verbose flag is set (if None, print only to the terminal)\n", + "\n", + " Returns\n", + " -------\n", + "\n", + " gaiadata1: an Astropy Table containing the catalogue of sources\n", + " \"\"\"\n", + " coord = SkyCoord(ra=ra, dec=dec, unit=(u.degree, u.degree),\n", + " frame='icrs')\n", + " radius = u.Quantity(radius/3600.0, u.deg)\n", + " Gaia.ROW_LIMIT = -1\n", + " Gaia.MAIN_GAIA_TABLE = \"gaiadr3.gaia_source\"\n", + " gaiadata1 = Gaia.query_object_async(coordinate=coord, width=radius*2, height=radius*2)\n", + " if verbose:\n", + " gphot = gaiadata1['phot_g_mean_mag']\n", + " gbphot = gaiadata1['phot_bp_mean_mag']\n", + " grphot = gaiadata1['phot_rp_mean_mag']\n", + " parallax = gaiadata1['parallax']\n", + " names = gaiadata1['DESIGNATION']\n", + " with warnings.catch_warnings():\n", + " warnings.simplefilter(\"ignore\")\n", + " gabs = gphot-5.*np.log10(1000./parallax)+5.\n", + " gcol = gbphot - grphot\n", + " if epoch is None:\n", + " epoch = 2016.0\n", + " gaiadata2 = apply_precession(gaiadata1, epoch-2016.0)\n", + " ra = gaiadata2['ra']\n", + " dec = gaiadata2['dec']\n", + " with warnings.catch_warnings():\n", + " warnings.simplefilter(\"ignore\")\n", + " for ind in range(len(ra)):\n", + " g_col = gcol[ind] if (gbphot[ind] < 90.) and (grphot[ind] < 90.) else 0.0\n", + " line = f\"{ra[ind]:12.8f} {dec[ind]:13.8f} {gphot[ind]:10.6f} {gbphot[ind]:10.6f} {grphot[ind]:10.6f} {parallax[ind]:10.3f} {gabs[ind]:10.4f} {g_col:10.4f} '{names[ind]}'\"\n", + " line = line.replace('nan', '0.0')\n", + " if filename:\n", + " with open(filename, 'w') as outfile:\n", + " print(' RA DEC g gbp grp parallax abs_g gbp_grp GAIA_Designation', file=outfile)\n", + " for ind in range(len(ra)):\n", + " g_col = gcol[ind] if (gbphot[ind] < 90.) and (grphot[ind] < 90.) else 0.0\n", + " line = f\"{ra[ind]:12.8f} {dec[ind]:13.8f} {gphot[ind]:10.6f} {gbphot[ind]:10.6f} {grphot[ind]:10.6f} {parallax[ind]:10.3f} {gabs[ind]:10.4f} {g_col:10.4f} '{names[ind]}'\"\n", + " line = line.replace('nan', '0.0')\n", + " print(line, file=outfile)\n", + " return gaiadata2\n", + "\n", + "\n", + "def apply_precession(catalog, deltat):\n", + "\n", + " \"\"\"\n", + " Apply precession to update a catalog of values. Uses astropy.SkyCoord \n", + " apply_space_motion to update the positions.\n", + "\n", + " Parameters:\n", + "\n", + " catalog: a numpy Table type variable with elements 'ra', 'dec', 'pmra', \n", + " and 'pmdec' as per the Gaia DR3 catalog; sky coordinates \n", + " must be in degrees and the proper motions must be in mas/year\n", + "\n", + " deltat: a float value, the time change in decimal years for the motion\n", + "\n", + " Returns\n", + " -------\n", + "\n", + " newcatalog: a copy of catalog with 'newra', 'newdec' elements with the \n", + " revised sky positions; if no proper motion data are \n", + " available then newra = ra and newdec = dec\n", + " \"\"\"\n", + " newcatalog = copy.deepcopy(catalog)\n", + " newcatalog['newra'] = catalog['ra']\n", + " newcatalog['newdec'] = catalog['dec']\n", + " with warnings.catch_warnings():\n", + " warnings.simplefilter(\"ignore\")\n", + " sky_coords = SkyCoord(catalog['ra'], catalog['dec'], \n", + " unit=(u.deg, u.deg),\n", + " pm_ra_cosdec=catalog['pmra'], \n", + " pm_dec=catalog['pmdec'],\n", + " obstime=Time('2016-01-01 00:00:00.0'))\n", + " newpos = sky_coords.apply_space_motion(dt=deltat*u.yr)\n", + " newra = newpos.ra.value\n", + " newdec = newpos.dec.value\n", + " newcatalog['newra'] = newra\n", + " newcatalog['newdec'] = newdec\n", + " inds = np.isnan(newra)\n", + " newcatalog['newra'][inds] = newcatalog['ra'][inds]\n", + " newcatalog['newdec'][inds] = newcatalog['dec'][inds]\n", + " return newcatalog" + ] + }, + { + "cell_type": "markdown", + "id": "d848cd59-b179-43d3-a56b-c0d803cd8de1", + "metadata": {}, + "source": [ + "### Get sky coordinates of a list of GAIA targets around the pointing of the pure parallel image" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "63028692", + "metadata": {}, + "outputs": [], + "source": [ + "# Now get such a list around pointing of NIRISS pure parallel image jw01571078001_03201_00001_nis_rate.fits\n", + "# downloaded from the MAST archive. \n", + "# Note: files downloaded or created here will be placed in the current working directory.\n", + "boxlink = 'https://stsci.box.com/shared/static/ydxn3hhndwup0qr85fuyqro6suufa6fx.fits'\n", + "boxfile = 'jw01571078001_03201_00002_nis_rate.fits'\n", + "try:\n", + " urllib.request.urlretrieve(boxlink, boxfile)\n", + "except Exception as e:\n", + " print(f\"Error downloading file: {e}\")\n", + "try:\n", + " imfile = pyfits.open(boxfile) \n", + " hdr1 = imfile[1].header\n", + " ra, dec = hdr1['CRVAL1'], hdr1['CRVAL2']\n", + " print(f'RA = {ra}, DEC = {dec}')\n", + "except Exception as e:\n", + " print(f\"Error opening FITS file: {e}\")\n", + "mygaia = query_gaia(ra, dec, 65., verbose=True, epoch=2022.9968, filename='gaiacoords.out')" + ] + }, + { + "cell_type": "markdown", + "id": "973801cd-bd8e-496d-b036-1ccd18533744", + "metadata": {}, + "source": [ + "### Convert those sky coordinates to (x, y) in the image according to its WCS, then show them on image" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "05d01381", + "metadata": {}, + "outputs": [], + "source": [ + "t = ascii.read('gaiacoords.out')\n", + "ra, dec = (t['RA'], t['DEC'])\n", + "ratefile = 'jw01571078001_03201_00002_nis_rate.fits'\n", + "instr = imfile[0].header['INSTRUME']\n", + "aperture = imfile[0].header['APERNAME']\n", + "siaf = pysiaf.Siaf(instr)\n", + "myaper = siaf[aperture]\n", + "mod = dm.open(ratefile)\n", + "am = pysiaf.utils.rotations.attitude_matrix(0, 0, mod.meta.pointing.ra_v1,\n", + " mod.meta.pointing.dec_v1,\n", + " mod.meta.pointing.pa_v3)\n", + "myaper.set_attitude_matrix(am)\n", + "x, y = myaper.sky_to_sci(ra, dec)\n", + "\n", + "xytab = Table([ra, dec, x, y], names=('RA', 'DEC', 'x', 'y'))\n", + "xytab['RA'].info.format = '.8f'\n", + "xytab['DEC'].info.format = '.8f'\n", + "xytab['x'].info.format = '.4f'\n", + "xytab['y'].info.format = '.4f'\n", + "print(xytab)\n", + "\n", + "ys, xs = imfile[1].data.shape\n", + "fig = plt.figure(figsize=(15, 15))\n", + "ax = fig.add_subplot(111)\n", + "xsize, ysize = (20, 20)\n", + "for i in range(len(t)):\n", + " ax.add_patch(patches.Ellipse(\n", + " (x[i], y[i]),\n", + " (xsize),\n", + " (ysize), fill=False, color='red'))\n", + "# Create image with bad pixels interpolated over (for display purposes)\n", + "kernel = Gaussian2DKernel(x_stddev=2)\n", + "fixed_image = interpolate_replace_nans(imfile[1].data, kernel)\n", + "ax.imshow(fixed_image, cmap='binary', origin='lower', extent=[0, xs-1, 0, ys-1], vmin=0.7, vmax=2)" + ] + }, + { + "cell_type": "markdown", + "id": "65f2e855-10c2-4555-9f1e-256b943ab0a6", + "metadata": {}, + "source": [ + "Note that the WCS is 'off' by more than a pixel, which is a problem for spectral extraction in the pipeline for WFSS data, since **that extraction relies blindly on the accuracy of the WCS**. \n", + "\n", + "### Now apply the update_parallel_wcs.py script to correct the WCS. Note the offset in pixel coordinates calculated by the script:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "544edfca", + "metadata": {}, + "outputs": [], + "source": [ + "# This cell assumes that you have script update_parallel_wcs.py in your working directory. \n", + "# Note: Outside of this notebook, this script should be run as follows on each _rate.fits file in the pure parallel visit:\n", + "# \n", + "# python update_parallel_wcs.py my_rate.fits \n", + "# \n", + "# After the script has been run on the input file(s), one can then (re-)run Stage 2 of the JWST calibration pipeline \n", + "# (`calwebb_image2` and/or `calwebb_spec2`) on those files, which will now result in correct WCSes in the data headers \n", + "# (and corrected locations of spectral extractions in case of pure parallel WFSS data).\n", + "#\n", + "# By default, `update_parallel_wcs.py` displays the input and output values of the CRVAL1/2 keywords when the script is run \n", + "# as seen below. One can avoid this by setting the optional parameter to anything other than `True`.\n", + "#\n", + "# The script keeps track of its executions using a log file called \"pure_parallel_wcs_logfile\" in the working directory.\n", + "\n", + "%run update_parallel_wcs.py jw01571078001_03201_00002_nis_rate.fits" + ] + }, + { + "cell_type": "markdown", + "id": "0d48372b-5c1e-45e7-9cda-8ca364f2170c", + "metadata": {}, + "source": [ + "### Calculate (x, y) for the GAIA targets in the updated fits file and print the updated (x, y) coordinates:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b00e16e9", + "metadata": {}, + "outputs": [], + "source": [ + "with pyfits.open('jw01571078001_03201_00002_nis_rate.fits') as imfile:\n", + " hdr0 = imfile[0].header\n", + " myaper = pysiaf.Siaf(hdr0['INSTRUME'])[hdr0['APERNAME']]\n", + " mod = dm.open('jw01571078001_03201_00002_nis_rate.fits')\n", + " am = pysiaf.utils.rotations.attitude_matrix(0, 0, mod.meta.pointing.ra_v1,\n", + " mod.meta.pointing.dec_v1,\n", + " mod.meta.pointing.pa_v3)\n", + " myaper.set_attitude_matrix(am)\n", + " newx, newy = myaper.sky_to_sci(ra, dec)\n", + " xytab = Table([x, y, newx, newy], names=('x', 'y', 'x_corr', 'y_corr'))\n", + " for col in xytab.colnames:\n", + " xytab[col].info.format = '.4f'\n", + " print(xytab)\n", + "\n", + "# In practice, one would run the script on all _rate.fits files for a given target, \n", + "# then run the stage-2 and stage-3 pipelines on the resulting images \n", + "# (specifically, Image2Pipeline for the direct images, followed by Image3Pipeline to create a combined \n", + "# image and source catalog, followed by the stage-2 and stage-3 pipelines for the grism images \n", + "# (Spec2Pipeline and Spec3Pipeline, respectively)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4e3eca3e-eca5-4343-911d-8c5fe37b5ee2", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "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.11.10" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/notebooks/cross_instrument/update_pure_parallel_wcs/README.md b/notebooks/cross_instrument/update_pure_parallel_wcs/README.md new file mode 100644 index 000000000..c43f7a0d2 --- /dev/null +++ b/notebooks/cross_instrument/update_pure_parallel_wcs/README.md @@ -0,0 +1,30 @@ +# Pure_Parallels + +A Python script that can be used to correct the WCS of pure parallel exposures. + +## Available Scripts + +### `update_parallel_wcs.py` + +Update CRVAL1 and CRVAL2 header keywords (which indicate the sky coordinates (RA and Dec) of the reference pixel) of Pure Parallel data. + +Run as follows: +``` +python update_parallel_wcs.py fitsfile +``` +where `fitsfile` is the input FITS file. The typical file that is input to this script is a file that was run through Stage 1 of the JWST calibration pipeline, i.e., a `_rate.fits` or `_rateints.fits` file. After the script has been run on the input file(s), one can then (re-)run Stage 2 of the JWST calibration pipeline (`calwebb_image2` and/or `calwebb_spec2`) on those files, which will now result in correct WCSes in the data headers (and corrected locations of spectral extractions in case of pure parallel WFSS data). + +By default, `update_parallel_wcs.py` displays the input and output values of the CRVAL1/2 keywords when the script is run. One can avoid this by setting the optional parameter to anything other than `True`. + +The script keeps track of its executions using a log file called "pure_parallel_wcs_logfile" in the current directory. + + +## Dependencies + +The script mentioned above has external dependencies. It requires the packages astropy, mastquery, numpy, and pysiaf to be installed. The script assumes that an environment with the packages mentioned above has been installed on your machine and that you have loaded into that environment. + +## Notebooks + +### `NIRISS_correct_pure_parallel_WCS.ipynb` + +This notebook illustrates the impact of the use of the `update_parallel_wcs.py` script. diff --git a/notebooks/cross_instrument/update_pure_parallel_wcs/requirements.txt b/notebooks/cross_instrument/update_pure_parallel_wcs/requirements.txt new file mode 100644 index 000000000..e6a0b2391 --- /dev/null +++ b/notebooks/cross_instrument/update_pure_parallel_wcs/requirements.txt @@ -0,0 +1,7 @@ +numpy>=1.25.2 +astropy>=5.3.4 +jwst>=1.15.1 +matplotlib>=3.8.0 +pysiaf>=0.23.3 +astroquery>=0.4.6 +mastquery>=1.7.1 diff --git a/notebooks/cross_instrument/update_pure_parallel_wcs/update_parallel_wcs.py b/notebooks/cross_instrument/update_pure_parallel_wcs/update_parallel_wcs.py new file mode 100644 index 000000000..a573820e7 --- /dev/null +++ b/notebooks/cross_instrument/update_pure_parallel_wcs/update_parallel_wcs.py @@ -0,0 +1,163 @@ +import os +import sys +import astropy.io.fits as pyfits + + +def log_comment(LOGFILE, comment, verbose=False, show_date=False, mode='a'): + + """ + Log a message to a file, optionally including a date tag + """ + import time + + if show_date: + msg = '# ({0})\n'.format(time.ctime()) + else: + msg = '' + + msg += '{0}\n'.format(comment) + + if LOGFILE is not None: + fp = open(LOGFILE, mode) + fp.write(msg) + fp.close() + + if verbose: + print(msg[:-1]) + + return msg + + +def update_pure_parallel_wcs(file, logfile="pure_parallel_wcs_logfile", + fix_vtype='PARALLEL_PURE', verbose=True): + + """ + Update pointing-related keywords of pure parallel exposures using the + pointing info from the FGS (and the prime exposures) from the MAST database + and `pysiaf` + + 1. Find the FGS log from a MAST query that starts before the pure parallel + exposure starts and ends after the exposure ends. + 2. Use the ``ra_v1, dec_v1, pa_v3`` values from the FGS log to set the + pointing attitude with `pysiaf`. + 3. Compute the sky position of the ``CRPIX`` reference pixel of ``file`` + with `pysiaf` and put that position in the ``CRVAL`` keywords. + + Parameters + ---------- + file : str + Filename of a pure-parallel exposure (typically a rate.fits file) + + fix_vtype : str + Only run if ``file[0].header['VISITYPE'] == fix_vtype`` + + verbose : bool + Status messaging + + Returns + ------- + status : None, True + Returns None if some problem is found + + """ + import pysiaf + import mastquery.jwst as jwstquery + + if not os.path.exists(file): + msg = "PureParallelUtils.update_pure_parallel_wcs: " + msg += f" {file} not found" + log_comment(logfile, msg, verbose=verbose) + return None + + with pyfits.open(file) as im: + h0 = im[0].header.copy() + h1 = im[1].header.copy() + if 'VISITYPE' not in im[0].header: + msg = "PureParallelUtils.update_pure_parallel_wcs: " + msg += f" VISITYPE not found in header {file}" + log_comment(logfile, msg, verbose=verbose) + return None + + # Only continue if this is a Pure Parallel exposure + vtype = h0['VISITYPE'] + + if vtype != fix_vtype: + msg = "PureParallelUtils.update_pure_parallel_wcs: " + msg += f" VISITYPE ({vtype}) != {fix_vtype}, skip" + log_comment(logfile, msg, verbose=verbose) + return None + + crval_init = h1['CRVAL1'], h1['CRVAL2'] + + # Get correct pointings from FGS logs. Allow the associated FGS exposure + # to start up to 0.02 day ~ 30 min before the science exposure. + dt = 0.02 + gs = jwstquery.query_guidestar_log( + mjd=(h0['EXPSTART']-dt, h0['EXPEND']+dt), + program=None, + exp_type=['FGS_FINEGUIDE'], + ) + + keep = (gs['expstart'] < h0['EXPSTART']) + keep &= (gs['expend'] > h0['EXPEND']) + + if keep.sum() == 0: + msg = f"PureParallelUtils.update_pure_parallel_wcs: par_file='{file}'" + msg += " couldn't find corresponding exposure in FGS logs" + log_comment(logfile, msg, verbose=verbose) + return None + + gs = gs[keep][0] + pos = (gs['ra_v1'], gs['dec_v1'], gs['pa_v3']) + attmat = pysiaf.utils.rotations.attitude(0.0, 0.0, *pos) + + # And apply the pointing to the parallel aperture and reference pixel + par_aper = pysiaf.Siaf(h0['INSTRUME'])[h0['APERNAME']] + par_aper.set_attitude_matrix(attmat) + + crpix = h1['CRPIX1'], h1['CRPIX2'] + crpix_init = par_aper.sky_to_sci(*crval_init) + + crval_fix = par_aper.sci_to_sky(*crpix) + + msg = ( + f"PureParallelUtils.update_pure_parallel_wcs: File: {file}\n" + f"PureParallelUtils.update_pure_parallel_wcs: FGS: {gs['fileName']}\n" + f"PureParallelUtils.update_pure_parallel_wcs: original crval: " + f"{crval_init[0]:.7f} {crval_init[1]:.7f}\n" + f"PureParallelUtils.update_pure_parallel_wcs: new crval: " + f"{crval_fix[0]:.7f} {crval_fix[1]:.7f}\n" + f"PureParallelUtils.update_pure_parallel_wcs: d(pix): " + f"{crpix[0] - crpix_init[0]:6.3f} {crpix[1] - crpix_init[1]:6.3f}" + ) + + _ = log_comment(logfile, msg, verbose=verbose) + + with pyfits.open(file, mode='update') as im: + im[1].header['RA_V1'] = gs['ra_v1'] + im[1].header['DEC_V1'] = gs['dec_v1'] + im[1].header['PA_V3'] = gs['pa_v3'] + im[1].header['CRVAL1'] = crval_fix[0] + im[1].header['CRVAL2'] = crval_fix[1] + im[1].header['PUREPWCS'] = True, 'WCS updated from query of FGS log file' + im[1].header['PUREPEXP'] = gs['fileName'], 'FGS log file' + + im.flush() + + return True + + +if __name__ == "__main__": + + if len(sys.argv) < 2: + print('Syntax: update_parallel_wcs.py fitsfile ') + print(' where fitsfile is typically a _rate.fits or _rateints.fits file') + print(' and is a boolean to decide whether or not') + print(' to print the results (the default is `True`)') + sys.exit() + fitsfile = str(sys.argv[1]) + verbose = True + if len(sys.argv) > 2: + if str(sys.argv[2]).upper() != 'TRUE': + verbose = False + update_pure_parallel_wcs(fitsfile, verbose=verbose)