-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathPyra_grid.py
115 lines (96 loc) · 3.74 KB
/
Pyra_grid.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
from pyralysis.io import DaskMS
import numpy as np
import astropy.units as u
from pyralysis.transformers import Gridder, HermitianSymmetry, DirtyMapper
from pyralysis.io import FITS
from astropy.units import Quantity
def completeuvplane_complex(V_half):
V_half = np.squeeze(V_half)
nu, nv = V_half.shape
Nside = max(nu, nv)
V_full = np.zeros((Nside, Nside)).astype(complex)
V_full[:, int(Nside // 2):] = V_half[:, 0:int(Nside // 2)]
Vmatrix = np.matrix(V_full)
V_H = Vmatrix.getH()
V_H = np.asarray(V_H)
V_H = np.flip(V_H, 0)
V_H = np.rot90(V_H)
V_full[:, 0:int(Nside // 2)] = V_H[:, 0:int(Nside // 2)]
#V_full_R = V_full.real
#V_full_I = V_full.imag
#V_full_amp = V_full_R**2 + V_full_I**2
#print(V_full_amp.shape)
#Vtools.View(V_full_amp)
return V_full
def completeuvplane(V_half):
V_half = np.squeeze(V_half)
nu, nv = V_half.shape
Nside = max(nu, nv)
V_full = np.zeros((Nside, Nside))
V_full[:, int(Nside // 2):] = V_half[:, 0:int(Nside // 2)]
V_H = V_full.transpose()
V_H = np.flip(V_H, 0)
V_H = np.rot90(V_H)
V_full[:, 0:int(Nside // 2)] = V_H[:, 0:int(Nside // 2)]
return V_full
def gridvis(file_ms,
imsize=2048,
data_column=None,
hermitian_symmetry=True,
dx=None,
wantpsf=False,
wantdirtymap=False):
print("processing: ", file_ms)
x = DaskMS(input_name=file_ms, chunks={'row': 1E8, 'chan': 30})
dataset = x.read(filter_flag_column=True, calculate_psf=True)
#dataset.field.mean_ref_dir
#dataset.psf[0].sigma
print("baselines in klambdas:",
dataset.min_baseline / dataset.spws.lambda_max,
dataset.max_baseline / dataset.spws.lambda_min)
print("baselines in m:", dataset.min_baseline, dataset.max_baseline)
if hermitian_symmetry:
h_symmetry = HermitianSymmetry(input_data=dataset)
h_symmetry.apply()
dx_theo = Quantity(dataset.theo_resolution)
dx_theo = dx_theo.to(u.arcsec)
print("theoretical formula for finest angular scale ", dx_theo)
print("recommended pixel size", dx_theo / 7.)
if dx == None:
print("using theoretical formula for pixel size")
dx = dx_theo / 7.
else:
print("sky image pixels: ", dx.to(u.arcsec))
# du = (1/(imsize*dx)).to(u.lambdas, equivalencies=lambdas_equivalencies())
print("imsize", imsize)
#gridder = Gridder(imsize=imsize,
# cellsize=dx,
# padding_factor=1.0,
# hermitian_symmetry=hermitian_symmetry)
dirty_mapper = DirtyMapper(
input_data=dataset,
imsize=imsize,
data_column=data_column,
padding_factor=1.,
cellsize=dx,
stokes="I", # "I,Q"
hermitian_symmetry=hermitian_symmetry)
dirty_images_natural = dirty_mapper.transform()
gridded_visibilities_nat = np.squeeze(
dirty_mapper.uvgridded_visibilities.compute())
gridded_weights_nat = np.squeeze(dirty_mapper.uvgridded_weights.compute())
if hermitian_symmetry:
gridded_visibilities_nat = completeuvplane_complex(
gridded_visibilities_nat)
gridded_weights_nat = completeuvplane(gridded_weights_nat)
if wantdirtymap:
#dirty_image_natural = dirty_images_natural[0].data[0].compute()
fits_io = FITS()
print("punching dirty map ", wantdirtymap)
fits_io.write(dirty_images_natural[0].data, output_name=wantdirtymap)
if wantpsf:
#dirty_beam_natural = dirty_images_natural[1].data[0].compute()
fits_io = FITS()
print("punching dirty beam ", wantpsf)
fits_io.write(dirty_images_natural[1].data, output_name=wantpsf)
return dx, gridded_visibilities_nat, gridded_weights_nat