Pol convention

CHANGELOG.md

@@ -4,6 +4,9 @@ All notable changes to this project will be documented in this file.
 ## [Unreleased]
 
 ### Added
+- New UVParameter `pol_convention` on `UVData` and `UVCal`. This specifies the convention
+assumed for converting linear to stokes polarizations -- either "sum" or "avg". Also
+added to `uvcalibrate` to apply from the `UVCal` to the `UVData`.
 - New `ignore_telescope_param_update_warnings_for` function that globally ignores
 warnings for specific telescopes.
 - New `.telescope` cached property on the `FastUVH5Meta` object that auto-creates a

src/pyuvdata/utils/uvcalibrate.py

@@ -3,6 +3,7 @@
 # Licensed under the 2-clause BSD License
 """Code to apply calibration solutions to visibility data."""
 import warnings
+from typing import Literal
 
 import numpy as np
 
@@ -13,14 +14,15 @@ def uvcalibrate(
     uvdata,
     uvcal,
     *,
-    inplace=True,
-    prop_flags=True,
-    d_term_cal=False,
-    flip_gain_conj=False,
-    delay_convention="minus",
-    undo=False,
-    time_check=True,
-    ant_check=True,
+    inplace: bool = True,
+    prop_flags: bool = True,
+    d_term_cal: bool = False,
+    flip_gain_conj: bool = False,
+    delay_convention: Literal["minus", "plus"] = "minus",
+    undo: bool = False,
+    time_check: bool = True,
+    ant_check: bool = True,
+    pol_convention: Literal["sum", "avg"] | None = None,
 ):
     """
     Calibrate a UVData object with a UVCal object.
@@ -61,6 +63,14 @@ def uvcalibrate(
         object have calibration solutions in the UVCal object. If this option is
         set to False, uvcalibrate will proceed without erroring and data for
         antennas without calibrations will be flagged.
+    pol_convention : str, {"sum", "avg"}, optional
+        The convention used for combining linear polarizations (e.g. XX and YY) into
+        pseudo-Stokes parameters (e.g. I, Q, U, V). Options 'sum' and 'avg' correspond
+        to I=XX+YY and I=(XX+YY)/2 respectively. If None, the convention is determined
+        from the UVCal object. If `pol_convention` is not specified on the UVCal object
+        and is not specified, an error is raised. If the convention specified here is
+        different from that in the UVCal object, the one specified here will take
+        precedence, and the necessary conversions will be made.
 
     Returns
     -------
@@ -79,6 +89,36 @@ def uvcalibrate(
             "calibrations"
         )
 
+    if pol_convention is None:
+        pol_convention = uvdata.pol_convention if undo else uvcal.pol_convention
+
+    if pol_convention is None:
+        if not undo:
+            warnings.warn(
+                message=(
+                    "pol_convention is not specified, so the resulting UVData object "
+                    "will have ambiguous convention. This is deprecated and will be an "
+                    "error in pyuvdata v3.2."
+                ),
+                category=DeprecationWarning,
+                stacklevel=2,
+            )
+        else:
+            warnings.warn(
+                "The pol_convention is not specified on the UVData, so it is impossible"
+                " to know if the correct convention is being applied to undo the "
+                "calibration. This behaviour is deprecated, and for now will result in "
+                "the same convention being assumed for the UVData and UVCal objects. "
+                "Please set the pol_convention attribute on the UVData object.",
+                category=DeprecationWarning,
+                stacklevel=2,
+            )
+
+    elif pol_convention not in ["sum", "avg"]:
+        raise ValueError(
+            f"pol_convention must be either 'sum' or 'avg'. Got {pol_convention}"
+        )
+
     if not inplace:
         uvdata = uvdata.copy()
 
@@ -417,5 +457,43 @@ def uvcalibrate(
         if uvcal_use.gain_scale is not None:
             uvdata.vis_units = uvcal_use.gain_scale
 
+    # Set pol convention properly
+    if uvcal_use.gain_scale is None:
+        # ignore all pol_convention stuff, because who knows???
+        warnings.warn(
+            "gain_scale is not set, so resulting uvdata will not have correct units"
+        )
+    else:
+        if not undo:
+            uvdata.pol_convention = pol_convention
+            if pol_convention != uvcal_use.pol_convention:
+                for i, pol in enumerate(uvdata.polarization_array):
+                    if pol > 0:  # pseudo-stokes
+                        if pol_convention == "sum":
+                            uvdata.data_array[..., i] *= 2
+                        else:
+                            uvdata.data_array[..., i] /= 2
+                    else:  # linear pols
+                        if pol_convention == "sum":
+                            uvdata.data_array[..., i] /= 2
+                        else:
+                            uvdata.data_array[..., i] *= 2
+        else:
+            uvdata.pol_convention = None
+            # here the pol_convention should be what is being assumed/specified
+            # on the uvdata
+            if pol_convention != uvcal_use.pol_convention:
+                for i, pol in enumerate(uvdata.polarization_array):
+                    if pol > 0:  # pseudo-stokes
+                        if pol_convention == "sum":
+                            uvdata.data_array[..., i] /= 2
+                        else:
+                            uvdata.data_array[..., i] *= 2
+                    else:  # linear pols
+                        if pol_convention == "sum":
+                            uvdata.data_array[..., i] *= 2
+                        else:
+                            uvdata.data_array[..., i] /= 2
+
     if not inplace:
         return uvdata

src/pyuvdata/uvcal/calh5.py

@@ -72,6 +72,7 @@ class FastCalH5Meta(hdf5_utils.HDF5Meta):
             "sky_catalog",
             "instrument",
             "version",
+            "pol_convention",
         }
     )
 
@@ -253,6 +254,7 @@ def _read_header(
             "ref_antenna_array",
             "scan_number_array",
             "sky_catalog",
+            "pol_convention",
         ]
 
         for attr in required_parameters:
@@ -773,6 +775,9 @@ def _write_header(self, header):
             header["phase_center_id_array"] = self.phase_center_id_array
         if self.Nphase is not None:
             header["Nphase"] = self.Nphase
+        if self.pol_convention is not None:
+            header["pol_convention"] = np.bytes_(self.pol_convention)
+
         if self.phase_center_catalog is not None:
             pc_group = header.create_group("phase_center_catalog")
             for pc, pc_dict in self.phase_center_catalog.items():

src/pyuvdata/uvcal/uvcal.py

@@ -597,6 +597,22 @@ def __init__(self):
             required=False,
         )
 
+        desc = (
+            "The convention used for combining linear polarizations (e.g. XX and YY) "
+            "into pseudo-Stokes parameters (e.g. I, Q, U, V). Options are 'sum' and "
+            "'avg', corresponding to I=XX+YY and I=(XX+YY)/2 respectively. This "
+            "parameter is not required, for backwards-compatibility reasons, but is "
+            "highly recommended."
+        )
+        self._pol_convention = uvp.UVParameter(
+            "pol_convention",
+            required=False,
+            description=desc,
+            form="str",
+            spoof_val="avg",
+            acceptable_vals=["sum", "avg"],
+        )
+
         super(UVCal, self).__init__()
 
         # Assign attributes to UVParameters after initialization, since UVBase.__init__

src/pyuvdata/uvdata/uvdata.py

@@ -583,6 +583,21 @@ def __init__(self):
             "filename", required=False, description=desc, expected_type=str
         )
 
+        desc = (
+            "The convention used for combining linear polarizations (e.g. XX and YY) "
+            "into pseudo-Stokes parameters (e.g. I, Q, U, V). Options are 'sum' and "
+            "'avg', corresponding to I=XX+YY and I=(XX+YY)/2 respectively. This "
+            "parameter is not required, and only makes sense for calibrated data."
+        )
+        self._pol_convention = uvp.UVParameter(
+            "pol_convention",
+            required=False,
+            description=desc,
+            form="str",
+            spoof_val="avg",
+            acceptable_vals=["sum", "avg"],
+        )
+
         self.__antpair2ind_cache = {}
         self.__key2ind_cache = {}
 
@@ -2199,6 +2214,13 @@ def check(
         )
         logger.debug("... Done UVBase Check")
 
+        # Check consistency between pol_convention and units of data
+        if self.vis_units == "uncalib" and self.pol_convention is not None:
+            raise ValueError(
+                "pol_convention is set but the data is uncalibrated. This "
+                "is not allowed."
+            )
+
         # then run telescope object check
         self.telescope.check(
             check_extra=check_extra, run_check_acceptability=run_check_acceptability

src/pyuvdata/uvdata/uvh5.py

@@ -109,6 +109,7 @@ class FastUVH5Meta(hdf5_utils.HDF5Meta):
             "eq_coeffs_convention",
             "phase_center_frame",
             "version",
+            "pol_convention",
         }
     )
 
@@ -606,6 +607,7 @@ def _read_header_with_fast_meta(
             "flex_spw_id_array",
             "flex_spw_polarization_array",
             "extra_keywords",
+            "pol_convention",
         ]:
             try:
                 setattr(self, attr, getattr(obj, attr))
@@ -1213,6 +1215,8 @@ def _write_header(self, header):
             header["blts_are_rectangular"] = self.blts_are_rectangular
         if self.time_axis_faster_than_bls is not None:
             header["time_axis_faster_than_bls"] = self.time_axis_faster_than_bls
+        if self.pol_convention is not None:
+            header["pol_convention"] = np.bytes_(self.pol_convention)
 
         # write out extra keywords if it exists and has elements
         if self.extra_keywords:

tests/conftest.py

@@ -58,6 +58,9 @@ def uvcalibrate_init_data_main():
             os.path.join(DATA_PATH, "zen.2458098.45361.HH.omni.calfits_downselected")
         )
 
+    uvcal.pol_convention = "avg"
+    uvcal.gain_scale = "Jy"
+
     yield uvdata, uvcal
 
 

tests/utils/test_uvcalibrate.py

@@ -28,6 +28,7 @@ def test_uvcalibrate_apply_gains_oldfiles(uvcalibrate_uvdata_oldfiles):
     # downselect to match each other in shape (but not in actual values!)
     uvd.select(frequencies=uvd.freq_array[:10])
     uvc.select(times=uvc.time_array[:3])
+    uvc.pol_convention = "avg"
 
     with pytest.raises(
         ValueError,
@@ -78,6 +79,8 @@ def test_uvcalibrate_delay_oldfiles(uvcalibrate_uvdata_oldfiles):
     # downselect to match
     uvc.select(times=uvc.time_array[3])
     uvc.gain_convention = "multiply"
+    uvc.pol_convention = "avg"
+    uvc.gain_scale = "Jy"
 
     freq_array_use = np.squeeze(uvd.freq_array)
     chan_with_use = uvd.channel_width
@@ -133,12 +136,30 @@ def test_uvcalibrate(uvcalibrate_data, flip_gain_conj, gain_convention, time_ran
     uvc.gain_convention = gain_convention
 
     if gain_convention == "divide":
-        assert uvc.gain_scale is None
+        # set the gain_scale to None to test handling
+        uvc.gain_scale = None
+        cal_warn_msg = (
+            "gain_scale is not set, so resulting uvdata will not have correct " "units"
+        )
+        cal_warn_type = UserWarning
+        undo_warn_msg = [
+            "The pol_convention is not specified on the UVData, so it is "
+            "impossible to know if the correct convention is being applied to "
+            "undo the calibration. This behaviour is deprecated, and for now "
+            "will result in the same convention being assumed for the UVData and "
+            "UVCal objects. Please set the pol_convention attribute on the UVData "
+            "object.",
+            "gain_scale is not set, so resulting uvdata will not have correct units",
+        ]
+        undo_warn_type = [DeprecationWarning, UserWarning]
     else:
-        # set the gain_scale to "Jy" to test that vis units are set properly
-        uvc.gain_scale = "Jy"
+        cal_warn_msg = ""
+        cal_warn_type = None
+        undo_warn_msg = ""
+        undo_warn_type = None
 
-    uvdcal = uvcalibrate(uvd, uvc, inplace=False, flip_gain_conj=flip_gain_conj)
+    with check_warnings(cal_warn_type, match=cal_warn_msg):
+        uvdcal = uvcalibrate(uvd, uvc, inplace=False, flip_gain_conj=flip_gain_conj)
     if gain_convention == "divide":
         assert uvdcal.vis_units == "uncalib"
     else:
@@ -171,15 +192,16 @@ def test_uvcalibrate(uvcalibrate_data, flip_gain_conj, gain_convention, time_ran
         )
 
     # test undo
-    uvdcal = uvcalibrate(
-        uvdcal,
-        uvc,
-        prop_flags=True,
-        ant_check=False,
-        inplace=False,
-        undo=True,
-        flip_gain_conj=flip_gain_conj,
-    )
+    with check_warnings(undo_warn_type, match=undo_warn_msg):
+        uvdcal = uvcalibrate(
+            uvdcal,
+            uvc,
+            prop_flags=True,
+            ant_check=False,
+            inplace=False,
+            undo=True,
+            flip_gain_conj=flip_gain_conj,
+        )
 
     np.testing.assert_array_almost_equal(uvd.get_data(key), uvdcal.get_data(key))
     assert uvdcal.vis_units == "uncalib"
@@ -250,6 +272,7 @@ def test_uvcalibrate_flag_propagation(uvcalibrate_data):
 
     assert exp_err == str(errinfo.value)
 
+    uvc_sub.gain_scale = "Jy"
     with pytest.warns(UserWarning) as warninfo:
         uvdcal = uvcalibrate(
             uvd, uvc_sub, prop_flags=True, ant_check=False, inplace=False
@@ -268,6 +291,7 @@ def test_uvcalibrate_flag_propagation(uvcalibrate_data):
 @pytest.mark.filterwarnings("ignore:Cannot preserve total_quality_array")
 def test_uvcalibrate_flag_propagation_name_mismatch(uvcalibrate_init_data):
     uvd, uvc = uvcalibrate_init_data
+    uvc.gain_scale = "Jy"
 
     # test flag propagation
     uvc.flag_array[0] = True
@@ -320,6 +344,7 @@ def test_uvcalibrate_extra_cal_antennas(uvcalibrate_data):
 
 def test_uvcalibrate_antenna_names_mismatch(uvcalibrate_init_data):
     uvd, uvc = uvcalibrate_init_data
+    uvc.gain_scale = "Jy"
 
     with pytest.raises(
         ValueError,
@@ -346,7 +371,7 @@ def test_uvcalibrate_antenna_names_mismatch(uvcalibrate_init_data):
 @pytest.mark.parametrize("time_range", [True, False])
 def test_uvcalibrate_time_mismatch(uvcalibrate_data, time_range):
     uvd, uvc = uvcalibrate_data
-
+    uvc.gain_scale = "Jy"
     if time_range:
         tstarts = uvc.time_array - uvc.integration_time / (86400 * 2)
         tends = uvc.time_array + uvc.integration_time / (86400 * 2)