casa / casa6 / 91bb1c7b4aa

Commits

Federico Montesino Pouzols authored 91bb1c7b4aa 20 Jan 2025

don't convert to J2000, write phasecenter vals using same frame as in the PHASE_DIR metadata, CAS-13616

casatasks/src/private/task_phaseshift.py

Modified

--
+from typing import Optional, Tuple, Union
-+
+from typing import Dict, Optional, Tuple, Union
-2
-3
+import numpy as np
-4
+from .mstools import write_history
-5
+from casatools import table, ms, mstransformer
-6
+from casatools import measures as me
-7
+from casatasks import casalog
-8
+from .parallel.parallel_data_helper import ParallelDataHelper
-9
-10
-11
+def phaseshift(
-133
+    try:
-134
+        tblocal.open(vis)
-135
+        colnames = tblocal.colnames()
-136
+    finally:
-137
+        tblocal.done()
-138
+    return colnames
-139
-140
-141
+def _update_field_subtable(outputvis: str, field: str, phasecenter: Union[str, dict]):
-142
+    """Update MS/FIELD subtable with shifted center(s)."""
-+
-144
+    try:
-145
+        tblocal = table()
-146
+        # modify FIELD table
-147
+        tblocal.open(outputvis + "/FIELD", nomodify=False)
-148
+        pcol = tblocal.getcol("PHASE_DIR")
-149
-+
+        field_frames = _find_field_ref_frames(tblocal)
-151
+        if isinstance(phasecenter, str):
--
+            thenewra_rad, thenewdec_rad = _convert_to_ra_dec_j2000(phasecenter)
-152
+            if field:
-153
+                try:
-154
+                    field_id = int(field)
-155
+                except ValueError as _exc:
-156
+                    fnames = tblocal.getcol("NAME")
-157
+                    field_id = np.where(fnames == field)[0][0]
-+
+                thenewra_rad, thenewdec_rad = _convert_to_ref_frame(phasecenter, field_frames[field_id])
-159
+                pcol[0][0][field_id] = thenewra_rad
-160
+                pcol[1][0][field_id] = thenewdec_rad
-161
+            else:
-162
+                for row in range(0, tblocal.nrows()):
-+
+                    thenewra_rad, thenewdec_rad = _convert_to_ref_frame(phasecenter, field_frames[row])
-164
+                    pcol[0][0][row] = thenewra_rad
-165
+                    pcol[1][0][row] = thenewdec_rad
-166
-167
+        elif isinstance(phasecenter, dict):
-168
+            for field_id, field_center in phasecenter.items():
--
+                thenewra_rad, thenewdec_rad = _convert_to_ra_dec_j2000(field_center)
-169
+                field_iidx = int(field_id)
-+
+                thenewra_rad, thenewdec_rad = _convert_to_ref_frame(field_center, field_frames[field_iidx])
-171
+                pcol[0][0][field_iidx] = thenewra_rad
-172
+                pcol[1][0][field_iidx] = thenewdec_rad
-173
-174
+        tblocal.putcol("PHASE_DIR", pcol)
-175
-176
+    except Exception as instance:
-177
+        casalog.post(f"*** Error '%s' updating FIELD subtable {instance}", "WARN")
-178
+        raise RuntimeError(str(instance))
-179
+    finally:
-180
+        tblocal.done()
-181
-182
--
+def _convert_to_ra_dec_j2000(phasecenter: str) -> Tuple[float, float]:
--
+    """Parse phase center string to obtain ra/dec (in rad)"""
-+
+def _find_field_ref_frames(tblocal: table) -> Dict[int, str]:
-+
+    """
-+
+    Given an open FIELD subtable, returns a dict of {field: reference_frame} for the PHASE_DIR
-+
+    column, where field is of type int, and reference_frame is of type str.
-+
-+
+    This handles:
-+
+    - simple metadata where the reference frame is in the keywords of the PHASE_DIR column (same
-+
+      ref frame for all fields)
-+
+    - variable (per-field) reference frame metadata, usually in a PhaseDir_Ref additional column
-+
+    """
-+
-+
+    dir_col = "PHASE_DIR"
-+
+    metainfo = tblocal.getcolkeyword(dir_col, "MEASINFO")
-+
+    nrows = tblocal.nrows()
-+
+    if "Ref" in metainfo:
-+
+        ref_frame = metainfo["Ref"]
-+
+        field_ref_frames = {field_id: ref_frame for field_id in range(0, nrows)}
-+
+    elif "Ref" not in metainfo and ("VarRefCol" in metainfo and "TabRefTypes" in metainfo and
-+
+                                    "TabRefCodes" in metainfo):
-+
+        col = metainfo["VarRefCol"]   # presumably PhaseDir_Ref
-+
+        ref_frame_codes = tblocal.getcol(col)
-+
+        ref_codes_to_frame_idx = {code: idx for idx, code in enumerate(metainfo["TabRefCodes"])}
-+
+        ref_frame_names = [metainfo["TabRefTypes"][ref_codes_to_frame_idx[code]] for code in
-+
+                           ref_frame_codes]
-+
+        field_ref_frames = {field_id: ref_frame_names[field_id] for field_id in np.arange(0, nrows)}
-+
+    else:
-+
+        raise RuntimeError("Error when retrieving reference frames from the metadata of column "
-+
+                           f"{dir_col}. The field 'Ref' is not present but could not find the "
-+
+                           "fields 'VarRefCol', 'TabRefTypes', and 'TabRefCodes'")
-+
-+
+    return field_ref_frames
-+
-+
-+
+def _convert_to_ref_frame(phasecenter: str, output_ref_frame: str) -> Tuple[float, float]:
-+
+    """
-+
+    Converts one phasecenter (presumably given as input to this task) to another reference
-+
+    frame (presumably the frame used in the (output) MS for the relevant field).
-+
-+
+    When applying phaseshift, the output MS has the same reference frames as the input MS, as
-+
+    propagated by mstransform.
-+
-+
+    Returns the v0,v1 (RA,Dec) values in units of radians and using the requested frame, ready to be
-+
+    written to a FIELD subtable.
-+
+    """
-227
+    def parse_phasecenter(center: str) -> Tuple[str, str, str]:
-228
+        """
--
+        Splits the:
-+
+        Parse phase center string to obtain ra/dec (in rad). Splits the:
-230
+        - (optional) frame,
-231
+        - v0 (typically RA),
-232
+        - v1 (typically Dec)
-233
+        in a phasecenter string.
-234
+        """
-235
+        dirstr = center.split(" ")
-236
+        if 3 == len(dirstr):
-237
+            dir_frame, dir_v0, dir_v1 = dirstr[0], dirstr[1], dirstr[2]
-238
+        elif 2 == len(dirstr):
-239
+            dir_frame = ""  # J2000 will be default in me.direction, etc.
-240
+            dir_v0, dir_v1 = dirstr[0], dirstr[1]
-241
+        else:
-242
+            raise AttributeError(f"Wrong phasecenter string: '{center}'. It must have 3 or 2 "
-243
+                                 "items separated by space(s): 'reference_frame RA Dec' or 'RA Dec'")
-244
-245
+        return dir_frame, dir_v0, dir_v1
-246
-247
+    try:
-248
+        melocal = me()
-249
+        dir_frame, dir_v0, dir_v1 = parse_phasecenter(phasecenter)
-+
-+
+        # Note: even if the input frame is the same as output_ref_frame, we need to ensure units of
-+
+        # radians for the FIELD subtable
-253
+        if dir_frame:
-254
+            thedir = melocal.direction(dir_frame, dir_v0, dir_v1)
-255
+        else:
-256
+            thedir = melocal.direction(v0=dir_v0, v1=dir_v1)
-257
+        if not thedir:
-258
+            raise RuntimeError(
-259
+                f"measures.direction() failed for phasecenter string: {phasecenter}"
-260
+            )
--
+        if dir_frame != "J2000":
--
+            # Convert to J2000
--
+            thedir = melocal.measure(thedir, "J2000")
-+
-+
+        if dir_frame != output_ref_frame:
-+
+            thedir = melocal.measure(thedir, output_ref_frame)
-264
+        thenewra_rad = thedir["m0"]["value"]
-265
+        thenewdec_rad = thedir["m1"]["value"]
-266
+    except Exception as instance:
-267
+        casalog.post(
-268
+            f"*** Error {instance} when interpreting parameter 'phasecenter': ",
-269
+            "SEVERE",
-270
+        )
-271
+        raise RuntimeError(str(instance))
-272
+    finally:
-273
+        melocal.done()

casa6

Commits

Jira Issues

Add shortcut