casa / casa6 / 601dacfebee

Commits

Federico Montesino Pouzols authored and Ville Suoranta committed 601dacfebee Merge 24 Jan 2025

Pull request #907: CAS-13616: Add multi-field specification of phasecenter for task_phaseshift

Merge in CASA/casa6 from CAS-13616 to master

* commit '38e1a8f2031e0cc25b5d6ce1ae8ceb2f6230bac9': (32 commits)
  update and extend tests for output center values in same frame as metadata, CAS-13616
  don't convert to J2000, write phasecenter vals using same frame as in the PHASE_DIR metadata, CAS-13616
  fix bug re empty phasecenter frame, add tests, CAS-13616
  clarify default phasecenter frame J2000, CAS-13616
  use typing.Tuple for Python 3.8, CAS-13616
  prefer f-strings to concats, CAS-13616
  more formatting and indentation fixes from code review, try to get rid of tabs, CAS-13616
  re-format with black, CAS-13616
  suggestions from code review, more type hints prevent small leak, CAS-13616
  fixes from code review, prevent small leak, CAS-13616
  improve wording, from code review, CAS-13616
  parse dicts, reorganize tvi config, shift depending on field, CAS-13616
  test cases with dict and field selection, share common output checks with other classes, CAS-13616
  check field subtable in all multified tests, check vis in unchanged fields, CAS-13616
  extend and reorganize output checks in phaseshift_test cases, CAS-13616
  check in/out FIELD subt, extend multified tests with better MS, CAS-13616
  fix _update_field_subt when single phasecenter and field selection used together, CAS-13616
  reduce tool try contexts and indentation levels, CAS-13616
  remove assumptions that break with reindexing off, CAS-13616
  remove duplicate phasecenter check, defer to TVI, CAS-13616
  ...

casatasks/src/private/task_phaseshift.py

Modified

-+
+from typing import Dict, Optional, Tuple, Union
-2
--
-+
+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(
--
+    vis=None, outputvis=None, keepmms=None, field=None,
--
+    spw=None, scan=None, intent=None, array=None,
--
+    observation=None, datacolumn=None, phasecenter=None
-+
+    vis: str,
-+
+    outputvis: str,
-+
+    keepmms: bool,
-+
+    field: Optional[str],
-+
+    spw: Optional[str],
-+
+    scan: Optional[str],
-+
+    intent: Optional[str],
-+
+    array: Optional[str],
-+
+    observation: Optional[str],
-+
+    datacolumn: Optional[str],
-+
+    phasecenter: Union[str, dict],
-23
+):
-24
+    """
-25
+    Changes the phase center for either short or large
-26
+    offsets/angles w.r.t. the original
-27
+    """
--
+    casalog.origin('phaseshift')
-+
+    casalog.origin("phaseshift")
-29
-30
+    if len(phasecenter) == 0:
--
+        raise ValueError('phasecenter parameter must be specified')
-+
+        raise ValueError("phasecenter parameter must be specified")
-32
+    # Initiate the helper class
-33
+    pdh = ParallelDataHelper("phaseshift", locals())
-34
-35
+    # Validate input and output parameters
-36
+    try:
-37
+        pdh.setupIO()
-38
+    except Exception as instance:
--
+        casalog.post(str(instance), 'ERROR')
-+
+        casalog.post(str(instance), "ERROR")
-40
+        raise RuntimeError(str(instance))
-41
-42
+    # Input vis is an MMS
-43
+    if pdh.isMMSAndNotServer(vis) and keepmms:
-44
+        if not pdh.validateInputParams():
--
+            raise RuntimeError('Unable to continue with MMS processing')
-+
+            raise RuntimeError("Unable to continue with MMS processing")
-46
--
+        pdh.setupCluster('phaseshift')
-+
+        pdh.setupCluster("phaseshift")
-48
-49
+        # Execute the jobs
-50
+        try:
-51
+            pdh.go()
-52
+        except Exception as instance:
--
+            casalog.post(str(instance), 'ERROR')
-+
+            casalog.post(str(instance), "ERROR")
-54
+            raise RuntimeError(str(instance))
-55
+        return
-56
--
+    # Create local copies of tools (has to be here, otherwise
--
+    # ParallelDataHelper has a porblem digest the locals
-+
+    # Actual task code starts here
-+
+    # Gather all the parameters in a dictionary.
-+
+    config = {}
-+
-+
+    config = pdh.setupParameters(
-+
+        inputms=vis,
-+
+        outputms=outputvis,
-+
+        field=field,
-+
+        spw=spw,
-+
+        array=array,
-+
+        scan=scan,
-+
+        intent=intent,
-+
+        observation=observation,
-+
+    )
-+
-+
+    colnames = _get_col_names(vis)
-+
+    # Check if CORRECTED column exists, when requested
-+
+    datacolumn = datacolumn.upper()
-+
+    if datacolumn == "CORRECTED":
-+
+        if "CORRECTED_DATA" not in colnames:
-+
+            casalog.post(
-+
+                "Input data column CORRECTED_DATA does not exist. Will use DATA", "WARN"
-+
+            )
-+
+            datacolumn = "DATA"
-+
-+
+    casalog.post(f"Will use datacolumn = {datacolumn}", "DEBUG")
-+
+    config["datacolumn"] = datacolumn
-+
-+
+    # Call MSTransform framework with tviphaseshift=True
-+
+    config["tviphaseshift"] = True
-+
+    config["reindex"] = False
-+
+    tviphaseshift_config = {"phasecenter": phasecenter}
-+
+    config["tviphaseshiftlib"] = tviphaseshift_config
-+
-+
+    # Configure the tool
-+
+    casalog.post(str(config), "DEBUG1")
-93
--
+    tblocal = table()
--
+    mslocal = ms()
-94
+    mtlocal = mstransformer()
--
+    melocal = me()
-+
+    try:
-+
+        mtlocal.config(config)
-97
--
+    # Actual task code starts here
-+
+        # Open the MS, select the data and configure the output
-+
+        mtlocal.open()
-+
-+
+        # Run the tool
-+
+        casalog.post("Shift phase center")
-+
+        mtlocal.run()
-+
+    finally:
-+
+        mtlocal.done()
-+
-+
+    # Write history to output MS, not the input ms.
-108
+    try:
--
+        dirstr = phasecenter.split(' ')
--
+        if not melocal.direction(dirstr[0], dirstr[1], dirstr[2]):
--
+            raise ValueError("Illegal phacecenter specification " + phasecenter)
--
+        try:
--
+            # Gather all the parameters in a dictionary.
--
+            config = {}
-+
+        mslocal = ms()
-+
+        param_names = phaseshift.__code__.co_varnames[: phaseshift.__code__.co_argcount]
-+
+        vars = locals()
-+
+        param_vals = [vars[p] for p in param_names]
-+
+        casalog.post("Updating the history in the output", "DEBUG1")
-+
+        write_history(
-+
+            mslocal, outputvis, "phaseshift", param_names, param_vals, casalog
-+
+        )
-+
+    except Exception as instance:
-+
+        casalog.post(f"*** Error {instance} updating HISTORY", "WARN")
-+
+        raise RuntimeError(str(instance))
-+
+    finally:
-+
+        mslocal.done()
-122
--
+            config = pdh.setupParameters(
--
+                inputms=vis, outputms=outputvis, field=field,
--
+                spw=spw, array=array, scan=scan, intent=intent,
--
+                observation=observation
--
+            )
-+
+    casalog.post(
-+
+        "Updating the FIELD subtable of the output MeasurementSet with shifted"
-+
+        " phase centers",
-+
+        "INFO",
-+
+    )
-+
+    _update_field_subtable(outputvis, field, phasecenter)
-129
--
+            # Check if CORRECTED column exists, when requested
--
+            datacolumn = datacolumn.upper()
--
+            if datacolumn == 'CORRECTED':
--
+                tblocal.open(vis)
--
+                if 'CORRECTED_DATA' not in tblocal.colnames():
--
+                    casalog.post(
--
+                        'Input CORRECTED_DATA does not exist. Will use DATA',
--
+                        'WARN'
--
+                    )
--
+                    datacolumn = 'DATA'
--
+                tblocal.close()
-130
--
+            casalog.post('Will use datacolumn = ' + datacolumn, 'DEBUG')
--
+            config['datacolumn'] = datacolumn
-+
+def _get_col_names(vis: str) -> np.ndarray:
-+
+    tblocal = table()
-+
+    try:
-+
+        tblocal.open(vis)
-+
+        colnames = tblocal.colnames()
-+
+    finally:
-+
+        tblocal.done()
-+
+    return colnames
-139
--
+            # Call MSTransform framework with tviphaseshift=True
--
+            config['tviphaseshift'] = True
--
+            tviphaseshift_config = {}
--
+            tviphaseshift_config['phasecenter'] = phasecenter
--
+            config['tviphaseshiftlib'] = dict(tviphaseshift_config)
-140
--
+            # Configure the tool
--
+            casalog.post(str(config), 'DEBUG1')
--
+            mtlocal.config(config)
-+
+def _update_field_subtable(outputvis: str, field: str, phasecenter: Union[str, dict]):
-+
+    """Update MS/FIELD subtable with shifted center(s)."""
-143
--
+            # Open the MS, select the data and configure the output
--
+            mtlocal.open()
-+
+    try:
-+
+        tblocal = table()
-+
+        # modify FIELD table
-+
+        tblocal.open(outputvis + "/FIELD", nomodify=False)
-+
+        pcol = tblocal.getcol("PHASE_DIR")
-149
--
+            # Run the tool
--
+            casalog.post('Shift phase center')
--
+            mtlocal.run()
-+
+        field_frames = _find_field_ref_frames(tblocal)
-+
+        if isinstance(phasecenter, str):
-+
+            if field:
-+
+                try:
-+
+                    field_id = int(field)
-+
+                except ValueError as _exc:
-+
+                    fnames = tblocal.getcol("NAME")
-+
+                    field_id = np.where(fnames == field)[0][0]
-+
+                thenewra_rad, thenewdec_rad = _convert_to_ref_frame(phasecenter, field_frames[field_id])
-+
+                pcol[0][0][field_id] = thenewra_rad
-+
+                pcol[1][0][field_id] = thenewdec_rad
-+
+            else:
-+
+                for row in range(0, tblocal.nrows()):
-+
+                    thenewra_rad, thenewdec_rad = _convert_to_ref_frame(phasecenter, field_frames[row])
-+
+                    pcol[0][0][row] = thenewra_rad
-+
+                    pcol[1][0][row] = thenewdec_rad
-166
--
+            mtlocal.done()
-+
+        elif isinstance(phasecenter, dict):
-+
+            for field_id, field_center in phasecenter.items():
-+
+                field_iidx = int(field_id)
-+
+                thenewra_rad, thenewdec_rad = _convert_to_ref_frame(field_center, field_frames[field_iidx])
-+
+                pcol[0][0][field_iidx] = thenewra_rad
-+
+                pcol[1][0][field_iidx] = thenewdec_rad
-173
--
+        except Exception as instance:
--
+            mtlocal.done()
--
+            casalog.post(str(instance), 'ERROR')
--
+            raise RuntimeError(str(instance))
-+
+        tblocal.putcol("PHASE_DIR", pcol)
-175
--
+        # Write history to output MS, not the input ms.
--
+        try:
--
+            param_names = phaseshift.__code__.co_varnames[
--
+                :phaseshift.__code__.co_argcount
--
+            ]
--
+            vars = locals()
--
+            param_vals = [vars[p] for p in param_names]
--
--
+            casalog.post('Updating the history in the output', 'DEBUG1')
--
+            write_history(
--
+                mslocal, outputvis, 'phaseshift', param_names,
--
+                param_vals, casalog
--
+            )
--
+        except Exception as instance:
--
+            casalog.post(
--
+                "*** Error " + str(instance)
--
+                + " updating HISTORY", 'WARN'
--
+            )
--
+            raise RuntimeError(str(instance))
-+
+    except Exception as instance:
-+
+        casalog.post(f"*** Error '%s' updating FIELD subtable {instance}", "WARN")
-+
+        raise RuntimeError(str(instance))
-+
+    finally:
-+
+        tblocal.done()
-181
--
+        # Update field table
--
+        try:
-182
--
+            # Parse phase center string to obtain ra/dec
--
+            dirstr = phasecenter.split(' ')
--
+            try:
--
+                thedir = melocal.direction(dirstr[0], dirstr[1], dirstr[2])
--
+                if (dirstr[0] != 'J2000'):
--
+                    # Convert to J2000
--
+                    thedir = melocal.measure(thedir, 'J2000')
--
+                thenewra_rad = thedir['m0']['value']
--
+                thenewdec_rad = thedir['m1']['value']
--
+            except Exception as instance:
--
+                casalog.post(
--
+                    "*** Error " + str(instance)
--
+                    + " when interpreting parameter \'phasecenter\': ",
--
+                    'SEVERE'
--
+                )
--
+                raise RuntimeError(str(instance))
--
--
+            # modify FIELD table
--
+            tblocal.open(outputvis + '/FIELD', nomodify=False)
--
+            pcol = tblocal.getcol('PHASE_DIR')
--
+            for row in range(0, tblocal.nrows()):
--
+                pcol[0][0][row] = thenewra_rad
--
+                pcol[1][0][row] = thenewdec_rad
--
+            tblocal.putcol('PHASE_DIR', pcol)
--
+        except Exception as instance:
--
+            casalog.post(
--
+                "*** Error \'%s\' updating FIELD" + str(instance), 'WARN'
-+
+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.
-+
+    """
-+
+    def parse_phasecenter(center: str) -> Tuple[str, str, str]:
-+
+        """
-+
+        Parse phase center string to obtain ra/dec (in rad). Splits the:
-+
+        - (optional) frame,
-+
+        - v0 (typically RA),
-+
+        - v1 (typically Dec)
-+
+        in a phasecenter string.
-+
+        """
-+
+        dirstr = center.split(" ")
-+
+        if 3 == len(dirstr):
-+
+            dir_frame, dir_v0, dir_v1 = dirstr[0], dirstr[1], dirstr[2]
-+
+        elif 2 == len(dirstr):
-+
+            dir_frame = ""  # J2000 will be default in me.direction, etc.
-+
+            dir_v0, dir_v1 = dirstr[0], dirstr[1]
-+
+        else:
-+
+            raise AttributeError(f"Wrong phasecenter string: '{center}'. It must have 3 or 2 "
-+
+                                 "items separated by space(s): 'reference_frame RA Dec' or 'RA Dec'")
-+
-+
+        return dir_frame, dir_v0, dir_v1
-+
-+
+    try:
-+
+        melocal = me()
-+
+        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
-+
+        if dir_frame:
-+
+            thedir = melocal.direction(dir_frame, dir_v0, dir_v1)
-+
+        else:
-+
+            thedir = melocal.direction(v0=dir_v0, v1=dir_v1)
-+
+        if not thedir:
-+
+            raise RuntimeError(
-+
+                f"measures.direction() failed for phasecenter string: {phasecenter}"
-260
+            )
--
+            raise RuntimeError(str(instance))
-261
-+
+        if dir_frame != output_ref_frame:
-+
+            thedir = melocal.measure(thedir, output_ref_frame)
-+
+        thenewra_rad = thedir["m0"]["value"]
-+
+        thenewdec_rad = thedir["m1"]["value"]
-+
+    except Exception as instance:
-+
+        casalog.post(
-+
+            f"*** Error {instance} when interpreting parameter 'phasecenter': ",
-+
+            "SEVERE",
-+
+        )
-+
+        raise RuntimeError(str(instance))
-272
+    finally:
--
+        if tblocal:
--
+            tblocal.done()
--
+            tblocal = None
--
+        if mslocal:
--
+            mslocal.done()
--
+            mslocal = None
--
+        if mtlocal:
--
+            mtlocal.done()
--
+            mtlocal = None
--
+        if melocal:
--
+            melocal.done()
--
+            melocal = None
-+
+        melocal.done()
-274
-+
+    return thenewra_rad, thenewdec_rad

casa6

Commits

Jira Issues

Add shortcut