From 70cf735e7bee9c320229432400c775f6da8f0a5d Mon Sep 17 00:00:00 2001 From: Aaron LI Date: Sat, 1 Jul 2017 01:54:22 +0800 Subject: Add astro/osm-randomize-coord.py Randomize the coordinates of the OSKAR sky model by adding an offset sampling from a Gaussian of given sigma. Try to reproduce the source subtraction position errors investigated by Datta et al. 2010 ApJ. --- astro/oskar/osm-randomize-coord.py | 163 +++++++++++++++++++++++++++++++++++++ 1 file changed, 163 insertions(+) create mode 100755 astro/oskar/osm-randomize-coord.py (limited to 'astro') diff --git a/astro/oskar/osm-randomize-coord.py b/astro/oskar/osm-randomize-coord.py new file mode 100755 index 0000000..452cd9d --- /dev/null +++ b/astro/oskar/osm-randomize-coord.py @@ -0,0 +1,163 @@ +#!/usr/bin/env python3 +# +# Copyright (c) 2017 Weitian LI +# MIT License +# + +""" +Create new randomized coordinates by adding random offset to the existing +OSKAR sky model (i.e., osm), and replace original coordinates with the +specified new ones. +""" + +import os +import argparse + +import numpy as np + + +class OskarSkyModel: + """ + OSKAR sky model + """ + def __init__(self, infile): + self.filename = infile + self.header = self.get_header(infile) + self.load_data(infile) + + @staticmethod + def get_header(infile): + """ + Read the OSM header lines + """ + header = [] + with open(infile) as f: + for line in f.readlines(): + if line[0] == "#": + header.append(line) + else: + break + print("Read OSM header:\n%s" % "".join(header)) + return header + + def load_data(self, infile): + try: + data = np.loadtxt(infile) + except ValueError: + data = np.loadtxt(infile, delimiter=",") + self.ra = data[:, 0] # [deg] + self.dec = data[:, 1] # [deg] + self.flux = data[:, 2] # [Jy] + self.number = len(self.flux) + print("Loaded OSM data from file: %s" % infile) + + def randomize_coord(self, sigma): + """ + Randomize the coordinates by adding an offset sampling + from a Gaussian of sigma. + """ + self.offset_sigma = sigma/3600.0 # [arcsec] -> [deg] + print("Random offset: %.1f [arcsec] == %.6f [deg]" % + (sigma, self.offset_sigma)) + self.ra += np.random.normal(loc=0.0, scale=self.offset_sigma, + size=self.number) + self.dec += np.random.normal(loc=0.0, scale=self.offset_sigma, + size=self.number) + print("Generated randomized coordinates") + + def replace_coord(self, coordfile): + """ + Replace the coordinates with the data from the given + coordinate file. + """ + try: + coord_new = np.loadtxt(coordfile) + except ValueError: + coord_new = np.loadtxt(coordfile, delimiter=",") + ra_new = coord_new[:, 0] + dec_new = coord_new[:, 1] + if self.number != len(ra_new): + raise RuntimeError("invalid coordinate file: %s" % coordfile) + self.ra = ra_new + self.dec = dec_new + print("Replaced coordinates") + + def save_data(self, data, outfile, clobber=False): + if os.path.exists(outfile) and not clobber: + raise OSError("file already exists: %s" % outfile) + with open(outfile, "wb") as fb: + for line in self.header: + fb.write(line.encode("utf-8")) + np.savetxt(fb, data) + + def save_coord(self, outfile, clobber=False): + data = np.column_stack([self.ra, self.dec]) + self.save_data(data=data, outfile=outfile, clobber=clobber) + print("Wrote coordinates to file: %s" % outfile) + + def save_osm(self, outfile, clobber=False): + data = np.column_stack([self.ra, self.dec, self.flux]) + self.save_data(data=data, outfile=outfile, clobber=clobber) + print("Wrote OSM to file: %s" % outfile) + + +def cmd_create(args): + """ + sub-command: create - create randomized coordinates + """ + osm = OskarSkyModel(args.infile) + osm.randomize_coord(sigma=args.sigma) + osm.save_coord(outfile=args.coordfile, clobber=args.clobber) + + +def cmd_replace(args): + """ + sub-command: create - create randomized coordinates + """ + osm = OskarSkyModel(args.infile) + osm.replace_coord(coordfile=args.coordfile) + if not args.nobackup: + backfile = args.infile + ".bak" + os.rename(args.infile, backfile) + print("Backed up OSM as: %s" % backfile) + osm.save_osm(outfile=args.infile, clobber=True) + + +def main(): + parser = argparse.ArgumentParser( + description="Randomize OSKAR sky model source coordinates") + subparsers = parser.add_subparsers(dest="subparser_name", + title="sub-commands", + help="additional help") + # sub-command: "create" + parser_create = subparsers.add_parser( + "create", help="create randomized coordinates") + parser_create.add_argument("-C", "--clobber", dest="clobber", + action="store_true", + help="overwrite existing output file") + parser_create.add_argument("-s", "--sigma", dest="sigma", + required=True, type=float, + help="random offset sigma [arcsec]") + parser_create.add_argument("-c", "--coord-file", dest="coordfile", + required=True, + help="output coordinate file") + parser_create.add_argument("infile", help="input OSKAR sky model") + parser_create.set_defaults(func=cmd_create) + # sub-command: "replace" + parser_replace = subparsers.add_parser( + "replace", help="replace coordinates of OSM") + parser_replace.add_argument("-B", "--no-backup", dest="nobackup", + action="store_true", + help="do NOT backup original OSM") + parser_replace.add_argument("-c", "--coord-file", dest="coordfile", + required=True, + help="file of new coordinates") + parser_replace.add_argument("infile", help="input OSKAR sky model") + parser_replace.set_defaults(func=cmd_replace) + # + args = parser.parse_args() + args.func(args) + + +if __name__ == "__main__": + main() -- cgit v1.2.2