Add uvsim/telescope.py and bin/make-ska1low-model

The `bin/make-ska1low-model` executable uses `telescope.py` to generate
the SKA1-low telescope model for OSKAR simulation usage.

1 files changed, 218 insertions, 0 deletions

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)