-
Thomas Jürges authoredThomas Jürges authored
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
holography_dataset.py 16.03 KiB
from .holography_specification import HolographySpecification
from .holography_dataset_definitions import *
from lofar.calibration.common.datacontainers.holography_observation import HolographyObservation
import logging
import numpy
import h5py
logger = logging.getLogger(__name__)
class HolographyDataset():
def __init__(self):
self.rcu_list = set() # array of ints
self.mode = None # integer
self.sas_ids = set() # array of strings
self.target_station_name = None # string
self.target_station_position = None # list of 3 floats
self.source_name = None # string
self.source_position = None # list of 3 floats
self.start_time = None # date time when the observation started in MJD in seconds (float)
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 = 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
# 4 polarization crosscorrelation for the 4 polarizations, the l and m coordinates, and
# the timestamp in mjd of the sample, and whether or not the data has been flagged
# numpy.dtype([('XX', numpy.float),
# ('YY', numpy.float),
# ('XY', numpy.float),
# ('YX', numpy.float),
# ('l', numpy.float),
# ('m', numpy.float),
# ('t', numpy.float),
# ('flag', numpy.bool)]
# )
def load_from_beam_specification_and_ms(self, station_name, list_of_hbs_ms_tuples):
"""
Loads the dataset from the specification files and the measurements for the given station
name
:param station_name: target station name
:param hb_specifications: a list containing (hbs, ms) per frequency
"""
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:
"""
self.data = dict()
self.flags = 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:
reference_station_names, cross_correlation =\
ho.ms_for_a_given_beamlet_number[beamlet].\
read_cross_correlation_time_flags_per_station_names(station_name,
self.reference_stations)
print(reference_station_names)
for reference_station_index, reference_station in\
enumerate(reference_station_names):
if reference_station not in self.data:
self.data[reference_station] = dict()
if reference_station not in self.flags:
self.flags[reference_station] = dict()
if frequency not in self.data[reference_station]:
self.data[reference_station][frequency] = dict()
if frequency not in self.flags[reference_station]:
self.flags[reference_station][frequency] = dict()
self.data[reference_station][frequency][beamlet] = \
cross_correlation[reference_station_index, :]
def __collect_preliminary_information(self, station_name, list_of_hbs_ho_tuples):
"""
This routines reads both the holography beam specifications files and the holography
observation to gather the list of rcus, the mode, the target station name and position,
the source name and position, the start and the end time, the rotation matrix to convert
from ra and dec to l and m, the antenna field positions, the list of the reference
stations, the frequencies, the ra and dec at which the beams point at.
All this data is essential to interpret the recorded holography beams cross correlations
:param list_of_hbs_ho_tuples: a list containing (hbs, ho) per frequency
:type list_of_hbs_ho_tuples: list[(HolographySpecification, HolographyObservation)]
"""
mode = set()
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)
if station_name in hbs.target_station_names:
beam_specifications = hbs.get_beam_specifications_per_station_name(station_name)
for beam_specification in beam_specifications:
self.rcu_list.update(beam_specification.rcus_involved)
mode.add(beam_specification.rcus_mode)
source_name.add(ho.source_name)
source_position.add(
(beam_specification.station_pointing['ra'],
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.sas_ids.add(ho.sas_id)
self.target_station_name = station_name
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)
self.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):
self.ra_dec[frequency_index, beamlet_index, :] = \
virtual_pointing[(frequency, beamlet)]
# 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]
first_ms = first_holography_observation.ms_for_a_given_beamlet_number.values()[0]
station_position, tile_offset, axes_coordinates = first_ms.\
get_station_position_tile_offsets_and_axes_coordinate_for_station_name(
station_name)
self.antenna_field_position = [list(station_position - antenna_offset)
for antenna_offset in tile_offset]
self.target_station_position = list(station_position)
self.rotation_matrix = axes_coordinates
if station_name not in target_stations:
logger.error('Station %s was not involved in the observation.'
' The target stations for this observation are %s',
station_name, target_stations)
raise Exception('Station %s was not involved in the observation.'
% station_name,)
if len(mode) > 1:
raise ValueError('Multiple RCUs mode are not supported')
else:
self.mode = mode.pop()
if len(source_position) > 1:
logger.error('Multiple source positions are not supported: %s', source_position)
raise ValueError('Multiple source positions are not supported')
else:
self.source_position = source_position.pop()
if len(source_name) > 1:
raise ValueError('Multiple source name are not supported')
else:
self.source_name = source_name.pop()
@staticmethod
def load_from_file(path):
"""
It reads a holography dataset from an HDF5 file and returns a
HolographyDataset class
:param path: path to file
:return: the read dataset
"""
f = h5py.File(path, "r")
result = HolographyDataset()
result.version = f.attrs[HDS_VERSION]
result.mode= f.attrs[HDS_MODE]
result.rcu_list = f.attrs[HDS_RCU_LIST]
result.sas_ids = f.attrs[HDS_SAS_ID]
result.target_station_name = f.attrs[HDS_TARGET_STATION_NAME]
result.target_station_position = f.attrs[HDS_TARGET_STATION_POSITION]
result.source_name = f.attrs[HDS_SOURCE_NAME]
result.source_position = f.attrs[HDS_SOURCE_POSITION]
(result.start_time, result.end_time) = f.attrs[HDS_OBSERVATION_TIME]
result.rotation_matrix = f.attrs[HDS_ROTATION_MATRIX]
result.antenna_field_position = f.attrs[HDS_ANTENNA_FIELD_POSITION]
result.reference_stations = f[HDS_REFERENCE_STATION]
result.frequencies = f[HDS_FREQUENCY]
result.ra_dec = f[HDS_RA_DEC]
data = f[HDS_DATA]
result.beamlets = [i for i in range(data.shape[2])]
for reference_station_index, reference_station in enumerate(result.reference_stations):
for frequency_index, frequency in result.frequencies:
for beamlet in result.beamlets:
result.data[reference_station][frequency][beamlet]["XX"] =\
data[reference_station_index, frequency_index, beamlet]["XX"]
result.data[reference_station][frequency][beamlet]["XY"] =\
data[reference_station_index, frequency_index, beamlet]["XY"]
result.data[reference_station][frequency][beamlet]["YX"] =\
data[reference_station_index, frequency_index, beamlet]["YX"]
result.data[reference_station][frequency][beamlet]["YY"] =\
data[reference_station_index, frequency_index, beamlet]["YY"]
result.data[reference_station][frequency][beamlet]["t"] =\
data[reference_station_index, frequency_index, beamlet]["t"]
result.data[reference_station][frequency][beamlet]["l"] =\
data[reference_station_index, frequency_index, beamlet]["l"]
result.data[reference_station][frequency][beamlet]["m"] =\
data[reference_station_index, frequency_index, beamlet]["m"]
result.data[reference_station][frequency][beamlet]["flag"] =\
data[reference_station_index, frequency_index, beamlet]["flag"]
result.data = data
f.close()
return result
def store_to_file(self, path):
"""
Stores the holography dataset at the given path
:param path: path to file
"""
# Prepare the HDF5 data structs.
f = h5py.File(path, "w")
try:
# Create the ATTRS
# Start with the version information
f.attrs[HDS_VERSION] = HOLOGRAPHY_DATA_SET_VERSION
# Create a numpy array that contains the RCU list as unicode strings.
# The type has to be str because nump/h5py cannot handle unicode string
# arrays.
f.attrs[HDS_RCU_LIST] = numpy.array(self.rcu_list,
dtype=h5py.special_dtype(vlen = str))
# RCU mode
f.attrs[HDS_MODE] = self.mode
# Moan... Again this needs to be stored like that.
f.attrs[HDS_SAS_ID] = numpy.array(self.sas_ids,
dtype=h5py.special_dtype(vlen = str))
f.attrs[HDS_TARGET_STATION_NAME] = self.target_station_name
f.attrs[HDS_TARGET_STATION_POSITION] = self.target_station_position
f.attrs[HDS_SOURCE_NAME] = self.source_name
f.attrs[HDS_SOURCE_POSITION] = self.source_position
f.attrs[HDS_OBSERVATION_TIME] = numpy.array([self.start_time, self.end_time])
f.attrs[HDS_ROTATION_MATRIX] = self.rotation_matrix
f.attrs[HDS_ANTENNA_FIELD_POSITION] = self.antenna_field_position
# Data tables
f[HDS_REFERENCE_STATION] = numpy.array(self.reference_stations,
dtype = h5py.special_dtype(vlen = str))
f[HDS_FREQUENCY] = self.frequencies
ra_dec_type = numpy.dtype([('RA', numpy.float64),
('DEC', numpy.float64)])
f[HDS_RA_DEC] = numpy.array(self.ra_dec, dtype = ra_dec_type)
dataset = f.create_dataset("data",
(len(self.reference_stations),
len(self.frequencies),
len(self.beamlets)),
dtype = h5py.special_dtype(vlen = data_sample_type))
for reference_station_index, reference_station in enumerate(self.reference_stations):
for frequency_index, frequency in enumerate(self.frequencies):
for beamlet in self.beamlets:
dataset[reference_station_index, frequency_index, beamlet]["XX"] = \
self.data[reference_station][frequency][beamlet]["XX"]
dataset[reference_station_index, frequency_index, beamlet]["XY"] = \
self.data[reference_station][frequency][beamlet]["XY"]
dataset[reference_station_index, frequency_index, beamlet]["YX"] = \
self.data[reference_station][frequency][beamlet]["YX"]
dataset[reference_station_index, frequency_index, beamlet]["YY"] = \
self.data[reference_station][frequency][beamlet]["YY"]
dataset[reference_station_index, frequency_index, beamlet]["t"] = \
self.data[reference_station][frequency][beamlet]["t"]
dataset[reference_station_index, frequency_index, beamlet]["l"] = \
self.data[reference_station][frequency][beamlet]["l"]
dataset[reference_station_index, frequency_index, beamlet]["m"] = \
self.data[reference_station][frequency][beamlet]["m"]
dataset[reference_station_index, frequency_index, beamlet]["flag"] = \
self.data[reference_station][frequency][beamlet]["flag"]
finally:
f.close()