astro/oskar: Add beamsize.py to summary beam sizes of multiple images

1 files changed, 59 insertions, 0 deletions

diff --git a/astro/oskar/beamsize.py b/astro/oskar/beamsize.py
new file mode 100755
index 0000000..8f0b018
--- /dev/null
+++ b/astro/oskar/beamsize.py
@@ -0,0 +1,59 @@
+#!/usr/bin/env python3
+#
+# Copyright (c) 2017 Weitian LI <weitian@aaronly.me>
+# MIT license
+#
+
+"""
+Read the beam parameters from one or more FITS images (e.g., imaged by
+WSClean), and briefly summarize the results.
+"""
+
+import argparse
+
+import numpy as np
+
+from astropy.io import fits
+
+
+def main():
+    parser = argparse.ArgumentParser(
+        description="Read and summarize the beam parameters")
+    parser.add_argument("infile", nargs="+", help="input FITS images")
+    args = parser.parse_args()
+
+    nfile = len(args.infile)
+    beam = np.zeros((nfile, 4))
+    for i, fn in enumerate(args.infile):
+        header = fits.open(fn)[0].header
+        bmajor = header["BMAJ"] * 3600  # [arcsec]
+        bminor = header["BMIN"] * 3600  # [arcsec]
+        bpa = header["BPA"]  # [deg]
+        bsize = np.pi * bmajor*bminor / (4*np.log(2))  # [arcsec^2]
+        beam[i, :] = (bsize, bmajor, bminor, bpa)
+
+    if nfile == 1:
+        print("size: %.1f [arcsec^2]" % beam[0, 0])
+        print("major: %.1f [arcsec]" % beam[0, 1])
+        print("minor: %.1f [arcsec]" % beam[0, 2])
+        print("position angle: %.1f [deg]" % beam[0, 3])
+        return 0
+
+    # Two or more input images ...
+    fnlen = max([len(fn) for fn in args.infile])
+    fmt_header = "%-{fnlen}s: %8s  %6s  %6s  %6s".format(fnlen=fnlen)
+    fmt_beam = "%-{fnlen}s: %8.0f  %6.1f  %6.1f  %6.1f".format(fnlen=fnlen)
+    print(fmt_header % ("filename", "size", "major", "minor", "pangle"))
+    print(fmt_header % (" ", "arcsec^2", "arcsec", "arcsec", "deg"))
+    print("------------------------------------------------------------"
+          "--------------------")
+    for i, fn in enumerate(args.infile):
+        print(fmt_beam % ((fn,) + tuple(beam[i, :])))
+    print("------------------------------------------------------------"
+          "--------------------")
+    print(fmt_beam % (("mean",) + tuple(beam.mean(axis=0))))
+    print(fmt_beam % (("std",) + tuple(beam.std(axis=0))))
+
+
+if __name__ == "__main__":
+    main()