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

docs/uvcal_tutorial.rst

@@ -42,6 +42,23 @@ gives the beginning and end of the time range. The local sidereal times follow a
 pattern, UVCal objects should have either an ``lst_array`` or an ``lst_range`` attribute
 set.
 
+Generating calibration solutions typically requires choosing a convention concerning how
+polarized sky emission is mapped to the linear polarizations of the instrument. For
+linear polarizations ``XX`` and ``YY``, the stokes ``I`` sky emission can be mapped to
+``I = (XX + YY)/2`` (the ``avg`` convention) or ``I = XX + YY`` (the ``sum``
+convention). This choice is generally encoded in the sky model to which the visibilities
+are calibrated. Different tools and simulators make different choices, generally following
+a standard choice for the field. In ``pyuvdata`` either of these choices are OK, but the
+choice should be recorded as the ``pol_convention`` parameter in both ``UVCal`` and
+``UVData`` objects. Since the ``pol_convention`` has always (at least implicitly) been
+chosen for calibration solutions, we suggest *always* specifying this parameter on the
+``UVCal`` object (though we do not enforce this, for backwards compatibility reasons).
+Only *calibrated* ``UVData`` objects make sense to have the ``pol_convention`` specified.
+To learn more about this parameter and how ``pyuvdata`` deals with it, please see the
+section below `UVCal: Calibrating UVData`_.
+
+
+
 For most users, the convenience methods for quick data access (see
 `UVCal: Quick data access`_) are the easiest way to get data for particular antennas.
 Those methods take the antenna numbers (i.e. numbers listed in ``telescope.antenna_numbers``)
@@ -260,13 +277,72 @@ UVCal: Calibrating UVData
 Calibration solutions in a :class:`pyuvdata.UVCal` object can be applied to a
 :class:`pyuvdata.UVData` object using the :func:`pyuvdata.utils.uvcalibrate` function.
 
+Generating calibration solutions typically requires choosing a convention concerning how
+polarized sky emission is mapped to the linear polarizations of the instrument. For
+linear polarizations ``XX`` and ``YY``, the stokes ``I`` sky emission can be mapped to
+``I = (XX + YY)/2`` (the ``avg`` convention) or ``I = XX + YY`` (the ``sum``
+convention). This choice is generally encoded in the sky model to which the visibilities
+are calibrated. Different tools and simulators make different choices, generally following
+a standard choice for the field. When calibrating a ``UVData`` object with a ``UVCal``
+object using :func:`pyuvdata.utils.uvcalibrate`, it is *required* to specify this
+convention. At this time, the convention can be specified either on the ``UVCal`` object
+itself, or as a parameter to :func:`pyuvdata.utils.uvcalibrate`. The chosen ``pol_convention``
+will then be applied to and stored on the resulting ``UVData`` object.
+
+There are a few non-trivial combinations of parameters concerning the ``pol_convention``
+that should be noted:
+
+* There are two parameters to :func:`pyuvdata.utils.uvcalibrate` that specify how
+  the convention should be handled: ``uvd_pol_convention`` and ``uvc_pol_convention``,
+  and these act differently depending on whether ``undo`` is True or False. The
+  ``uvc_pol_convention`` is only ever meant to specify what convention the ``UVCal``
+  object actually uses, and is therefore unnecessary if ``UVCal.pol_convention`` is
+  specified (regardless of whether calibrating or uncalibrating). On the other hand,
+  the ``uvd_pol_convention`` specifies the *desired* convention on the resulting
+  ``UVData`` object if calibrating, and otherwise specifies the actual convention on
+  the ``UVData`` object (if uncalibrating, and this convention is not already specified
+  on the object itself).
+* Regardless of the value of ``undo``, the convention that is inferred for the
+  calibration solutions is determined as follows:
+  * If neither ``uvc_pol_convention`` nor ``uvcal.pol_convention`` are specified, a
+    a warning is raised (since the resulting calibrated data is not well-determined),
+    and it is *assumed* that the solutions have the same convention as the ``UVData``
+    (i.e. the desired convention in the case of calibration, or the actual convention
+    in the case of uncalibration). If these are also not specified, no convention
+    corrections are applied, and the result is ambiguous.
+  * If both ``uvc_pol_convention`` and ``uvcal.pol_convention`` are specified and are
+    different, an error is raised.
+* When **calibrating** in :func:`pyuvdata.utils.uvcalibrate` (i.e. ``undo=False``):
+  * If ``uvdata.pol_convention`` is specified, an error is raised, because you are
+    trying to calibrate already-calibrated data.
+  * The convention applied to the resulting ``UVData`` object is inferred in the
+    following precedence: (i) the value of ``uvd_pol_convention``, (ii) whatever is
+    specified as the convention of the ``UVCal`` object (either via ``uvc_pol_convention``
+    or ``uvcal.pol_convention``, see above), (iii) if still unspecified, no convention
+    will be used and a warning will be raised. This was always the behaviour in earlier
+    versions of ``pyuvdata`` (pre-v3).
+* When **un-calibrating** with :func:`pyuvdata.utils.uvcalibrate` (i.e. ``undo=True``):
+  * If both ``uvd_pol_convention`` and ``uvdata.pol_convention`` are defined and
+    are different, an error is raised.
+  * If neither are set, a warning is raised, since the resulting un-calibrated values
+    may not be the same as the original values before calibration (since a different
+    convention could have been used to calibrate originally than is being used to
+    de-calibrate). However, calibration will continue, assuming that the ``UVData``
+    object has the same convention as the ``UVCal`` object used to de-calibrate.
+* It is not supported to have ``pol_convention`` set on ``UVCal``, but *not*
+  ``gain_scale``. A ``pol_convention`` only makes sense in the context of having a
+  scale for the gains.
+* Mis-matching ``uvd_pol_convention`` and ``uvc_pol_convention`` is perfectly fine: any
+  necessary corrections in the calibration will be made to obtain the correct desired
+  convention.
 
 a) Calibration of UVData by UVCal
 *********************************
 .. code-block:: python
 
   >>> # We can calibrate directly using a UVCal object
   >>> import os
+  >>> import numpy as np
   >>> from pyuvdata import UVData, UVCal, utils
   >>> from pyuvdata.data import DATA_PATH
   >>> uvd = UVData.from_file(
@@ -281,7 +357,14 @@ a) Calibration of UVData by UVCal
   >>> uvc.telescope.antenna_names = np.array(
   ...     [name.replace("ant", "HH") for name in uvc.telescope.antenna_names]
   ... )
+  >>> # We should also set the gain_scale and pol_convention, which was not set
+  >>> # in this old file. In old HERA files, like this one, the pol_convention
+  >>> # was implicitly "avg" but in new files it is explicitly "sum"
+  >>> uvc.gain_scale = "Jy"
+  >>> uvc.pol_convention = "avg"
   >>> uvd_calibrated = utils.uvcalibrate(uvd, uvc, inplace=False)
+  >>> print(uvd_calibrated.pol_convention)
+  avg
 
   >>> # We can also un-calibrate using the same UVCal
   >>> uvd_uncalibrated = utils.uvcalibrate(uvd_calibrated, uvc, inplace=False, undo=True)

src/pyuvdata/utils/uvcalibrate.py

@@ -3,24 +3,107 @@
 # 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
 
 from .pol import POL_TO_FEED_DICT, jnum2str, parse_jpolstr, polnum2str, polstr2num
 
 
+def _get_pol_conventions(
+    uvdata,
+    uvcal,
+    undo: bool,
+    uvc_pol_convention: Literal["sum", "avg"] | None,
+    uvd_pol_convention: Literal["sum", "avg"] | None,
+):
+    if uvc_pol_convention is None and uvcal.pol_convention is None:
+        warnings.warn(
+            message=(
+                "pol_convention is not specified on the UVCal object, and "
+                "uvc_pol_convention was not specified. This is deprecated, and will"
+                " be an error in pyuvdata 3.2, since it is impossible to properly "
+                "infer. For now, we are tentatively assuming "
+                "that the UVCal and UVData objects (implicitly) have the same "
+                "convention."
+            ),
+            category=DeprecationWarning,
+            stacklevel=2,
+        )
+        uvc_pol_convention = uvd_pol_convention or uvdata.pol_convention
+    elif uvc_pol_convention is None:
+        uvc_pol_convention = uvcal.pol_convention
+    elif uvc_pol_convention != uvcal.pol_convention:
+        raise ValueError(
+            "uvc_pol_convention is set, and different than uvcal.pol_convention."
+        )
+
+    if undo:
+        if uvd_pol_convention is None and uvdata.pol_convention is None:
+            warnings.warn(
+                message=(
+                    "pol_convention is not specified on the UVData object, and "
+                    "uvd_pol_convention was not specified. This is deprecated, and will"
+                    " be an error in pyuvdata 3.2, since it is impossible to know how "
+                    "to properly un-calibrate. For now, we are tentatively assuming "
+                    "that the UVCal and UVData objects (implicitly) have the same "
+                    "convention."
+                ),
+                category=DeprecationWarning,
+                stacklevel=2,
+            )
+            uvd_pol_convention = uvc_pol_convention
+        elif uvd_pol_convention is None:
+            uvd_pol_convention = uvdata.pol_convention
+        elif uvd_pol_convention != uvdata.pol_convention:
+            raise ValueError(
+                "Both uvd_pol_convention and uvdata.pol_convention were specified with"
+                f"different values: {uvd_pol_convention} and {uvdata.pol_convention}."
+            )
+    else:
+        if uvdata.pol_convention is not None:
+            raise ValueError("You are trying to calibrate already-calibrated data.")
+        if uvd_pol_convention is None:
+            uvd_pol_convention = uvc_pol_convention
+
+        if uvd_pol_convention is None:
+            # Both uvc and uvd have no pol_convention specified
+            warnings.warn(
+                message=(
+                    "Neither uvd_pol_convention nor uvc_pol_convention are 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,
+            )
+
+    if uvd_pol_convention not in ["sum", "avg", None]:
+        raise ValueError(
+            f"uvd_pol_convention must be 'sum' or 'avg'. Got {uvd_pol_convention}"
+        )
+    if uvc_pol_convention not in ["sum", "avg", None]:
+        raise ValueError(
+            f"uvc_pol_convention must be 'sum' or 'avg'. Got {uvc_pol_convention}"
+        )
+
+    return uvc_pol_convention, uvd_pol_convention
+
+
 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,
+    uvc_pol_convention: Literal["sum", "avg"] | None = None,
+    uvd_pol_convention: Literal["sum", "avg"] | None = None,
 ):
     """
     Calibrate a UVData object with a UVCal object.
@@ -61,6 +144,16 @@ 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 for how instrumental polarizations (e.g. XX and YY) are
+        converted to Stokes parameters. Options are 'sum' and 'avg', corresponding
+        to I=XX+YY and I=(XX+YY)/2 (for linear instrumental polarizations)
+        respectively. If None, the convention is determined from the UVCal object.
+        If `pol_convention` is not specified and is not set on the UVCal object,
+        a deprecation warning is raised (will be an error in the future). 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 +172,25 @@ def uvcalibrate(
             "calibrations"
         )
 
+    if uvcal.gain_scale is None:
+        warnings.warn(
+            "gain_scale is not set, so there is no way to know what the resulting units"
+            " are. For now, we assume that `gain_scale` matches whatever is on the "
+            "UVData object (i.e. we do not change its units). Furthermore, all "
+            "corrections concerning the pol_convention will be ignored.",
+            category=UserWarning,
+            stacklevel=2,
+        )
+    elif undo and uvcal.gain_scale != uvdata.vis_units:
+        raise ValueError(
+            "Cannot undo calibration if gain_scale is not the same as the units on "
+            "the UVData object."
+        )
+
+    uvc_pol_convention, uvd_pol_convention = _get_pol_conventions(
+        uvdata, uvcal, undo, uvc_pol_convention, uvd_pol_convention
+    )
+
     if not inplace:
         uvdata = uvdata.copy()
 
@@ -417,5 +529,29 @@ def uvcalibrate(
         if uvcal_use.gain_scale is not None:
             uvdata.vis_units = uvcal_use.gain_scale
 
+    # Set pol convention properly
+    if uvcal.gain_scale is not None:
+        _apply_pol_convention_corrections(
+            uvdata, undo, uvd_pol_convention, uvc_pol_convention
+        )
+
     if not inplace:
         return uvdata
+
+
+def _apply_pol_convention_corrections(
+    uvdata, undo, uvd_pol_convention, uvc_pol_convention
+):
+    if uvd_pol_convention != uvc_pol_convention:
+        correction = np.ones(uvdata.Npols)
+        correction[uvdata.polarization_array > 0] *= 2
+        correction[uvdata.polarization_array < 0] /= 2
+
+        if (undo and uvd_pol_convention == "sum") or (
+            not undo and uvd_pol_convention == "avg"
+        ):
+            correction = 1 / correction
+
+        uvdata.data_array *= correction
+
+    uvdata.pol_convention = None if undo else uvd_pol_convention

src/pyuvdata/uvcal/calfits.py

@@ -207,6 +207,8 @@ def write_calfits(
             prihdr["DIFFUSE"] = self.diffuse_model
         if self.gain_scale is not None:
             prihdr["GNSCALE"] = self.gain_scale
+        if self.pol_convention is not None:
+            prihdr["POLCONV"] = self.pol_convention
         if self.Ntimes > 1:
             prihdr["INTTIME"] = median_int_time
         else:
@@ -599,13 +601,16 @@ def read_calfits(
 
             self.gain_convention = hdr.pop("GNCONVEN")
             self.gain_scale = hdr.pop("GNSCALE", None)
+            self.pol_convention = hdr.pop("POLCONV", None)
             self.cal_type = hdr.pop("CALTYPE")
 
             # old files might have a freq range for gain types but we don't want them
             if self.cal_type == "delay":
                 self.freq_range = np.array(
                     [list(map(float, hdr.pop("FRQRANGE").split(",")))]
                 )
+            else:
+                hdr.pop("FRQRANGE", None)
 
             self.cal_style = hdr.pop("CALSTYLE")
             if self.cal_style == "sky":

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/fhd_cal.py

@@ -234,6 +234,8 @@ def read_fhd_cal(
 
         self._set_sky()
         self.gain_convention = "divide"
+        self.gain_scale = "Jy"
+        self.pol_convetions = "sum"
         self._set_gain()
 
         # currently don't have branch info. may change in future.

src/pyuvdata/uvcal/ms_cal.py

@@ -234,6 +234,7 @@ def read_ms_cal(
 
         self.sky_catalog = main_keywords.get("pyuvdata_sky_catalog", None)
         self.gain_scale = main_keywords.get("pyuvdata_gain_scale", None)
+        self.pol_convention = main_keywords.get("pyuvdata_polconv", None)
         self.observer = main_keywords.get("pyuvdata_observer", None)
         if "pyuvdata_cal_style" in main_keywords:
             self.cal_style = main_keywords["pyuvdata_cal_style"]
@@ -491,6 +492,8 @@ def write_ms_cal(self, filename, clobber=False):
 
             if self.gain_scale is not None:
                 ms.putkeyword("pyuvdata_gain_scale", self.gain_scale)
+            if self.pol_convention is not None:
+                ms.putkeyword("pyuvdata_polconv", self.pol_convention)
 
             if self.observer is not None:
                 ms.putkeyword("pyuvdata_observer", self.observer)