Commit 4c437262 authored by Philipp Arras's avatar Philipp Arras
Browse files

Add demo with ms interface

parent 44103b6f
# This program 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.
#
# This program 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 this program. If not, see <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2019 Max-Planck-Society
from os.path import join
from time import time
import matplotlib.pyplot as plt
import nifty_gridder as ng
import numpy as np
from casacore.tables import table
# Assumptions:
# - Only one field
# - Only one spectral window
# - Flag both LL and RR if one is flagged
ms = 'supernovashell.55.7+3.4.spw0.ms'
t = table(ms, readonly=True)
uvw = t.getcol("UVW") # [:, uvw]
vis = np.array(t.getcol("DATA"), dtype=np.complex128) # [:, ch, corr]
wgt = t.getcol("WEIGHT")
# Flag if one correlation is flagged
flags = np.any(np.array(t.getcol('FLAG'), np.bool), axis=2) # [:, ch]
if len(set(t.getcol('FIELD_ID'))) != 1:
raise RuntimeError
if len(set(t.getcol('DATA_DESC_ID'))) != 1:
raise RuntimeError
t.close()
print('# Rows: {}'.format(vis.shape[0]))
print('# Channels: {}'.format(vis.shape[1]))
print('# Correlations: {}'.format(vis.shape[2]))
print("{} % flagged".format(np.sum(flags)/flags.size*100))
t = table(join(ms, 'SPECTRAL_WINDOW'), readonly=True)
freq = t.getcol('CHAN_FREQ')[0]
t.close()
# Select either RR+LL or XX+YY
t = table(join(ms, 'POLARIZATION'), readonly=True)
pol = list(t.getcol('CORR_TYPE')[0])
t.close()
if set(pol) <= set([5, 6, 7, 8]):
ind = [pol.index(5), pol.index(8)]
else:
ind = [pol.index(9), pol.index(12)]
vis = np.sum(vis[:, :, ind], axis=2)
wgt = 1/np.sum(1/wgt, axis=1)
# WEIGHT -> WEIGHT_SPECTRUM
wgt = np.repeat(wgt[:, None], len(freq), axis=1)
wgt[flags] = 0
npixdirty = 4096
DEG2RAD = np.pi/180
pixsize = 2.3/npixdirty*DEG2RAD
nthreads = 4
epsilon = 1e-4
t0 = time()
print('Start gridding...')
dirty = ng.full_gridding(uvw, freq, vis, wgt, npixdirty, npixdirty, pixsize,
pixsize, epsilon, nthreads)
print('Done')
t = time() - t0
print("{} s".format(t))
print("{} visibilities/nthreads/s".format(np.sum(wgt != 0)/nthreads/t))
plt.imshow(dirty.real)
plt.show()
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment