SSB-42: reading function for holography dataset

.gitattributes

@@ -26,6 +26,9 @@ CAL/CalibrationCommon/test/t_holography_dataset_class.sh -text
 CAL/CalibrationCommon/test/t_holography_ms_class.py -text
 CAL/CalibrationCommon/test/t_holography_ms_class.run -text
 CAL/CalibrationCommon/test/t_holography_ms_class.sh -text
+CAL/CalibrationCommon/test/t_holography_observation.py -text
+CAL/CalibrationCommon/test/t_holography_observation.run -text
+CAL/CalibrationCommon/test/t_holography_observation.sh -text
 CEP/Calibration/BBSControl/include/BBSControl/CommandHandlerEstimator.h -text
 CEP/Calibration/BBSControl/include/BBSControl/CommandHandlerReducer.h -text
 CEP/Calibration/BBSControl/include/BBSControl/OptionParser.h -text

CAL/CalibrationCommon/lib/CMakeLists.txt

@@ -21,7 +21,8 @@ python_install(
     __init__.py
     datacontainers/holography_dataset.py
     datacontainers/holography_specification.py
-    datacontainers/measurementset.py
+    datacontainers/holography_observation.py
+    datacontainers/holography_measurementset.py
     datacontainers/__init__.py
     coordinates.py
     utils.py

CAL/CalibrationCommon/lib/datacontainers/holography_dataset.py

 from .holography_specification import HolographySpecification
 from lofar.calibration.common.datacontainers.holography_observation import HolographyObservation
 import logging
+import numpy
 
 logger = logging.getLogger(__name__)
 
@@ -17,10 +18,12 @@ class HolographyDataset():
         self.end_time = None  # date time when the observation started in MJD in seconds (float)
         self.rotation_matrix = None  # array(3,3), translation matrix for
         # (RA, DEC) <-> (l, m) conversion
+
+        self.beamlets = None # list of beamlet numbers
         self.antenna_field_position = []  # coordinates of the antenna position in the target
         # station
         self.reference_stations = set()  # set of reference station names
-        self.frequencies = set()  # set of frequencies
+        self.frequencies = list()  # list of frequencies
         self.ra_dec = None # array(Nfrequency, Nbeamlets, 2) contains the ra_dec of which a beam
         # points at given a frequency and a beamlet number
         self.data = None # array(NreferenceStations, Nfrequencies, Nbeamlets) that contains the
@@ -36,6 +39,36 @@ class HolographyDataset():
 
         """
         self.__collect_preliminary_information(station_name, list_of_hbs_ms_tuples)
+        self.__read_data(station_name, list_of_hbs_ms_tuples)
+
+
+    def __read_data(self, station_name, list_of_hbs_ms_tuples):
+        """
+
+        :param station_name:
+        :param list_of_hbs_ms_tuples:
+        :type list_of_hbs_ms_tuples: list[(HolographySpecification, HolographyObservation)]
+        :return:
+        """
+
+        n_frequencies = len(self.frequencies)
+        n_beamlets = len(self.beamlets)
+        n_reference_stations = len(self.reference_stations)
+
+        data = dict()
+        for hbs, ho in list_of_hbs_ms_tuples:
+            if station_name in hbs.target_station_names:
+                frequency = ho.frequency
+
+
+
+                for beamlet in self.beamlets:
+                    cross_correlation, timestamp, flags = ho.ms_for_a_given_beamlet_number[beamlet].\
+                        read_cross_correlation_time_flags_per_station_names(station_name,
+                                                                 self.reference_stations)
+
+                    #for reference_station in enumerate(self.reference_stations):
+                    #data[(reference_station, )]
 
     def __collect_preliminary_information(self, station_name, list_of_hbs_ho_tuples):
         """
@@ -53,6 +86,12 @@ class HolographyDataset():
         source_name = set()
         source_position = set()
         target_stations = set()
+        reference_stations = set()
+        beamlets = set()
+        virtual_pointing = dict()
+        frequencies = set()
+        start_mjd = None
+        end_mjd = None
         for hbs, ho in list_of_hbs_ho_tuples:
 
             target_stations.update(hbs.target_station_names)
@@ -70,20 +109,46 @@ class HolographyDataset():
                          beam_specification.station_pointing['dec'],
                          beam_specification.station_pointing['coordinate_system']
                         ))
+                    if start_mjd is None or start_mjd > ho.start_mjd:
+                        start_mjd = ho.start_mjd
+
+                    if end_mjd is None or end_mjd < ho.end_mjd:
+                        end_mjd = ho.end_mjd
+
+                    frequencies.add(ho.frequency)
 
-                    self.start_time = ho.start_mjd
-                    self.end_time = ho.end_mjd
-                    self.frequencies.add(ho.frequency)
                     self.sas_ids.add(ho.sas_id)
 
                     self.target_station_name = station_name
-                    self.reference_stations.update(hbs.reference_station_names)
-
+                    reference_stations.update(hbs.reference_station_names)
+                    try:
+                        single_beamlet = int(beam_specification.beamlets)
+                    except ValueError as e:
+                        logger.exception('Target station specification incorrect')
+                        raise e
+
+                    beamlets.add(single_beamlet)
+
+                    virtual_pointing [(ho.frequency, single_beamlet)] =\
+                        (beam_specification.virtual_pointing['ra'],
+                         beam_specification.virtual_pointing['dec'])
+            else:
+                continue
+        self.frequencies = sorted(frequencies)
+        self.beamlets = sorted(beamlets)
+        self.start_time = start_mjd
+        self.end_time = end_mjd
 
+        self.reference_stations = list(reference_stations)
 
+        n_frequencies = len(self.frequencies)
+        n_beamlets = len(beamlets)
+        ra_dec = numpy.zeros((n_frequencies, n_beamlets, 2))
+        for frequency_index, frequency in enumerate(self.frequencies):
+            for beamlet_index, beamlet in enumerate(self.beamlets):
+                ra_dec[frequency_index, beamlet_index, :] = \
+                    virtual_pointing[(frequency, beamlet)]
 
-            else:
-                continue
         # reads the target station position and the coordinate of its axes
         # and does this only once since the coordinate will not change
         first_holography_observation = list_of_hbs_ho_tuples[0][1]

CAL/CalibrationCommon/lib/datacontainers/holography_measurementset.py

@@ -2,6 +2,7 @@ from casacore.tables import table as MS_Table
 
 import os
 import re
+import numpy
 
 
 class HolographyMeasurementSet(object):
@@ -118,27 +119,60 @@ class HolographyMeasurementSet(object):
         return self.__extract_source_name_from_pointing()
 
 
-    def read_cross_correlation_per_station_names(self, reference, target):
+    def read_cross_correlation_time_flags_per_station_names(self, target_station, reference_stations):
+        """
+        Read the crosscorrelation for a given station name and a list of reference stations
+
+        :param reference_stations:
+        :type reference_stations: list[str]
+        :param target_station:
+        :type target_station: str
+        :return:
+        """
         data_table = self.get_data_table()
         antennas_table = self.get_antenna_table()
 
         try:
-            antenna_name_id_map = {name:i for i, name in enumerate(antennas_table.getcol('NAME'))}
-            antenna1_list = data_table.getcol('ANTENNA1')
-            antenna2_list = data_table.getcol('ANTENNA2')
-            timestamp = data_table.getcol('TIME')
-            cross_correlation = data_table.getcol('DATA')
-            reference_antenna_id = antenna_name_id_map[reference]
-            target_antenna_id = antenna_name_id_map[target]
-
-            selected_data = [index for index, (a_i, a_j) in enumerate(zip(antenna1_list, antenna2_list))
-                             if a_i == reference_antenna_id and a_j == target_antenna_id]
+
+            antenna_name_index = antennas_table.index('NAME')
+            reference_antennas_indexes = [antenna_name_index[reference_station]
+                                          for reference_station in reference_stations]
+            target_antenna_index = antenna_name_index[target_station]
+            TAQL_query_syntax = '(ANTENNA1 = $target_antenna_index ' \
+                                'OR' \
+                                ' ANTENNA2 = $target_antenna_index)' \
+                                ' AND ' \
+                                '(ANTENNA1 IN $reference_antennas_indexes ' \
+                                'OR' \
+                                ' ANTENNA2 IN $reference_antennas_indexes)'
+
+            table = data_table.query(TAQL_query_syntax)
+
+            timestamps = list(table.getcol('TIME'))
+            flags = table.getcol('FLAG_ROW')
+            crosscorrelations = numpy.squeeze(table.getcol('DATA'))
+
+            n_reference_stations = len(reference_stations)
+
+            antenna1 = list(table.getcol('ANTENNA1'))[:n_reference_stations]
+            antenna2 = list(table.getcol('ANTENNA2'))[:n_reference_stations]
+
+
+
+            timestamps = timestamps[::n_reference_stations]
+            n_timestamps = len(timestamps)
+            n_polarizations = crosscorrelations.shape[-1]
+
+            ordered_crosscorrelations = crosscorrelations.reshape([n_reference_stations,
+                                                                   n_timestamps,
+                                                                   n_polarizations], order='F')
+            flags = flags.reshape([n_reference_stations, n_timestamps], order='F')
+
         finally:
             data_table.close()
             antennas_table.close()
 
-
-        return (timestamp[selected_data], cross_correlation[selected_data])
+        return (crosscorrelations, timestamps, flags)
 
     def __repr__(self):
         return 'MeasurementSet(%d) located in %s for sas_id %d and sub_band_id %d' % (id(self),

CAL/CalibrationCommon/lib/utils.py

@@ -3,17 +3,6 @@ import os
 from glob import glob
 
 
-def list_bsf_files_in_path(path):
-    bsf_files_name_pattern = 'Holog-*.txt'
-    bsf_files_path_pattern = os.path.join(path, bsf_files_name_pattern)
-    matched_path_list = glob(bsf_files_path_pattern)
-    matched_file_path_list = list(filter(lambda path_i : os.path.isfile(path_i), matched_path_list))
-    matched_file_name_list = list(map(lambda path_i: os.path.basename(path_i),
-                                      matched_file_path_list))
-    holography_beam_specification = HolographySpecification.create_hs_list_from_name_list_and_path(
-        matched_file_name_list, path)
-    return holography_beam_specification
-
 
 def is_observation_in_range(start, end, from_datetime, to_datetime):
     """
@@ -35,7 +24,7 @@ def is_observation_in_range(start, end, from_datetime, to_datetime):
 
 
 def match_holography_beam_specification_file_with_observation(path):
-    bsf_files = list_bsf_files_in_path(path)
+    bsf_files = HolographySpecification.list_bsf_files_in_path(path)
     observation_list = HolographyObservation.list_observations_in_path(path)
     matched_observation_bsf_pair = []
     for bsf_file in bsf_files:

CAL/CalibrationCommon/test/CMakeLists.txt

@@ -22,6 +22,6 @@ include(LofarCTest)
 lofar_add_test(t_holography_ms_class)
 lofar_add_test(t_holography_beam_specification_class)
 lofar_add_test(t_holography_dataset_class)
-
+lofar_add_test(t_holography_observation)
 
 

CAL/CalibrationCommon/test/t_holography_dataset_class.py

 import logging
 import unittest
-from lofar.calibration.common.datacontainers import HolographyDataset
+from lofar.calibration.common.datacontainers import HolographyDataset, HolographySpecification
+from lofar.calibration.common.datacontainers import HolographyObservation
 
 logger = logging.getLogger('t_holography_dataset_class')
+path_to_test_data = '/data/test/HolographyObservation'
+
 
 class TestHolographyDatasetClass(unittest.TestCase):
     def test_one(self):
-        self.assertTrue(True)
\ No newline at end of file
+        holist = HolographyObservation.list_observations_in_path(path_to_test_data)
+        hbsflist = HolographySpecification.list_bsf_files_in_path(path_to_test_data)
+        for hbsf in hbsflist:
+            hbsf.read_file()
+        ho_per_ms = [(hbsf, ho) for hbsf, ho in zip(hbsflist, holist)]
+
+        holography_dataset = HolographyDataset()
+
+        holography_dataset.load_from_beam_specification_and_ms('CS013HBA0', ho_per_ms)
+        self.assertEqual(holography_dataset.source_name, '3C 147')
+
+if __name__ == '__main__':
+    logging.basicConfig(format='%(name)s : %(message)s')
+    unittest.main()
\ No newline at end of file

CAL/CalibrationCommon/test/t_holography_observation.py

+import logging
+import unittest
+from lofar.calibration.common.datacontainers import HolographyObservation
+
+logger = logging.getLogger('t_holography_dataset_class')
+path_to_test_data = '/data/test/HolographyObservation'
+
+
+class TestHolographyDatasetClass(unittest.TestCase):
+    def test_holography_observation_path_listing(self):
+        list = HolographyObservation.list_observations_in_path(path_to_test_data)
+
+        ## Assert correct reading of the files
+        self.assertEqual(len(list), 1)
+        HolographyObservation_freq_one = list.pop()
+        self.assertEqual(HolographyObservation_freq_one.sas_id, 661022)
+        self.assertEqual(HolographyObservation_freq_one.source_name, '3C 147')
+        self.assertEqual(HolographyObservation_freq_one.sub_band, 76)
+        self.assertEqual(HolographyObservation_freq_one.start_mjd, 5038710120.0)
+        self.assertEqual(HolographyObservation_freq_one.end_mjd, 5038710149.540695)
+        self.assertAlmostEqual((HolographyObservation_freq_one.end_datetime -
+                         HolographyObservation_freq_one.start_datetime).seconds,
+                               29)
+
+if __name__ == '__main__':
+    unittest.main()
\ No newline at end of file

CAL/CalibrationCommon/test/t_holography_observation.run

+#!/bin/bash
+
+# Copyright (C) 2018    ASTRON (Netherlands Institute for Radio Astronomy)
+# P.O. Box 2, 7990 AA Dwingeloo, The Netherlands
+#
+# This file is part of the LOFAR software suite.
+# The LOFAR software suite is free software: you can redistribute it and/or
+# modify it under the terms of the GNU General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# The LOFAR software suite is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.    See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License along
+# with the LOFAR software suite. If not, see <http://www.gnu.org/licenses/>.
+
+# Run the unit test
+source python-coverage.sh
+python_coverage_test "*CAL*" t_holography_observation.py

CAL/CalibrationCommon/test/t_holography_observation.sh

+#!/bin/sh
+
+# Copyright (C) 2018    ASTRON (Netherlands Institute for Radio Astronomy)
+# P.O. Box 2, 7990 AA Dwingeloo, The Netherlands
+#
+# This file is part of the LOFAR software suite.
+# The LOFAR software suite is free software: you can redistribute it and/or
+# modify it under the terms of the GNU General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# The LOFAR software suite is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.    See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License along
+# with the LOFAR software suite. If not, see <http://www.gnu.org/licenses/>.
+
+./runctest.sh t_holography_observation
\ No newline at end of file