How to get RMs from FITS images

Given appropriate metadata, it is possible to obtain ionospheric RMs from FITS images. Further, Van Eck (2021) derived a method to obtain a time-integrated correction suitable for applying to image data. There is also Python implementation `FRion <https://github.com/CIRADA-Tools/FRion>`__, which we have partly re-implemented here in spinifex for tight integration with the underlying ionospheric modelling routines.

Van Eck (2021) Equation 3 provides the time-integrated effect on linear polarisation as:

\[\tilde\Theta(\lambda^2) = \frac{1}{t_1 - t_0} \int_{t_0}^{t_1} e^{2i\lambda^2\phi(t)}dt.\]

Where \(t_0\) and \(t_1\) are the start and end times, respecitively, \(\lambda^2\) are the observed wavelengths-squared, and \(\phi\) is the time-dependent Faraday rotation. \(\tilde\Theta\) is a complex quantity whose amplitude can be interpreted as the depolarisation caused by the ionosphere, and the phase as the change of polarisation angle.

For this example, we’ll generate a mock FITS file using astropy.

from __future__ import annotations

from pathlib import Path
from pprint import pprint

import matplotlib.pyplot as plt
import numpy as np
from astropy import units as u
from astropy.io import fits
from astropy.visualization import quantity_support, time_support
from spinifex.image_tools import (
    get_freq_from_fits,
    get_integrated_rm_from_fits,
    get_rm_from_fits,
)

_ = quantity_support()
_ = time_support()

example_fits_file = Path("example.fits")

header = fits.Header()
header["NAXIS"] = 4
header["NAXIS1"] = 10
header["NAXIS2"] = 10
header["NAXIS3"] = 6
header["NAXIS4"] = 1
header["CTYPE1"] = "RA---SIN"
header["CRVAL1"] = 0
header["CRPIX1"] = 5
header["CDELT1"] = -(1 / 3600)
header["CUNIT1"] = "deg"
header["CTYPE2"] = "DEC--SIN"
header["CRVAL2"] = 0
header["CRPIX2"] = 5
header["CDELT2"] = 1 / 3600
header["CUNIT2"] = "deg"
header["CTYPE3"] = "FREQ"
header["CRVAL3"] = 1.4e9
header["CRPIX3"] = 1
header["CDELT3"] = 1e6
header["CUNIT3"] = "Hz"
header["CTYPE4"] = "STOKES"
header["CRVAL4"] = 1
header["CRPIX4"] = 1
header["CDELT4"] = 1
header["CUNIT4"] = ""
header["BUNIT"] = "Jy/beam"
header["DATE-OBS"] = (
    "2019-04-25T12:45:52.893302"  # pulled from a random SPICE-RACS file
)
header["DURATION"] = 900
header["TELESCOP"] = "ASKAP"
data = np.ones((1, 6, 10, 10), dtype=np.float32)
hdu = fits.PrimaryHDU(data, header)

hdu.writeto(example_fits_file, overwrite=True)

# Get the time-dependent RM - this returns the same object as the other `get_rm` functions
rm = get_rm_from_fits(
    fits_path=example_fits_file,
    timestep=1
    * u.min,  # timestep is optional, defaults to 15 minutes - but that's our total duration here!
)


fig, ax = plt.subplots()
ax.errorbar(
    rm.times.datetime,
    rm.rm * u.rad / u.m**2,
    yerr=rm.rm_error * u.rad / u.m**2,
    fmt=".",
)
ax.set(
    xlabel="Time",
    ylabel=f"RM / {u.rad / u.m**2:latex_inline}",
    title="Time-dependent RM",
)

INFO fits_tools - get_metadata_from_fits: Getting telescope location from `TELESCOP` key
INFO models - parse_iono_kwargs: Using ionospheric model <function get_density_ionex_single_layer at 0x7430ee7cd440> with options server=<Servers.CHAPMAN: 'chapman'> prefix='uqr' url_stem=None time_resolution=None solution='final' output_directory=None correct_uqrg_rms=True height=<Quantity 350. km> remove_midnight_jumps=True
INFO get_ipp - _get_ipp_simple: Calculating ionospheric piercepoints
INFO get_rm - _get_rm: Calculating rotation measure
INFO download - download_file_http: Downloading from http://chapman.upc.es/tomion/rapid/2019/115_190425.15min/uqrg1150.19i.Z
INFO download - download_file_http: Downloading from http://chapman.upc.es/tomion/rapid/2019/116_190426.15min/uqrg1160.19i.Z

[Text(0.5, 0, 'Time'),
 Text(0, 0.5, 'RM / $\\mathrm{rad\\,m^{-2}}$'),
 Text(0.5, 1.0, 'Time-dependent RM')]

# Get the time-integrated product
integrated_rm = get_integrated_rm_from_fits(
    fits_path=example_fits_file,
    timestep=1
    * u.min,  # timestep is optional, defaults to 15 minutes - but that's our total duration here!
)
# We'll use our helper function to get the frequency axis
frequencies = get_freq_from_fits(example_fits_file)
fig, ax = plt.subplots()
ax2 = ax.twinx()
ax.plot(frequencies.to(u.GHz), np.abs(integrated_rm.theta), label=r"$|\Theta|$")
ax.set(
    ylim=(0, 1.1),
    ylabel="Fractional polarisation",
    xlabel=f"Frequency / {u.GHz:latex_inline}",
    title="Time-integrated polarisation effect",
)
ax2.plot(
    frequencies.to(u.GHz),
    (np.angle(integrated_rm.theta) * u.rad).to(u.deg),
    color="tab:orange",
    label=r"$\arg(\Theta)$",
)
ax2.set(ylim=(-90, 90), ylabel="Polarisation angle")
fig.legend()

INFO fits_tools - get_metadata_from_fits: Getting telescope location from `TELESCOP` key
INFO models - parse_iono_kwargs: Using ionospheric model <function get_density_ionex_single_layer at 0x7430ee7cd440> with options server=<Servers.CHAPMAN: 'chapman'> prefix='uqr' url_stem=None time_resolution=None solution='final' output_directory=None correct_uqrg_rms=True height=<Quantity 350. km> remove_midnight_jumps=True
INFO get_ipp - _get_ipp_simple: Calculating ionospheric piercepoints
INFO get_rm - _get_rm: Calculating rotation measure
INFO download - download_or_copy_url: File /home/docs/checkouts/readthedocs.org/user_builds/spinifex/checkouts/latest/docs/source/examples/ionex_files/uqrg1150.19i.Z already exists. Skipping download.
INFO download - download_or_copy_url: File /home/docs/checkouts/readthedocs.org/user_builds/spinifex/checkouts/latest/docs/source/examples/ionex_files/uqrg1160.19i.Z already exists. Skipping download.
WARNING fits_tools - get_freq_from_fits: Could not find `SPECSYS` in header, assuming `TOPOCENT`
WARNING fits_tools - get_freq_from_fits: Could not find `SPECSYS` in header, assuming `TOPOCENT`

<matplotlib.legend.Legend at 0x7430de6916a0>

# The integrated RM object is a namedtuple with the following fields
# These are mostly the same as the RM object, but averaged over time
pprint(integrated_rm._asdict())

{'azimuth': np.float64(-145.46225393250046),
 'b_parallel': np.float64(19656.944490267193),
 'electron_density': np.float64(11.586633172722578),
 'elevation': np.float64(-58.69630540677612),
 'height': np.float64(450.0153044178698),
 'theta': array([0.99804308-0.06252483j, 0.99804866-0.06243572j,
       0.99805422-0.0623468j , 0.99805976-0.06225807j,
       0.99806528-0.06216953j, 0.99807078-0.06208118j]),
 'time': <Time object: scale='utc' format='mjd' value=58598.536723302124>}