Collection of scripts and small programs used by the EVN Support Scientists at JIVE during the regular data processing of EVN observations.
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

133 lines
5.5 KiB

#!/usr/bin/env python3
4 years ago
Flag visibilities with weights below the provided threshold.
Usage: msdata threshold
msdata : str MS data set containing the data to be flagged.
threshold : float Visibilities with a weight below the specified
value will be flagged. Must be positive.
Version: 3.0
Date: Apr 2019
4 years ago
Written by Benito Marcote (
version 3.0 changes
- Refactoring code (thanks to Harro).
version 2.0 changes
- Major revision. Now it does not modify the weights anymore. Instead, it
flags those data with weights below the given threshold by modifying the
FLAG table.
- Small change in print messages to show '100%' instead of '1e+02%' in certain
version 1.4 changes
- Now it also reports the percentage or data that were different from
zero and will be flagged (not only the total data as before).
version 1.3 changes
- Minor fixes (prog name in optparse info).
4 years ago
version 1.2 changes
- Minor fixes.
version 1.1 changes
- Added option -v that allows you to just get how many visibilities will
be flagged (but without actually flagging the data).
4 years ago
from pyrap import tables as pt
import numpy as np
import sys
__version__ = 3.0
4 years ago
help_msdata = 'Measurement set containing the data to be corrected.'
help_threshold = 'Visibilities with a weight below this value will be flagged. Must be positive.'
help_v = 'Only checks the visibilities to flag (do not flag the data).'
4 years ago
usage = "%(prog)s [-h] <measurement set> <weight threshold>"
description="""Flag visibilities with weights below the provided threshold.
import argparse
parser = argparse.ArgumentParser(description=description, prog='', usage=usage)
parser.add_argument('msdata', type=str, help=help_msdata)
parser.add_argument('threshold', type=float, help=help_threshold)
parser.add_argument('--version', action='version', version='%(prog)s {}'.format(__version__))
parser.add_argument("-v", "--verbose", default=True, action="store_false" , help=help_v)
4 years ago
arguments = parser.parse_args()
4 years ago
#print('The arguments ', arguments)
verbose = arguments.verbose
4 years ago
msdata = arguments.msdata[:-1] if arguments.msdata[-1]=='/' else arguments.msdata
threshold = arguments.threshold
except ImportError:
4 years ago
usage = "%prog [-h] [-v] <measurement set> <weight threshold>"
4 years ago
description="""Flag visibilities with weights below the provided threshold.
# Compatibility with Python 2.7 in eee
import optparse
parser = optparse.OptionParser(usage=usage, description=description, prog='', version='%prog 1.3')
parser.add_option("-v", action="store_false", dest="verbose", default=True, help=help_v)
theparser = parser.parse_args()
verbose = theparser[0].verbose
arguments = theparser[1]
#arguments = parser.parse_args()[1]
4 years ago
if len(arguments) != 2:
4 years ago
print('Two arguments must be provided: [-h] [-v] <measurement set> <weight threshold>')
print('Use -h to get help.')
4 years ago
msdata = arguments[0][:-1] if arguments[0][-1]=='/' else arguments[0]
threshold = float(arguments[1])
assert threshold > 0.0
def chunkert(f, l, cs, verbose=True):
while f<l:
n = min(cs, l-f)
yield (f, n)
f = f + n
percent = lambda x, y: (float(x)/float(y)) * 100.0
with pt.table(msdata, readonly=False, ack=False) as ms:
total_number = 0
flagged_before, flagged_after = (0, 0)
flagged_nonzero, flagged_nonzero_before, flagged_nonzero_after = (0, 0, 0)
# WEIGHT: (nrow, npol)
# WEIGHT_SPECTRUM: (nrow, npol, nfreq)
# flags[weight < threshold] = True
weightcol = 'WEIGHT_SPECTRUM' if 'WEIGHT_SPECTRUM' in ms.colnames() else 'WEIGHT'
transpose = (lambda x:x) if weightcol == 'WEIGHT_SPECTRUM' else (lambda x: x.transpose((1, 0, 2)))
for (start, nrow) in chunkert(0, len(ms), 5000):
# shape: (nrow, npol, nfreq)
flags = transpose(ms.getcol("FLAG", startrow=start, nrow=nrow))
total_number += np.product(flags.shape)
# count how much data is already flagged
flagged_before += np.sum(flags)
# extract weights and compute new flags based on threshold
weights = ms.getcol(weightcol, startrow=start, nrow=nrow)
# how many non-zero did we flag
flagged_nonzero_before = np.logical_and(flags, weights > 0)
# join with existing flags and count again
flags = np.logical_or(flags, weights < threshold)
flagged_after += np.sum(flags)
flagged_nonzero_after = np.logical_and(flags, weights > 0)
flagged_nonzero += np.sum(np.logical_xor(flagged_nonzero_before, flagged_nonzero_after))
# one thing left to do: write the updated flags to disk
#flags = ms.putcol("FLAG", flags.transpose((1, 0 , 2)), startrow=start, nrow=nrow)
if verbose:
flags = ms.putcol("FLAG", transpose(flags), startrow=start, nrow=nrow)
print("Got {0:11} visibilities".format(total_number))
print("Got {0:11} visibilities to flag using threshold {1}\n".format(flagged_after-flagged_before,
print("{0:.2f}% total vis. flagged ({2:.2f}% to flag in this execution).\n{1:.2f}% data with non-zero weights flagged.\n".format(percent(flagged_after, total_number), percent(np.sum(flagged_nonzero_after), total_number), percent(flagged_after-flagged_before, total_number)))
4 years ago
if verbose:
print('Flags have not been applied.')
4 years ago