Skip to content
Snippets Groups Projects
Commit 585442fb authored by Mattia Mancini's avatar Mattia Mancini
Browse files

SSB-42: dataset class work in progress

parent 1b734370
No related branches found
No related tags found
1 merge request!44Merge back holography to master
......@@ -2,7 +2,7 @@
Authors: [Thomas Jürges](mailto:jurges@astron.nl), [Mattia Mancini](mailto:mancini@astron.nl)
HOLOGRAPHY_DATA_SET_VERSION
Date: 2018-09-10
......
......@@ -33,6 +33,8 @@ class HolographyDataset():
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
self.flags = None # array(NReferenceStations, Nfrequencies, Nbeamlet) that contains the
# an array of boolean that define if a data in time has been flagged
def load_from_beam_specification_and_ms(self, station_name, list_of_hbs_ms_tuples):
"""
......@@ -55,24 +57,26 @@ class HolographyDataset():
:return:
"""
n_frequencies = len(self.frequencies)
n_beamlets = len(self.beamlets)
n_reference_stations = len(self.reference_stations)
data = defaultdict()
self.data = defaultdict(defaultdict)
self.flags = defaultdict(defaultdict)
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].\
cross_correlation, reference_station_names, 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, )]
for reference_station_index, reference_station in\
enumerate(reference_station_names):
print(reference_station_index)
self.data[reference_station][frequency][beamlet] = \
cross_correlation[reference_station_index, :, :]
self.flags[reference_station][frequency][beamlet] = \
flags[reference_station_index, :]
def __collect_preliminary_information(self, station_name, list_of_hbs_ho_tuples):
"""
......@@ -201,26 +205,26 @@ class HolographyDataset():
f = h5py.File(path, "r")
result = HolographyDataset()
result.version = f.attr["Version"]
result.version = f.attrs["Version"]
result.mode= f.attr["mode"]
result.mode= f.attrs["mode"]
result.rcu_list = f.attr["RCU list"]
result.rcu_list = f.attrs["RCU list"]
result.sas_ids = f.attr["SAS id"]
result.sas_ids = f.attrs["SAS id"]
result.target_station_name = f.attr["Target station name"]
result.target_station_name = f.attrs["Target station name"]
result.target_station_position = f.attr["Target station position"]
result.target_station_position = f.attrs["Target station position"]
result.source_name = f.attr["Source name"]
result.source_name = f.attrs["Source name"]
result.source_position = f.attr["Source position"]
(result.start_time, result.end_time) = f.attr["Observation time"]
result.source_position = f.attrs["Source position"]
(result.start_time, result.end_time) = f.attrs["Observation time"]
result.rotation_matrix = f.attr["Rotation matrix"]
result.rotation_matrix = f.attrs["Rotation matrix"]
result.antenna_field_position = f.attr["Antenna field position"]
result.antenna_field_position = f.attrs["Antenna field position"]
result.reference_stations = f["Reference station"]
......@@ -229,17 +233,28 @@ class HolographyDataset():
result.ra_dec = f["RA DEC"]
data = f["data"]
for reference_station in result.reference_stations:
for frequency in result.frequencies:
for beamlet in self.beamlets:
result.data[reference_station][frequency][beamlet]["XX"] = data[reference_station][frequency][beamlet]["XX"]
result.data[reference_station][frequency][beamlet]["XY"] = data[reference_station][frequency][beamlet]["XY"]
result.data[reference_station][frequency][beamlet]["YX"] = data[reference_station][frequency][beamlet]["YX"]
result.data[reference_station][frequency][beamlet]["YY"] = data[reference_station][frequency][beamlet]["YY"]
result.data[reference_station][frequency][beamlet]["t"] = data[reference_station][frequency][beamlet]["t"]
result.data[reference_station][frequency][beamlet]["l"] = data[reference_station][frequency][beamlet]["l"]
result.data[reference_station][frequency][beamlet]["m"] = data[reference_station][frequency][beamlet]["m"]
result.data[reference_station][frequency][beamlet]["flag"] = data[reference_station][frequency][beamlet]["flag"]
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()
......@@ -254,64 +269,80 @@ class HolographyDataset():
# Prepare the HDF5 data structs.
f = h5py.File(path, "w")
# Create the ATTRS
# Start with the version information
f.attrs["Version"] = self.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.attr["RCU list"] = numpy.array(self.rcu_list, dtype = h5py.special_dtype(vlen = str))
# RCU mode
f.attr["mode"] = self.mode
# Moan... Again this needs to be stored like that.
f.attr["SAS id"] = numpy.array(self.sas_ids, dtype = h5py.special_dtype(vlen = str))
f.attr["Target station name"] = self.target_station_name
f.attr["Target station position"] = self.target_station_position
f.attr["Source name"] = self.source_name
f.attr["Source position"] = self.source_position
f.attr["Observation time"] = numpy.array([self.start_time, self.end_time])
f.attr["Rotation matrix"] = self.rotation_matrix
f.attr["Antenna field position"] = self.antenna_field_position
# Data tables
f["Reference station"] = numpy.array(self.reference_stations, dtype = h5py.special_dtype(vlen = str))
f["frequency"] = self.frequencies
f["RA DEC"] = self.ra_dec
sample_type = numpy.dtype([
('XX', numpy.complex64),
('XY', numpy.complex64),
('YX', numpy.complex64),
('YY', numpy.complex64),
('t', numpy.float64),
('l', numpy.float64),
('m', numpy.float64),
('flag', numpy.bool)])
dataset = f.create_dataset("data", (len(self.reference_stations), len(self.frequencies), len(self.beamlets)), dtype = h5py.special_dtype(vlen = sample_type))
for reference_station in self.reference_stations:
for frequency in self.frequencies:
for beamlet in self.beamlets:
dataset[reference_station][frequency][beamlet]["XX"] = self.data[reference_station][frequency][beamlet]["XX"]
dataset[reference_station][frequency][beamlet]["XY"] = self.data[reference_station][frequency][beamlet]["XY"]
dataset[reference_station][frequency][beamlet]["YX"] = self.data[reference_station][frequency][beamlet]["YX"]
dataset[reference_station][frequency][beamlet]["YY"] = self.data[reference_station][frequency][beamlet]["YY"]
dataset[reference_station][frequency][beamlet]["t"] = self.data[reference_station][frequency][beamlet]["t"]
dataset[reference_station][frequency][beamlet]["l"] = self.data[reference_station][frequency][beamlet]["l"]
dataset[reference_station][frequency][beamlet]["m"] = self.data[reference_station][frequency][beamlet]["m"]
dataset[reference_station][frequency][beamlet]["flag"] = self.data[reference_station][frequency][beamlet]["flag"]
f.close()
try:
# Create the ATTRS
# Start with the version information
f.attrs["Version"] = self.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["RCU list"] = numpy.array(self.rcu_list,
dtype=h5py.special_dtype(vlen = str))
# RCU mode
f.attrs["mode"] = self.mode
# Moan... Again this needs to be stored like that.
f.attrs["SAS id"] = numpy.array(self.sas_ids,
dtype=h5py.special_dtype(vlen = str))
f.attrs["Target station name"] = self.target_station_name
f.attrs["Target station position"] = self.target_station_position
f.attrs["Source name"] = self.source_name
f.attrs["Source position"] = self.source_position
f.attrs["Observation time"] = numpy.array([self.start_time, self.end_time])
f.attrs["Rotation matrix"] = self.rotation_matrix
f.attrs["Antenna field position"] = self.antenna_field_position
# Data tables
f["Reference station"] = numpy.array(self.reference_stations,
dtype=h5py.special_dtype(vlen = str))
f["frequency"] = self.frequencies
f["RA DEC"] = self.ra_dec
sample_type = numpy.dtype([
('XX', numpy.complex64),
('XY', numpy.complex64),
('YX', numpy.complex64),
('YY', numpy.complex64),
('t', numpy.float64),
('l', numpy.float64),
('m', numpy.float64),
('flag', numpy.bool)])
dataset = f.create_dataset("data",
(len(self.reference_stations),
len(self.frequencies),
len(self.beamlets)),
dtype = h5py.special_dtype(vlen = 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()
......@@ -137,7 +137,7 @@ class HolographyMeasurementSet(object):
antenna_name_index = antennas_table.index('NAME')
baseline_selection = ','.join(reference_stations)
baseline_selection += '&&' + target_station
baseline_selection += '&' + target_station
TAQL_query_syntax = 'mscal.baseline($baseline_selection)'
......@@ -150,32 +150,34 @@ class HolographyMeasurementSet(object):
'mscal.ant2name() as ANTENNA_NAME2',
sortlist='TIME, ANTENNA1, ANTENNA2')
print(table.getcol('ANTENNA_NAME1'))
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]
antenna1 = table.getcol('ANTENNA1')
antenna2 = table.getcol('ANTENNA2')
reference_stations_names = [a2 if a1 == target_station else a1
for a1, a2 in zip(antenna1, antenna2)]
n_reference_stations = len(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')
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 (crosscorrelations, timestamps, flags)
return (crosscorrelations, reference_stations_names, timestamps, flags)
def __repr__(self):
return 'MeasurementSet(%d) located in %s for sas_id %d and sub_band_id %d' % (id(self),
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment