inspect_flagging_dataloss.cwl

class: CommandLineTool
cwlVersion: v1.1
id: inspect_flagging_dataloss
baseCommand:
  - python3
  - script.py
inputs:
  - id: input
    type: Directory
    inputBinding:
      position: 0
      shellQuote: false
  - id: uvw_sampling
    type: int?
    default: 100
    inputBinding:
      position: 1
outputs:
  - id: output
    type: 'File'
    outputBinding:
      glob: metrics.json
  - id: flags_output
    type: 'File'
    outputBinding:
      glob: flags.h5
label: inspect_flagging_dataloss
requirements:
  - class: ShellCommandRequirement
  - class: InitialWorkDirRequirement
    listing:
     - entryname: script.py
       entry: |
         import sys
         from casacore.tables import table
         import itertools
         import json
         import h5py
         import numpy as np
         import re
         from scipy.constants import c as LIGHT_SPEED
         input_ms_path = sys.argv[1]
         sampling = int(sys.argv[2])
         main_table = table(input_ms_path, 'r')

         antennas = table(input_ms_path + '/ANTENNA', 'r')
         spec_window = table(input_ms_path + '/SPECTRAL_WINDOW', 'r')
         observation = table(input_ms_path + '/OBSERVATION', 'r')

         sas_id, sap, sb = re.search('L(\d*).*SAP(\d{3}).*SB(\d{3}).*',
                   observation.getcell('LOFAR_FILENAME', 0)).groups()
         antenna_names = antennas.getcol('NAME')
         antenna_positions = antennas.getcol('POSITION')

         ref_frequency = spec_window.getcell('REF_FREQUENCY', 0)

         TO_UVW = ref_frequency / LIGHT_SPEED  # 1 / (c/nu)

         flags = np.array(main_table.getcol('FLAG'))
         vis = np.array(main_table.getcol('DATA'))

         uvw = np.array(main_table.getcol('UVW') * TO_UVW)

         antenna1 = np.array(main_table.getcol('ANTENNA1'))
         antenna2 = np.array(main_table.getcol('ANTENNA2'))
         time = np.array(main_table.getcol('TIME'))
         n_antennas = antennas.nrows()
         n_baselines = int((n_antennas + 1) * n_antennas * .5)
         n_times = flags.shape[0] // n_baselines
         dataloss = np.array(vis == 0, dtype=np.int)
         urange = np.linspace(-np.max(uvw[:, 0]), np.max(uvw[:, 0]), sampling, endpoint=True)
         vrange = np.linspace(-np.max(uvw[:, 1]), np.max(uvw[:, 1]), sampling, endpoint=True)
         
         dataloss_matrix, (urange, vrange, wrange) = np.histogramdd(uvw, bins=[urange, vrange, 1], weights=np.nanmean(dataloss, axis=(1,2)))
         dataloss_matrix += np.histogramdd(-uvw, bins=[urange, vrange, 1], weights=np.nanmean(dataloss, axis=(1,2)))[0]
         
         flag_matrix, (urange, vrange, wrange) = np.histogramdd(uvw, bins=[urange, vrange, 1], weights=np.nanmean(flags, axis=(1,2)))
         flag_matrix += np.histogramdd(-uvw, bins=[urange, vrange, 1], weights=np.nanmean(flags, axis=(1,2)))[0]
         
         coverage, (urange, vrange, wrange) = np.histogramdd(uvw, bins=[urange, vrange, 1])
         coverage += np.histogramdd(-uvw, bins=[urange, vrange, 1])[0]

         flags = flags.reshape(n_times, n_baselines, *flags.shape[1:], order='C')
         dataloss = dataloss.reshape(n_times, n_baselines, *flags.shape[2:], order='C')
         uvw = uvw.reshape(n_times, n_baselines, *uvw.shape[1:], order='C')
         timestamp = time.reshape(n_times, n_baselines, order='C')[:, 0]

         baselines = [
          "%s,%s" % (antenna_names[i], antenna_names[j])
          for i, j in itertools.product(range(n_antennas), range(n_antennas)) if i >= j]

         dataloss_as_function_of_time = np.nanmean(dataloss, axis=(1, 2, 3))
         dataloss_as_function_of_baseline = np.nanmean(dataloss, axis=(0, 2, 3))


         flags_as_function_of_time = np.nanmean(flags, axis=(1,2,3))
         flags_as_function_of_baseline = np.nanmean(flags, axis=(0,2,3))

         metrics = dict(baselines=baselines,
                        dataloss_vs_time=dataloss_as_function_of_time.tolist(),
                        dataloss_vs_baseline=dataloss_as_function_of_baseline.tolist(),
                        flags_vs_time=flags_as_function_of_time.tolist(),
                        flags_vs_baseline=flags_as_function_of_baseline.tolist(),
                        timestamp=timestamp.tolist(),
                        sap=sap,
                        subband=sb
                        )

         with open('metrics.json', 'w') as f_stream:
           json.dump(metrics, f_stream)

         with h5py.File('flags.h5', 'w') as flags_file:
           flags_file['/FLAGS'] = flag_matrix
           flags_file['/COVERAGE'] = coverage
           flags_file['/DATALOSS'] = dataloss_matrix

           flags_file['/'].attrs['REF_FREQUENCY'] = ref_frequency
           flags_file['/'].attrs['BASELINES'] = baselines
           flags_file['/'].attrs['PROJECT'] = observation.getcell('PROJECT', 0)
           flags_file['/'].attrs['TELESCOPE'] = observation.getcell('TELESCOPE_NAME', 0)
           flags_file['/'].attrs['SAS_ID'] = sas_id
           flags_file['/'].attrs['SAP'] = sap
           flags_file['/'].attrs['SUBBAND'] = sb

           flags_file['/URANGE'] = urange
           flags_file['/VRANGE'] = vrange
           flags_file['/WRANGE'] = wrange


  - class: DockerRequirement
    dockerPull: git.astron.nl:5000/ldv/imaging_compress_pipeline:v0.5.10
  - class: NetworkAccess
    networkAccess: true