SPR1-1207: Support wsclean sources

katpoint/catalogue.py

@@ -138,11 +138,12 @@ class Catalogue(object):
     and many *xephem* targets are stored in EDB database files. Editing these
     files to make each line a valid :class:`Target` description string is
     cumbersome, especially in the case of TLE files which are regularly updated.
-    Two special methods simplify the loading of targets from these files::
+    Three special methods simplify the loading of targets from these files::
 
         cat = katpoint.Catalogue()
         cat.add_tle(file('gps-ops.txt'))
         cat.add_edb(file('hipparcos.edb'))
+        cat.add_wsclean(file('<prefix>-sources.txt'))
 
     Whenever targets are added to the catalogue, a tag or list of tags may be
     specified. The tags can also be given as a single string of
@@ -155,6 +156,7 @@ class Catalogue(object):
         cat.add_tle(file('glo-ops.txt'), tags=['glonass', 'satellite'])
         cat.add(file('source_list.csv'), tags='calibrator')
         cat.add_edb(file('hipparcos.edb'), tags='star')
+        cat.add_wsclean(file('<prefix>-sources.txt'), tags='wsc')
 
     Finally, targets may be removed from the catalogue. The most recently added
     target with the specified name is removed from the targets list as well as
@@ -569,6 +571,100 @@ def add_edb(self, lines, tags=None):
             targets.append('xephem,' + line.replace(',', '~'))
         self.add(targets, tags)
 
+    def add_wsclean(self, lines, tags=None):
+        """Add WSClean format (BlackBoard Self-cal DPPP) targets to catalogue.
+
+        Examples of catalogue construction can be found in the :class:`Catalogue`
+        documentation.
+
+        Parameters
+        ----------
+        lines : sequence of strings
+            List of lines containing a target per line (may also be file object)
+        tags : string or sequence of strings, optional
+            Tag or list of tags to add to targets (strings will be split on
+            whitespace)
+
+        Examples
+        --------
+        Here is an example of adding WCS targets to a catalogue:
+
+        >>> from katpoint import Catalogue
+        >>> cat = Catalogue()
+        >>> cat.add_wsclean(file('wsclean-001-sources.txt'), tags='cal')
+        >>> lines = ['Format = Name, Type, Ra, Dec, I, SpectralIndex, LogarithmicSI,
+                    ReferenceFrequency=125584411.621094, MajorAxis, MinorAxis, Orientation',
+                    's0c0,POINT,08:28:05.152,39.35.08.511,0.000748810650400475,\
+                    [-0.00695379313004673,-0.0849693907803257],false,125584411.621094,,,',
+                    's1c1,GAUSSIAN,07:51:09.24,42.32.46.177,0.000660490865128381,\
+                    [0.00404869217508666,-0.011844732049232],false,125584411.621094,\
+                    83.6144111272856,83.6144111272856,0']
+        >>> cat.add_wsclean(lines)
+        """
+
+        targets = []
+        hd = {}
+        for line in lines:
+            if (line[0] == '#') or (len(line.strip()) == 0):
+                continue
+            if line.startswith('Format ='):
+                # if the first line is a format specifier (header), create a dictionary with keys
+                # the field names and default values if they exist
+                line = ''.join(line.split(' '))
+                for item in line[len('Format='):].split(','):
+                    item = item.replace("'", "").strip().split('=')
+                    if len(item) > 1:
+                        hd[item[0]] = item[1]
+                    else:
+                        hd[item[0]] = None
+                # remove SpectralIndex from the dictionary (see below)
+                hd.pop('SpectralIndex')
+                continue
+
+            # check that the format specifier dictionary has been populated
+            if not hd:
+                raise ValueError("WSCLean format not specified in header line")
+
+            # extract (a variable number of) spectral index coefficients from the middle of the
+            # WSClean format string.
+            si = line[line.find('[') + 1:line.find(']')] .replace(',', ' ')
+            line = line[0:line.find('[') - 1] + line[line.find(']') + 1:]
+
+            # create a dictionary from all fields except for si
+            try:
+                wsc_dict = {k: v for k, v in zip(list(hd.keys()), line.strip().split(','))}
+            except KeyError:
+                raise KeyError("malformed or nonexistent header/format specifier in source list")
+
+            # set default values
+            for item in wsc_dict:
+                if not wsc_dict[item]:
+                    if hd[item]:
+                        wsc_dict[item] = hd[item]
+
+            # add back in the SpectralIndex
+            if wsc_dict['LogarithmicSI']:  # spectral index is specified identically to katpoint
+                wsc_dict['SpectralIndex'] = si
+            else:  # TODO: convert between logarithmic and polynomial si (curr. uses just Stokes I)
+                wsc_dict['SpectralIndex'] = si.partition(',')[0]
+
+            # Add wsc source type ('point' | 'gaussian') as an extra tag
+            if wsc_dict['Type'] == 'POINT':
+                tags = 'point'
+            elif wsc_dict['Type'] == 'GAUSSIAN':
+                tags = 'gaussian'
+
+            # convert re-ordered dict to tilde-separated string
+            n_line = ''
+            for key, value in wsc_dict.items():
+                n_line += f'{value}~ '
+            n_line = n_line[:-2]
+
+            targets.append(f'wsclean {tags}, {wsc_dict["Ra"]}, {wsc_dict["Dec"]}, '
+                           f'{wsc_dict["SpectralIndex"]}, {n_line}')
+
+        self.add(targets)
+
     def remove(self, name):
         """Remove target from catalogue.
 

katpoint/target.py

@@ -56,7 +56,7 @@ class Target(object):
     the list may be empty. The <tags> field contains a space-separated list of
     descriptive tags for the target. The first tag is mandatory and indicates
     the body type of the target, which should be one of (*azel*, *radec*, *gal*,
-    *tle*, *special*, *star*, *xephem*).
+    *tle*, *special*, *star*, *xephem*, *wsclean*).
 
     The longitudinal and latitudinal fields are only relevant to *azel*, *radec*
     and *gal* targets, in which case they contain the relevant coordinates. The
@@ -85,11 +85,17 @@ class Target(object):
 
     For *tle* bodies, the final field in the description string should contain
     the three lines of the TLE. If the name list is empty, the target name is
-    taken from the TLE instead. The *xephem* body contains a string in XEphem
+    taken from the TLE instead.
+
+    The *xephem* body contains a string in XEphem
     EDB database format as the final field, with commas replaced by tildes. If
     the name list is empty, the target name is taken from the XEphem string
     instead.
 
+    For *wsclean* bodies, the final field in the description string should
+    contain the wsclean component list format. If the name list is empty, the
+    target name is taken from the WSClean body instead.
+
     When specifying a description string, the rest of the target parameters are
     ignored, except for the default antenna and flux frequency (which do not
     form part of the description string).
@@ -283,6 +289,14 @@ def description(self):
                 fields = [tags]
             fields += [edb_string]
 
+        elif self.body_type == 'wsclean':
+            # TODO: comprehensive wsclean specification
+            wsc_string = '~'.join([wsc_field.strip() for wsc_field in self.body.writedb().split(',')])
+            wsc_name = wsc_string[:wsc_string.index('~')]
+            if wsc_name == names:
+                fields = [tags]
+            fields += [wsc_string]
+
         return ', '.join(fields)
 
     def add_tags(self, tags):
@@ -979,7 +993,7 @@ def construct_target_params(description):
         raise ValueError("Target description '%s' must have at least two fields" % description)
     # Check if first name starts with body type tag, while the next field does not
     # This indicates a missing names field -> add an empty name list in front
-    body_types = ['azel', 'radec', 'gal', 'tle', 'special', 'star', 'xephem']
+    body_types = ['azel', 'radec', 'gal', 'tle', 'special', 'star', 'xephem', 'wsclean']
     if np.any([fields[0].startswith(s) for s in body_types]) and \
        not np.any([fields[1].startswith(s) for s in body_types]):
         fields = [''] + fields
@@ -1001,6 +1015,9 @@ def construct_target_params(description):
     while len(fields[-1]) == 0:
         fields.pop()
 
+    # Extract flux model if it is available
+    flux_model = FluxDensityModel(fields[4]) if (len(fields) > 4) and (len(fields[4].strip(' ()')) > 0) else None
+
     # Create appropriate PyEphem body based on body type
     if body_type == 'azel':
         if len(fields) < 4:
@@ -1106,12 +1123,44 @@ def construct_target_params(description):
         elif edb_type == 'P':
             tags.insert(1, 'special')
 
+    elif body_type == 'wsclean':
+        wsc_string = fields[-1].replace('~', ',')
+        wsc_name_field = wsc_string.partition(',')[0]
+        wsc_names = [name.strip() for name in wsc_name_field.split('|')]
+
+        if preferred_name:
+            wsc_string = wsc_string.replace(wsc_name_field, preferred_name)
+        else:
+            preferred_name = wsc_names[0]
+        if preferred_name != wsc_names[0]:
+            aliases.append(wsc_names[0])
+        for extra_name in wsc_names[1:]:
+            if not (extra_name in aliases) and not (extra_name == preferred_name):
+                aliases.append(extra_name)
+
+        wsc_string_l = wsc_string.split(',')
+        try:
+            wsc_ra_field = wsc_string_l[2]
+            wsc_dec_field = wsc_string_l[3]
+            wsc_flux_field = wsc_string_l[4]
+        except ValueError:
+            raise ValueError("WSClean target description string contains unknown key(s)")
+
+        body = ephem.FixedBody()
+        ra, dec = angle_from_hours(wsc_ra_field), angle_from_degrees(wsc_dec_field)
+        if preferred_name:
+            body.name = preferred_name
+        else:
+            body.name = "Ra: %s Dec: %s" % (ra, dec)
+        body._ra = ra
+        body._dec = dec
+
+        #  Extract wsc flux model
+        flux_model = FluxDensityModel(f'300 4000 {wsc_flux_field.strip()}')
+
     else:
         raise ValueError("Target description '%s' contains unknown body type '%s'" % (description, body_type))
 
-    # Extract flux model if it is available
-    flux_model = FluxDensityModel(fields[4]) if (len(fields) > 4) and (len(fields[4].strip(' ()')) > 0) else None
-
     return body, tags, aliases, flux_model
 
 # --------------------------------------------------------------------------------------------------