Collection of scripts and small programs used by the EVN Support Scientists at JIVE during the regular data processing of EVN observations.
Vous ne pouvez pas sélectionner plus de 25 sujets Les noms de sujets doivent commencer par une lettre ou un nombre, peuvent contenir des tirets ('-') et peuvent comporter jusqu'à 35 caractères.

172 lignes
6.4 KiB

  1. #! /usr/bin/python3
  2. """
  3. Invert the subband for specified antennas.
  4. Usage: test_invert_subband.py msdata antennas [-sb SUBBANDS]
  5. Options:
  6. msdata : str MS data set containing the data to be swapged.
  7. antennas : str Antenna or list of antennas that require to be
  8. swapped. Use the two-letter name (e.g. Ef, Mc,
  9. or Jb2). In case of more than one station, please
  10. use either string 'Ef, Mc' or a non-spaced str:
  11. Ef,Mc,Ys.
  12. Version: 1.0
  13. Date: Nov 2018
  14. Written by Benito Marcote (marcote@jive.eu)
  15. """
  16. import sys
  17. import copy
  18. import time
  19. import argparse
  20. from enum import IntEnum
  21. import datetime as dt
  22. import numpy as np
  23. from pyrap import tables as pt
  24. usage = "%(prog)s [-h] [-v] [-t1 STARTTIME] [-t2 ENDTIME] <measurement set> <antenna>"
  25. description="""Invert the subband for specified antennas.
  26. Fixes the problem when a subband is flipped (increasing frequency instead of decreasing along the
  27. different channels. Observed it in Jb2 during 2018 Session 1 in spectral line observations.
  28. """
  29. help_msdata = 'Measurement Set containing the data to be corrected.'
  30. help_antenna = 'Name of the antenna to be corrected as it appears in the MS (case insensitive).'
  31. help_t1 = 'Start time of the data that need to be corrected. By default the beginning of the observation.'\
  32. +'In Aips format: YYYY/MM/DD/hh:mm:ss or YYYY/DOY/hh:mm:ss'
  33. help_t2 = 'Ending time of the data that need to be corrected. By default the ending of the observation.'\
  34. +'In Aips format: YYYY/MM/DD/hh:mm:ss or YYYY/DOY/hh:mm:ss'
  35. parser = argparse.ArgumentParser(description=description, prog='polswap.py', usage=usage)
  36. parser.add_argument('msdata', type=str, help=help_msdata)
  37. parser.add_argument('antenna', type=str, help=help_antenna)
  38. parser.add_argument('-t1', '--starttime', default=None, type=str, help=help_t1)
  39. parser.add_argument('-t2', '--endtime', default=None, type=str, help=help_t2)
  40. parser.add_argument('-v', '--version', action='version', version='%(prog)s 1.0')
  41. arguments = parser.parse_args()
  42. msdata = arguments.msdata[:-1] if arguments.msdata[-1]=='/' else arguments.msdata
  43. def cli_progress_bar(current_val, end_val, bar_length=40):
  44. percent = current_val/end_val
  45. hashes = '#'*int(round(percent*bar_length))
  46. spaces = ' '*(bar_length-len(hashes))
  47. sys.stdout.write("\rProgress: [{0}] {1}%".format(hashes+spaces, int(round(percent*100))))
  48. sys.stdout.flush()
  49. def chunkert(f, l, cs, verbose=True):
  50. while f<l:
  51. n = min(cs, l-f)
  52. yield (f, n)
  53. f = f + n
  54. class Stokes(IntEnum):
  55. """The Stokes types defined as in the enum class from casacore code.
  56. """
  57. Undefined = 0 # Undefined value
  58. I = 1 # standard stokes parameters
  59. Q = 2
  60. U = 3
  61. V = 4
  62. RR = 5 # circular correlation products
  63. RL = 6
  64. LR = 7
  65. LL = 8
  66. XX = 9 # linear correlation products
  67. XY = 10
  68. YX = 11
  69. YY = 12
  70. RX = 13 # mixed correlation products
  71. RY = 14
  72. LX = 15
  73. LY = 16
  74. XR = 17
  75. XL = 18
  76. YR = 19
  77. YL = 20
  78. PP = 21 # general quasi-orthogonal correlation products
  79. PQ = 22
  80. QP = 23
  81. QQ = 24
  82. RCircular = 25 # single dish polarization types
  83. LCircular = 26
  84. Linear = 27
  85. Ptotal = 28 # Polarized intensity ((Q^2+U^2+V^2)^(1/2))
  86. Plinear = 29 # Linearly Polarized intensity ((Q^2+U^2)^(1/2))
  87. PFtotal = 30 # Polarization Fraction (Ptotal/I)
  88. PFlinear = 31 # linear Polarization Fraction (Plinear/I)
  89. Pangle = 32 # linear polarization angle (0.5 arctan(U/Q)) (in radians)
  90. def atime2datetime(atime):
  91. """Converts a string with the form YYYY/MM/DD/hh:mm:ss or YYYY/DOY/hh:mm:ss to MJD"""
  92. if atime.count('/') == 3:
  93. # Format: YYYY/MM/DD/hh:mm:ss
  94. return dt.datetime.strptime(atime, '%Y/%m/%d/%H:%M:%S')
  95. if atime.count('/') == 2:
  96. # Format: YYYY/DOY/hh:mm:ss
  97. return dt.datetime.strptime(atime, '%Y/%j/%H:%M:%S')
  98. else:
  99. raise ValueError('Date format must be YYYY/MM/DD/hh:mm:ss or YYYY/DOY/hh:mm:ss')
  100. with pt.table(msdata, readonly=False, ack=False) as ms:
  101. # Get CORR_TYPE, the polarization entries in the data
  102. changes = None
  103. with pt.table(ms.getkeyword('ANTENNA'), readonly=True, ack=False) as ms_ant:
  104. antenna_number = [i.upper() for i in ms_ant.getcol('NAME')].index(arguments.antenna.upper())
  105. with pt.table(ms.getkeyword('OBSERVATION'), readonly=True, ack=False) as ms_obs:
  106. time_range = (dt.datetime(1858, 11, 17, 0, 0, 2) + \
  107. ms_obs.getcol('TIME_RANGE')*dt.timedelta(seconds=1))[0]
  108. # Get the timerange to apply to polswap
  109. if arguments.starttime is not None:
  110. datetimes_start = atime2datetime(arguments.starttime)
  111. else:
  112. datetimes_start = time_range[0] - dt.timedelta(seconds=1)
  113. if arguments.endtime is not None:
  114. datetimes_end = atime2datetime(arguments.endtime)
  115. else:
  116. datetimes_end = time_range[1] + dt.timedelta(seconds=1)
  117. columns = ('DATA', 'FLOAT_DATA', 'FLAG', 'SIGMA_SPECTRUM', 'WEIGHT_SPECTRUM',
  118. 'WEIGHT', 'SIGMA')
  119. # shapes of DATA, FLOAT_DATA, FLAG, SIGMA_SPECTRUM, WEIGHT_SPECTRUM: (nrow, npol, nfreq)
  120. # shapes of WEIGHT, SIGMA: (nrow, npol)
  121. columns = [a_col for a_col in columns if a_col in ms.colnames()]
  122. print('\nThe following columns will be modified: {}.\n'.format(', '.join(columns)))
  123. for (start, nrow) in chunkert(0, len(ms), 5000):
  124. cli_progress_bar(start, len(ms), bar_length=40)
  125. for changei, antpos in zip(changes, ('ANTENNA1','ANTENNA2')):
  126. ants = ms.getcol(antpos, startrow=start, nrow=nrow)
  127. datetimes = dt.datetime(1858, 11, 17, 0, 0, 2) + \
  128. ms.getcol('TIME', startrow=start, nrow=nrow)*dt.timedelta(seconds=1)
  129. cond = np.where((ants == antenna_number) & (datetimes > datetimes_start) & (datetimes < datetimes_end))
  130. if len(cond[0]) > 0:
  131. for a_col in columns:
  132. ms_col = ms.getcol(a_col, startrow=start, nrow=nrow)
  133. if len(ms_col.shape) == 3:
  134. ms_col[cond,] = ms_col[cond,][:,::-1,:]
  135. elif len(ms_col.shape) == 2:
  136. ms_col[cond,] = ms_col[cond,][:,::-1]
  137. else:
  138. raise ValueError('Unexpected dimensions for {} column.'.format(a_col))
  139. ms.putcol(a_col, ms_col, startrow=start, nrow=nrow)
  140. print('\n{} modified correctly.'.format(msdata))