diff options
| -rwxr-xr-x | bin/make-ska1low-model | 63 | ||||
| -rw-r--r-- | fg21sim/uvsim/__init__.py | 2 | ||||
| -rw-r--r-- | fg21sim/uvsim/telescope.py | 218 |
3 files changed, 283 insertions, 0 deletions
diff --git a/bin/make-ska1low-model b/bin/make-ska1low-model new file mode 100755 index 0000000..9ef0034 --- /dev/null +++ b/bin/make-ska1low-model @@ -0,0 +1,63 @@ +#!/usr/bin/env python3 +# +# Copyright (c) 2017 Weitian LI <weitian@aaronly.me> +# MIT license +# + +""" +Make SKA1-low telescope model for OSKAR[1] simulation usage. + +[1] OSKAR: https://github.com/OxfordSKA/OSKAR +""" + +import os +import sys +import argparse +import logging + +from pkg_resources import resource_exists, resource_filename + +from fg21sim.configs import configs +from fg21sim.utils import setup_logging +from fg21sim.uvsim import telescope + + +def main(): + layoutfile = "ska1low.layout.txt" + layoutfile_required = True + if resource_exists(telescope.__package__, layoutfile): + layoutfile = resource_filename(telescope.__package__, layoutfile) + layoutfile_required = False + else: + layoutfile = None + + parser = argparse.ArgumentParser( + description="Make SKA1-low telescope model for OSKAR") + parser.add_argument("-d", "--debug", dest="debug", action="store_true", + help="show verbose debug information") + parser.add_argument("-C", "--clobber", dest="clobber", + action="store_true", + help="overwrite the existing output files") + parser.add_argument("-l", "--layout-file", dest="layoutfile", + required=layoutfile_required, + default=layoutfile, + help="SKA1-low layout configuration " + + "(default: %s)" % layoutfile) + parser.add_argument("-o", "--outdir", dest="outdir", + required=True, + help="output telescope model directory") + args = parser.parse_args() + + loglevel = "DEBUG" if args.debug else "INFO" + setup_logging(dict_config=configs.logging, level=loglevel) + tool = os.path.basename(sys.argv[0]) + logger = logging.getLogger(tool) + logger.info("COMMAND: {0}".format(" ".join(sys.argv))) + + ska1low = telescope.SKA1Low(args.layoutfile) + ska1low.generate_stations() + ska1low.make_oskar_model(args.outdir, clobber=args.clobber) + + +if __name__ == "__main__": + main() diff --git a/fg21sim/uvsim/__init__.py b/fg21sim/uvsim/__init__.py new file mode 100644 index 0000000..e3f5b75 --- /dev/null +++ b/fg21sim/uvsim/__init__.py @@ -0,0 +1,2 @@ +# Copyright (c) 2017 Weitian LI <weitian@aaronly.me> +# MIT license diff --git a/fg21sim/uvsim/telescope.py b/fg21sim/uvsim/telescope.py new file mode 100644 index 0000000..5fe8281 --- /dev/null +++ b/fg21sim/uvsim/telescope.py @@ -0,0 +1,218 @@ +# Copyright (c) 2017 Weitian LI <weitian@aaronly.me> +# MIT license + +""" +Radio interferometer layout configurations. +""" + +import os +import logging +import shutil + +import numpy as np +import pandas as pd + +from .wgs84 import geodetic2enu + + +logger = logging.getLogger(__name__) + + +class SKA1Low: + """ + Process SKA1-low layout data and generate the telescope model + for OSKAR. + + Parameters + ---------- + infile : str + Path to the SKA1-low layout data file + stn_antennas : int, optional + Number of antenna elements per station (default: 256) + stn_diameter : float, optional + Diameter of each station (unit: [m]) + ant_min_sep : float, optional + Minimum separation between two antennas (unit: [m]) + + Reference + --------- + [1] SKA-TEL-SKO-0000422, revision 02, 2016-05-31, Table 1 + http://astronomers.skatelescope.org/wp-content/uploads/2016/09/SKA-TEL-SKO-0000422_02_SKA1_LowConfigurationCoordinates-1.pdf + [2] OSKAR: telescope model + http://www.oerc.ox.ac.uk/~ska/oskar2/OSKAR-Telescope-Model.pdf + """ + def __init__(self, infile, stn_antennas=256, stn_diameter=45.0, + ant_min_sep=1.5): + self.infile = infile + self.stn_antennas = stn_antennas + self.stn_diameter = stn_diameter + self.ant_min_sep = ant_min_sep + self.data = pd.read_csv(infile, sep="\s+", comment="#", + index_col="Label") + logger.info("Read telescope layout data from: %s" % infile) + self.position_wgs84 = np.array(self.data.loc["CENTER", :]) + logger.info("Telescope center coordinate: (%f, %f)" % + tuple(self.position_wgs84)) + self.labels = self.make_station_labels(self.data.index[1:]) + # (longitudes, latitudes) + self.layouts_wgs84 = np.array(self.data.iloc[1:, :]) + logger.info("Number of stations: %d" % len(self.layouts_wgs84)) + + def generate_stations(self): + """ + Generate the antenna elements layouts for each station. + """ + layouts = [] + N = len(self.labels) + logger.info("Number of antennas per station: %d" % + self.stn_antennas) + logger.info("Station diameter: %.2f [m]" % self.stn_diameter) + logger.info("Station antennas minimum separation: %.2f [m]" % + self.ant_min_sep) + logger.info("Generating antenna elements layouts ...") + for i, label in enumerate(self.labels): + logger.debug("Generate layout for [#%d/%d] station: %s" % + (i+1, N, label)) + x, y, __ = self.rand_uniform_2d( + n=self.stn_antennas, r_max=self.stn_diameter/2.0, + min_sep=self.ant_min_sep) + layouts.append((x, y)) + self.stn_layouts = layouts + logger.info("DONE generate station layouts.") + + def make_oskar_model(self, outdir, clobber=False): + """ + Create the telescope model for OSKAR. + """ + if os.path.exists(outdir): + if clobber: + shutil.rmtree(outdir) + logger.warning("Removed existing model: %s" % outdir) + else: + raise FileExistsError("Output directory already exists: " % + outdir) + os.mkdir(outdir) + logger.info("Created telescope model at: %s" % outdir) + # Write position + fposition = os.path.join(outdir, "position.txt") + open(fposition, "w").writelines([ + "# SKA1-low layout: %s\n" % self.infile, + "# Telescope center position (WGS84)\n", + "# longitude[deg] latitude[deg]\n", + "%.8f %.8f\n" % tuple(self.position_wgs84) + ]) + logger.info("Wrote telescope position: %s" % fposition) + # Write layout of stations + p0 = [self.position_wgs84[0], self.position_wgs84[1], 0.0] + layouts = np.array([geodetic2enu((lon, lat, 0.0), p0) + for lon, lat in self.layouts_wgs84]) + layouts[:, 2] = 0.0 # set `up` to 0.0 + flayout = os.path.join(outdir, "layout.txt") + header = ["SKA1-low layout: %s" % self.infile, + "All stations layout", + "East[m] North[m] Up[m]"] + np.savetxt(flayout, layouts, header="\n".join(header)) + logger.info("Wrote station layouts: %s" % flayout) + # Write stations + N = len(self.labels) + for i, label in enumerate(self.labels): + stn_dir = os.path.join(outdir, label) + os.mkdir(stn_dir) + fstation = os.path.join(stn_dir, "layout.txt") + header = [ + "Antenna elements layout", + "Station label: %s" % label, + "Number of antennas: %d" % self.stn_antennas, + "Station diameter: %.2f [m]" % self.stn_diameter, + "Antenna minimum separation: %.2f [m]" % self.ant_min_sep, + "X[m] Y[m]" + ] + np.savetxt(fstation, np.column_stack(self.stn_layouts[i]), + header="\n".join(header)) + logger.debug("Wrote layout for [#%d/%d] station: %s" % + (i+1, N, fstation)) + logger.info("DONE wrote telescope model: %s" % outdir) + + @staticmethod + def make_station_labels(labels, base="stn"): + """ + Make the labels for each station, which will also be used + as the sub-directory names for the output telescope model. + """ + N = len(labels) + ndigits = int(np.log10(N)) + 1 + fmt = "{base}.%(id)0{ndigits}d.%(label)s".format( + base=base, ndigits=ndigits) + stnlabels = [fmt % {"id": i+1, "label": l} + for i, l in enumerate(labels)] + return stnlabels + + @staticmethod + def rand_uniform_2d(n, r_max, min_sep, r_min=None): + """ + Generate 2D random points with a minimum separation within + a radius range (i.e., annulus/circle). + + Credit: + * https://github.com/OxfordSKA/SKA1-low-layouts : + layouts/utilities/layout.py - Layout.rand_uniform_2d() + """ + grid_size = min(100, int(np.ceil(r_max * 2.0) / min_sep)) + grid_cell = r_max * 2.0 / grid_size + scale = 1.0 / grid_cell + + x, y = np.zeros(n), np.zeros(n) + grid = { + "start": np.zeros((grid_size, grid_size), dtype=np.int), + "end": np.zeros((grid_size, grid_size), dtype=np.int), + "count": np.zeros((grid_size, grid_size), dtype=np.int), + "next": np.zeros(n, dtype=int) + } + + num_tries, max_tries, total_tries = 0, 0, 0 + for j in range(n): + done = False + while not done: + xt, yt = np.random.rand(2) * 2*r_max - r_max + rt = (xt**2 + yt**2) ** 0.5 + if rt + (min_sep/2.0) > r_max: + num_tries += 1 + elif r_min and rt - (min_sep/2.0) < r_min: + num_tries += 1 + else: + jx = int(round(xt + r_max) * scale) + jy = int(round(yt + r_max) * scale) + x0 = max(0, jx-2) + x1 = min(grid_size, jx+3) + y0 = max(0, jy-2) + y1 = min(grid_size, jy+3) + # Find the minimum spacing between the trial point + # and other points + d_min = r_max * 2.0 + for ky in range(y0, y1): + for kx in range(x0, x1): + if grid["count"][ky, kx] > 0: + i_other = grid["start"][ky, kx] + for kh in range(grid["count"][ky, kx]): + dx = xt - x[i_other] + dy = yt - y[i_other] + d_other = (dx**2 + dy**2) ** 0.5 + d_min = min(d_min, d_other) + i_other = grid["next"][i_other] + if d_min >= min_sep: + x[j], y[j] = xt, yt + if grid["count"][jy, jx] == 0: + grid["start"][jy, jx] = j + else: + grid["next"][grid["end"][jy, jx]] = j + grid["end"][jy, jx] = j + grid["count"][jy, jx] += 1 + max_tries = max(max_tries, num_tries) + total_tries += num_tries + num_tries = 0 + done = True + else: + num_tries += 1 + + info = {"max_tries": max_tries, "total_tries": total_tries} + return (x, y, info) |