aboutsummaryrefslogtreecommitdiffstats
path: root/astro/oskar/wsclean.py
blob: 18e4b4c5113f355939b1731a6df35c97d2f8b40f (plain)
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
116
117
118
119
120
121
122
123
124
125
126
127
#!/usr/bin/env python3
#
# Copyright (c) 2017 Aaron LI
# MIT license
#
# Create image from OSKAR simulated visibility data using `WSClean`.
# WSClean: https://sourceforge.net/p/wsclean/
#
# 2017-09-01
#


import os
import argparse
import subprocess
import time


def wsclean(args):
    t1 = time.perf_counter()
    cmd = ["wsclean"] + args
    print("CMD: %s" % " ".join(cmd))
    subprocess.check_call(cmd)
    t2 = time.perf_counter()
    print("-----------------------------------------------------------")
    print("WSClean Elapsed time: %.1f [min]" % ((t2-t1)/60))
    print("-----------------------------------------------------------")


def main():
    parser = argparse.ArgumentParser(description="Run WSClean")
    parser.add_argument("--args", dest="args", nargs="*",
                        help="additional arguments for WSClean")
    parser.add_argument("--dirty", dest="dirty", action="store_true",
                        help="only create dirty images (by setting niter=0)")
    parser.add_argument("--update-model", dest="update_model",
                        action="store_true",
                        help="update the MODEL_DATA column in MS")
    parser.add_argument("--weight-briggs", dest="briggs",
                        type=float, default=0.0,
                        help="Briggs weight parameter")
    parser.add_argument("--niter", dest="niter", type=int, default=100000,
                        help="maximum number of CLEAN iterations")
    parser.add_argument("--gain", dest="gain", type=float, default=0.1,
                        help="CLEAN gain for each minor iteration")
    parser.add_argument("--mgain", dest="mgain", type=float, default=0.85,
                        help="CLEAN gain for major iterations")
    parser.add_argument("--size", dest="size", type=int,
                        required=True,
                        help="output image size (pixel number on a side)")
    parser.add_argument("--pixelsize", dest="pixelsize", type=float,
                        required=True,
                        help="output image pixel size [arcsec]")
    parser.add_argument("--taper-gaus", dest="taper_gaus", type=float,
                        help="taper the weights with a Gaussian function " +
                        "to reduce the contribution of long baselines. " +
                        "Gaussian beam size in [arcsec].")
    parser.add_argument("--fit-spec-order", dest="fit_spec_order", type=int,
                        help="do joined-channel CLEAN by fitting the " +
                        "spectra with [order] polynomial in normal-space")
    #
    exgrp = parser.add_mutually_exclusive_group()
    exgrp.add_argument("--threshold-auto", dest="threshold_auto",
                       type=float, default=1.5,
                       help="estimate noise level and stop at <sigma>*<std>")
    exgrp.add_argument("--threshold", dest="threshold", type=float,
                       help="stopping CLEAN threshold [Jy]")
    #
    parser.add_argument("--name", dest="name", required=True,
                        help="filename prefix for the output files")
    parser.add_argument("--ms", nargs="+", help="input visibility MSs")
    args = parser.parse_args()

    nms = len(args.ms)  # i.e., number of MS == number of channels

    cmdargs = [
        "-verbose",
        "-log-time",
        "-pol", "XX",  # OSKAR "Scalar" simulation only give "XX" component
        "-make-psf",  # always make the PSF, even no cleaning performed
    ]

    if args.dirty:
        cmdargs += ["-niter", str(0)]  # make dirty image only
    else:
        cmdargs += ["-niter", str(args.niter)]

    cmdargs += ["-weight", "briggs", str(args.briggs)]
    cmdargs += ["-gain", str(args.gain)]
    cmdargs += ["-mgain", str(args.mgain)]
    cmdargs += ["-size", str(args.size), str(args.size)]
    cmdargs += ["-scale", "{0}asec".format(args.pixelsize)]

    if args.fit_spec_order:
        cmdargs += ["-joinchannels", "-channelsout", str(nms),
                    "-fit-spectral-pol", str(args.fit_spec_order+1)]

    if args.update_model:
        cmdargs += ["-update-model-required"]
    else:
        cmdargs += ["-no-update-model-required"]

    if args.threshold:
        cmdargs += ["-threshold", str(args.threshold)]
    else:
        cmdargs += ["-auto-threshold", str(args.threshold_auto)]

    if args.taper_gaus:
        cmdargs += ["-taper-gaussian", str(args.taper_gaus)]

    # additional WSClean arguments
    if args.args:
        cmdargs += args.args

    nameprefix = args.name.rstrip("-_")
    cmdargs += ["-name", nameprefix]
    cmdargs += args.ms

    wsclean(cmdargs)

    if args.dirty:
        # Remove the output "-image" since it is identical to "-dirty"
        os.remove(nameprefix+"-image.fits")


if __name__ == "__main__":
    main()