Add polangle

arrakis/imager.py

@@ -17,6 +17,7 @@
 from astropy.io import fits
 from astropy.stats import mad_std
 from astropy.table import Table
+import astropy.units as u
 from fitscube import combine_fits
 from fixms.fix_ms_corrs import fix_ms_corrs
 from fixms.fix_ms_dir import fix_ms_dir
@@ -32,6 +33,7 @@
     field_idx_from_ms,
     field_name_from_ms,
     wsclean,
+    get_pol_axis,
 )
 from arrakis.utils.pipeline import logo_str, workdir_arg_parser
 
@@ -377,7 +379,11 @@ def image_beam(
 
 @task(name="Make Cube")
 def make_cube(
-    pol: str, image_set: ImageSet, common_beam_pkl: Path, aux_mode: Optional[str] = None
+    pol: str,
+    image_set: ImageSet,
+    common_beam_pkl: Path,
+    pol_angle: u.Quantity,
+    aux_mode: Optional[str] = None,
 ) -> Tuple[Path, Path]:
     """Make a cube from the images"""
     logger = get_run_logger()
@@ -392,6 +398,12 @@ def make_cube(
     new_header = hdu_list[0].header
     data_cube = hdu_list[0].data
 
+    # Add pol angle to header
+    new_header["INSTRUMENT_RECEPTOR_ANGLE"] = (
+        pol_angle.to(u.deg).value,
+        "Orig. pol. axis rotation angle in degrees",
+    )
+
     tmp_header = new_header.copy()
     # Need to swap NAXIS 3 and 4 to make LINMOS happy - booo
     for a, b in ((3, 4), (4, 3)):
@@ -726,8 +738,10 @@ def main(
             ms_fix = fix_ms_askap_corrs(
                 ms=ms_fix, data_column="DATA", corrected_data_column=data_column
             )
+            pol_angle = get_pol_axis(ms_fix, col="INSTRUMENT_RECEPTOR_ANGLE")
         else:
             ms_fix = ms
+            pol_angle = get_pol_axis(ms_fix, col="RECEPTOR_ANGLE")
         # Image with wsclean
         image_set = image_beam.submit(
             ms=ms_fix,
@@ -790,6 +804,7 @@ def main(
                     pol=pol,
                     image_set=sm_image_set,
                     common_beam_pkl=common_beam_pkl,
+                    pol_angle=pol_angle,
                     aux_mode=aux_mode,
                     wait_for=[sm_image_set],
                 )

arrakis/linmos.py

@@ -12,6 +12,9 @@
 from typing import Dict, List, Optional, Tuple
 from typing import NamedTuple as Struct
 
+from astropy.io import fits
+import astropy.units as u
+import numpy as np
 import pandas as pd
 import pymongo
 from astropy.utils.exceptions import AstropyWarning
@@ -160,6 +163,16 @@ def genparset(
     """
     logger.setLevel(logging.INFO)
 
+    pol_angles_list: List[float] = []
+    for im in image_paths.images:
+        _pol_angle: float = fits.getheader(im)["INSTRUMENT_RECEPTOR_ANGLE"]
+        pol_angles_list.append(_pol_angle)
+    pol_angles: u.Quantity = pol_angles_list * u.deg
+
+    alpha: u.Quantity = pol_angles[0]
+
+    assert np.allclose(pol_angles, alpha), "Polarisation angles are not the same!"
+
     image_string = f"[{','.join([im.resolve().with_suffix('').as_posix() for im in image_paths.images])}]"
     weight_string = f"[{','.join([im.resolve().with_suffix('').as_posix() for im in image_paths.weights])}]"
 
@@ -188,6 +201,7 @@ def genparset(
         parset += f"""
 linmos.primarybeam      = ASKAP_PB
 linmos.primarybeam.ASKAP_PB.image = {holofile.resolve().as_posix()}
+linmos.primarybeamASKAP_PB.alpha = {alpha.to(u.rad).value}
 linmos.removeleakage    = true
 """
     else:

arrakis/utils/msutils.py

@@ -3,9 +3,11 @@
 
 import copy
 import warnings
+from pathlib import Path
 from typing import Optional
 
 from astropy.utils.exceptions import AstropyWarning
+import astropy.units as u
 from casacore.tables import table
 from spectral_cube.utils import SpectralCubeWarning
 
@@ -15,6 +17,46 @@
 warnings.simplefilter("ignore", category=AstropyWarning)
 
 
+def get_pol_axis(
+    ms: Path, feed_idx: Optional[int] = None, col: str = "RECEPTOR_ANGLE"
+) -> u.Quantity:
+    """Get the polarization axis from the ASKAP MS. Checks are performed
+    to ensure this polarisation axis angle is constant throughout the observation.
+
+
+    Args:
+        ms (Path): The path to the measurement set that will be inspected
+        feed_idx (Optional[int], optional): Specify the entery in the FEED
+        table of `ms` to return. This might be required when a subset of a
+        measurement set has been extracted from an observation with a varying
+        orientation.
+        col (str, optional): The column to extract the polarization angle from.
+
+    Returns:
+        astropy.units.Quantity: The rotation of the PAF throughout the observing.
+    """
+    _known_cols = ("RECEPTOR_ANGLE", "INSTRUMENT_RECEPTOR_ANGLE")
+    if col not in _known_cols:
+        raise ValueError(f"Unknown column {col=}, please use one of {_known_cols}")
+    with table((ms / "FEED").as_posix(), readonly=True, ack=False) as tf:
+        ms_feed = tf.getcol(col) * u.rad
+        # PAF is at 45deg to feeds
+        # 45 - feed_angle = pol_angle
+        pol_axes = -(ms_feed - 45.0 * u.deg)
+
+    if feed_idx is None:
+        assert (ms_feed[:, 0] == ms_feed[0, 0]).all() & (
+            ms_feed[:, 1] == ms_feed[0, 1]
+        ).all(), f"The {col} changes with time, please check the MS"
+
+        pol_ang = pol_axes[0, 0].to(u.deg)
+    else:
+        logger.debug(f"Extracting the third-axis orientation for {feed_idx=}")
+        pol_ang = pol_axes[feed_idx, 0].to(u.deg)
+
+    return pol_ang
+
+
 def beam_from_ms(ms: str) -> int:
     """Work out which beam is in this MS"""
     with table(ms, readonly=True, ack=False) as t:

pyproject.toml

@@ -63,7 +63,7 @@ RMTable = ">=1.2.1"
 RM-Tools = ">=1.4.2"
 PolSpectra = ">=1.1.0"
 setuptools = "*"
-fixms = ">=0.2"
+fixms = ">=0.2.5"
 fitscube = ">=0.3"
 psycopg2-binary = "*"
 sqlalchemy = "*"