Skip to content

Easy spreadsheet copy and paste #104

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
lmagnus wants to merge 2 commits into
base: master
Choose a base branch
from
+36 −29
Open

Easy spreadsheet copy and paste #104

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
48d8973
Added a printing section to the end of the baseline cal script to all…
lmagnus Apr 21, 2016
de74a3f
attempt to get a spreadhseet formatable output from bl cal
lmagnus Jun 10, 2016
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
65 changes: 36 additions & 29 deletions reduction/baseline_cal_reduction.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,8 @@
#

import optparse

import warnings
warnings.filterwarnings('ignore')
import numpy as np
import matplotlib.pyplot as plt

Expand Down Expand Up @@ -53,7 +54,7 @@

katpoint.logger.setLevel(30)

print "\nLoading and processing data...\n"
#print "\nLoading and processing data...\n"
data = katdal.open(args, ref_ant=opts.ref_ant, time_offset=opts.time_offset)
center_freq = data.freqs[data.shape[1] // 2]

Expand Down Expand Up @@ -100,8 +101,8 @@
max_sigma_delay = delay_period / np.sqrt(12) / np.sqrt(num_chans - 1)

ant_list = ', '.join([(ant.name + ' (*ref*)' if ind == ref_ant_ind else ant.name) for ind, ant in enumerate(data.ants)])
print 'antennas (%d): %s [pol %s]' % (len(data.ants), ant_list, opts.pol)
print 'baselines (%d): %s' % (num_bls, ' '.join([('%d-%d' % (indA, indB)) for indA, indB in baseline_inds]))
#print 'antennas (%d): %s [pol %s]' % (len(data.ants), ant_list, opts.pol)
#print 'baselines (%d): %s' % (num_bls, ' '.join([('%d-%d' % (indA, indB)) for indA, indB in baseline_inds]))

# Create an extended baseline difference matrix mapping antenna parameters to baselines
# This array has shape (B P, 4), where B is number of baselines and P = 4 A is total number of parameters
Expand All @@ -120,13 +121,13 @@
for scan_ind, state, target in data.scans():
num_ts = data.shape[0]
if state != 'track':
print "scan %3d (%4d samples) skipped '%s'" % (scan_ind, num_ts, state)
#print "scan %3d (%4d samples) skipped '%s'" % (scan_ind, num_ts, state)
continue
if num_ts < 2:
print "scan %3d (%4d samples) skipped - too short" % (scan_ind, num_ts)
#print "scan %3d (%4d samples) skipped - too short" % (scan_ind, num_ts)
continue
if target.name in excluded_targets:
print "scan %3d (%4d samples) skipped - excluded '%s'" % (scan_ind, num_ts, target.name)
#print "scan %3d (%4d samples) skipped - excluded '%s'" % (scan_ind, num_ts, target.name)
continue
# Extract visibilities for scan as an array of shape (T, F, B)
vis = data.vis[:]
Expand Down Expand Up @@ -170,8 +171,8 @@
# Rearrange vis phase to have shape (B F, T), which will form one column in fringe plot
scan_phase.append(np.angle(vis).T.reshape(-1, num_ts))
scan_rel_sigma_delay = delay_stats_sigma.mean(axis=0) / max_sigma_delay
print "scan %3d (%4d samples) %s '%s'" % \
(scan_ind, num_ts, ' '.join([('%.3f' % rel_sigma) for rel_sigma in scan_rel_sigma_delay]), target.name)
#print "scan %3d (%4d samples) %s '%s'" % \
# (scan_ind, num_ts, ' '.join([('%.3f' % rel_sigma) for rel_sigma in scan_rel_sigma_delay]), target.name)
if not augmented_targetdir:
raise RuntimeError('No usable scans found (are you tracking any targets?)')
# Concatenate per-baseline arrays into a single array for data set
Expand Down Expand Up @@ -199,15 +200,15 @@
good = sigma_delay < opts.max_sigma * max_sigma_delay
A = A[:, good]
b = b[good]
print '\nFitting %d parameters to %d data points (discarded %d)...' % (A.shape[0], len(b), len(sigma_delay) - len(b))
#print '\nFitting %d parameters to %d data points (discarded %d)...' % (A.shape[0], len(b), len(sigma_delay) - len(b))
if len(b) == 0:
raise ValueError('No solution possible, as all data points were discarded')
# Solve linear least-squares problem using SVD (see NRinC, 2nd ed, Eq. 15.4.17)
U, s, Vrt = np.linalg.svd(A.transpose(), full_matrices=False)
params = np.dot(Vrt.T, np.dot(U.T, b) / s)
# Also obtain standard errors of parameters (see NRinC, 2nd ed, Eq. 15.4.19)
sigma_params = np.sqrt(np.sum((Vrt.T / s[np.newaxis, :]) ** 2, axis=1))
print 'Condition number = %.3f' % (s[0] / s[-1],)
#print 'Condition number = %.3f' % (s[0] / s[-1],)

# Reshape parameters to be per antenna, also inserting a row of zeros for reference antenna
ant_params = params.reshape(-1, 4)
Expand All @@ -232,24 +233,30 @@
new_resid = unwrapped_group_delay - new_predicted_delay

# Output results
#for n, ant in enumerate(data.ants):
# print "\nAntenna", ant.name, ' (*REFERENCE*)' if n == ref_ant_ind else ''
# print "------------"
# print "E (m): %7.3f +- %.5f (was %7.3f)%s" % \
# (positions[n, 0], sigma_positions[n, 0], old_positions[n, 0],
# ' *' if np.abs(positions[n, 0] - old_positions[n, 0]) > 3 * sigma_positions[n, 0] else '')
# print "N (m): %7.3f +- %.5f (was %7.3f)%s" % \
# (positions[n, 1], sigma_positions[n, 1], old_positions[n, 1],
# ' *' if np.abs(positions[n, 1] - old_positions[n, 1]) > 3 * sigma_positions[n, 1] else '')
# print "U (m): %7.3f +- %.5f (was %7.3f)%s" % \
# (positions[n, 2], sigma_positions[n, 2], old_positions[n, 2],
# ' *' if np.abs(positions[n, 2] - old_positions[n, 2]) > 3 * sigma_positions[n, 2] else '')
# print "Cable length (m): %7.3f +- %.5f (was %7.3f)%s" % \
# (cable_lengths[n], sigma_cable_lengths[n], old_cable_lengths[n],
# ' *' if np.abs(cable_lengths[n] - old_cable_lengths[n]) > 3 * sigma_cable_lengths[n] else '')
# print "Receiver delay (ns): %7.3f +- %.3f (was %7.3f)%s" % \
# (receiver_delays[n] * 1e9, sigma_receiver_delays[n] * 1e9, old_receiver_delays[n] * 1e9,
# ' *' if np.abs(receiver_delays[n] - old_receiver_delays[n]) > 3 * sigma_receiver_delays[n] else '')

# print lines for easy incoporation into a spreadsheet
#print "\nFor spreadsheet logging"
for n, ant in enumerate(data.ants):
print "\nAntenna", ant.name, ' (*REFERENCE*)' if n == ref_ant_ind else ''
print "------------"
print "E (m): %7.3f +- %.5f (was %7.3f)%s" % \
(positions[n, 0], sigma_positions[n, 0], old_positions[n, 0],
' *' if np.abs(positions[n, 0] - old_positions[n, 0]) > 3 * sigma_positions[n, 0] else '')
print "N (m): %7.3f +- %.5f (was %7.3f)%s" % \
(positions[n, 1], sigma_positions[n, 1], old_positions[n, 1],
' *' if np.abs(positions[n, 1] - old_positions[n, 1]) > 3 * sigma_positions[n, 1] else '')
print "U (m): %7.3f +- %.5f (was %7.3f)%s" % \
(positions[n, 2], sigma_positions[n, 2], old_positions[n, 2],
' *' if np.abs(positions[n, 2] - old_positions[n, 2]) > 3 * sigma_positions[n, 2] else '')
print "Cable length (m): %7.3f +- %.5f (was %7.3f)%s" % \
(cable_lengths[n], sigma_cable_lengths[n], old_cable_lengths[n],
' *' if np.abs(cable_lengths[n] - old_cable_lengths[n]) > 3 * sigma_cable_lengths[n] else '')
print "Receiver delay (ns): %7.3f +- %.3f (was %7.3f)%s" % \
(receiver_delays[n] * 1e9, sigma_receiver_delays[n] * 1e9, old_receiver_delays[n] * 1e9,
' *' if np.abs(receiver_delays[n] - old_receiver_delays[n]) > 3 * sigma_receiver_delays[n] else '')
print "%s\t%s\t%s\t%7.3f\t%7.3f\t%7.3f\t%7.3f\t%7.3f\t%7.3f\t%s"%(data.name.split('/')[-1].split('.')[0], \
Copy link
Copy Markdown
Contributor

@ludwigschwardt ludwigschwardt Apr 22, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The base part of the filename might be better extracted via os.path.splitext(os.path.basename('path'))[0] but that is a minor concern.

opts.pol,ant.name,positions[n, 0],positions[n, 1],positions[n, 2],cable_lengths[n],receiver_delays[n] * 1e9,s[0] / s[-1],data.ref_ant)

scan_lengths = [len(ts) for ts in scan_timestamps]
scan_bl_starts = num_bls * np.cumsum([0] + scan_lengths)[:-1]
Expand Down Expand Up @@ -343,4 +350,4 @@ def extract_scan_segments(x):
ax.set_yticklabels([])
plt.title('Antenna positions and source directions')

plt.show()
#plt.show()