Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

adding update_pure_parallel_wcs #241

Merged
merged 16 commits into from
Sep 17, 2024
Merged

adding update_pure_parallel_wcs #241

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
16 commits
Select commit Hold shift + click to select a range
ae89489
adding update_pure_parallel_wcs
goudfroo Sep 3, 2024
f4aaf85
adding one sentence to NIRISS_correct_pure_parallel_WCS.ipynb
goudfroo Sep 3, 2024
2482e94
adding update_parallel_wcs.py
goudfroo Sep 4, 2024
59f36e6
adding missing whitespace in code
goudfroo Sep 4, 2024
ee8e069
some indents were not multiples of 4
goudfroo Sep 4, 2024
8cef9df
correcting one more indent
goudfroo Sep 4, 2024
542e268
adding one black line (gosh...)
goudfroo Sep 5, 2024
42a9660
few pep8 errors in script addressed
goudfroo Sep 5, 2024
af4aedd
few more pep8 errors in script addressed
goudfroo Sep 5, 2024
d3f69b4
implementing suggestions of reviewer
goudfroo Sep 17, 2024
4f283f6
adding mastquery to requirements
goudfroo Sep 17, 2024
33ce565
adding mastquery to requirements
goudfroo Sep 17, 2024
bbd0666
clearing output cells
goudfroo Sep 17, 2024
aa2cf8e
correcting typo
goudfroo Sep 17, 2024
7522b24
trying to clear one import
goudfroo Sep 17, 2024
fc79c4b
Add new nb to the toc
haticekaratay Sep 17, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions _toc.yml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
374 changes: 374 additions & 0 deletions notebooks/cross_instrument/update_pure_parallel_wcs/NIRISS_correct_pure_parallel_WCS.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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 <verbosity>\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 <verbosity> 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
}
Loading
Loading