aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
-rwxr-xr-xbin/make-ska1low-model63
-rw-r--r--fg21sim/uvsim/__init__.py2
-rw-r--r--fg21sim/uvsim/telescope.py218
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)