oskar: Delete obsolete scripts runOSKAR.py & correct-pb.py

* runOSKAR.py: superseded by "fits2skymodel.py" and "run_oskar.py"
* correct-pb.py: it is generally not recommended to correct for
  primary beam

2 files changed, 0 insertions, 652 deletions

diff --git a/astro/oskar/correct-pb.py b/astro/oskar/correct-pb.py
deleted file mode 100755
index cac20ff..0000000
--- a/astro/oskar/correct-pb.py
+++ /dev/null
@@ -1,54 +0,0 @@
-#!/usr/bin/env python3
-#
-# Copyright (c) Weitna LI <weitian@aaronly.me>
-# MIT License
-#
-# Correct for primary beam by dividing the image by the (simulated)
-# primary beam pattern.
-#
-# 2017-06-20
-#
-
-import os
-import sys
-import argparse
-
-from astropy.io import fits
-
-
-def main():
-    parser = argparse.ArgumentParser(
-        description="Correct for the primary beam pattern")
-    parser.add_argument("-C", "--clobber", dest="clobber",
-                        action="store_true",
-                        help="overwrite existing files")
-    parser.add_argument("-p", "--primary-beam", dest="pb",
-                        required=True,
-                        help="file of the primary beam pattern")
-    parser.add_argument("infile", help="input image to be corrected for")
-    parser.add_argument("outfile", nargs="?",
-                        help="output pb-corrected image (default: add " +
-                        "'pbcorr.fits' suffix)")
-    args = parser.parse_args()
-
-    if args.outfile:
-        outfile = args.outfile
-    else:
-        outfile = os.path.splitext(args.infile)[0] + ".pbcorr.fits"
-
-    with fits.open(args.infile) as f:
-        imgin = f[0].data
-        header = f[0].header
-    pb = fits.open(args.pb)[0].data
-    imgout = imgin / pb
-    header.add_history(" ".join(sys.argv))
-    hdu = fits.PrimaryHDU(data=imgout, header=header)
-    try:
-        hdu.writeto(outfile, overwrite=args.clobber)
-    except TypeError:
-        hdu.writeto(outfile, clobber=args.clobber)
-    print("Wrote pb-corrected image: %s" % outfile)
-
-
-if __name__ == "__main__":
-    main()
diff --git a/astro/oskar/runOSKAR.py b/astro/oskar/runOSKAR.py
deleted file mode 100755
index 57d9b7f..0000000
--- a/astro/oskar/runOSKAR.py
+++ /dev/null
@@ -1,598 +0,0 @@
-#!/usr/bin/env python3
-#
-# Copyright (c) 2017 Weitian LI <liweitianux@live.com>
-# MIT license
-#
-# 2017-04-07
-#
-# Change logs:
-# 2017-06-08:
-#   * Add more command line arguments
-#   * Set "dryrun" to True
-#   * Do not output OSKAR binary visibility file (output MS only)
-# 2017-05-05:
-#   * Fix ConfigParser getfloat error
-#   * Workaround old astropy which uses "clobber" instead of "overwrite"
-#
-
-"""
-Run OSKAR to simulate the visibilities from the sky model specified
-by a FITS image.
-
-
-Credits
--------
-[1] GitHub: OxfordSKA/OSKAR
-    https://github.com/OxfordSKA/OSKAR
-[2] GitHub: OxfordSKA/EoR - Emma_files/sim_tidy.py
-    https://github.com/OxfordSKA/EoR/blob/master/Emma_files/sim_tidy.py
-"""
-
-import os
-import sys
-import subprocess
-import configparser
-import argparse
-import logging
-
-import numpy as np
-import astropy.io.fits as fits
-import astropy.constants as ac
-from astropy.wcs import WCS
-
-
-logging.basicConfig(level=logging.INFO)
-logger = logging.getLogger(os.path.basename(sys.argv[0]))
-
-
-class Settings:
-    """
-    OSKAR settings manager.
-    """
-    def __init__(self, infile):
-        self.infile = infile
-        self.config = configparser.ConfigParser(interpolation=None)
-        self.config.read(infile)
-        logger.info("Read in configuration file: %s" % infile)
-        self.init_oskar_settings()
-        self.update_oskar_settings(self.config)
-
-    @property
-    def my_settings(self):
-        return self.config["my"]
-
-    @property
-    def dryrun(self):
-        return self.my_settings.getboolean("dryrun", fallback=True)
-
-    @property
-    def clobber(self):
-        return self.my_settings.getboolean("clobber", fallback=False)
-
-    @property
-    def quiet(self):
-        return self.my_settings.getboolean("quiet", fallback=False)
-
-    @property
-    def oskar_bin(self):
-        oskar = self.my_settings.get("oskar_bin",
-                                     fallback="oskar_sim_interferometer")
-        return os.path.expanduser(oskar)
-
-    @property
-    def output_settings_fn(self):
-        """
-        String format pattern for the output OSKAR settings file.
-        """
-        default = "settings/sim_interferometer_{freq:.2f}.ini"
-        return self.my_settings.get("output_settings_fn", fallback=default)
-
-    @property
-    def output_skymodel_fn(self):
-        """
-        String format pattern for the output OSKAR sky model file.
-        """
-        default = "skymodel/skymodel_{freq:.2f}.txt"
-        return self.my_settings.get("output_skymodel_fn", fallback=default)
-
-    @property
-    def output_skyfits_fn(self):
-        """
-        String format pattern for the output FITS slice of the sky model.
-        """
-        default = "skymodel/skymodel_{freq:.2f}.fits"
-        return self.my_settings.get("output_skyfits_fn", fallback=default)
-
-    @property
-    def output_ms_fn(self):
-        """
-        String format pattern for the output simulated visibility
-        data in MeasurementSet format.
-        """
-        default = "visibility/visibility_{freq:.2f}.ms"
-        return self.my_settings.get("output_ms_fn", fallback=default)
-
-    @property
-    def output_vis_fn(self):
-        """
-        String format pattern for the output simulated visibility
-        data in OSKAR binary format.
-        """
-        default = "visibility/visibility_{freq:.2f}.oskar"
-        return self.my_settings.get("output_vis_fn", fallback=default)
-
-    @property
-    def telescope_model(self):
-        """
-        Telescope model used for visibility simulations.
-        """
-        return self.my_settings["telescope_model"]
-
-    @property
-    def input_cube(self):
-        """
-        Input FITS spectral cube.
-        """
-        return self.my_settings["input_cube"]
-
-    @property
-    def image_size(self):
-        """
-        Width/X and height/Y of the input FITS image (unit: pixel)
-        """
-        size = self.my_settings["image_size"].split(",")
-        return (int(size[0]), int(size[1]))
-
-    @property
-    def image_pixsize(self):
-        """
-        Pixel size of the input FITS image (unit: arcsec)
-        """
-        return self.my_settings.getfloat("image_pixsize")
-
-    @property
-    def frequency(self):
-        """
-        Frequency of the input image. (unit: MHz)
-
-        NOTE: required if the above input FITS file is not a cube, but
-              a 2D image.
-        """
-        return self.my_settings.getfloat("frequency")
-
-    @property
-    def bandwidth(self):
-        """
-        Bandwidth of the input image. (unit: MHz)
-        """
-        return self.my_settings.getfloat("bandwidth")
-
-    @property
-    def ra0(self):
-        """
-        R.A. of the center of the input sky field.
-        unit: deg
-        """
-        return self.my_settings.getfloat("ra0", fallback=0.0)
-
-    @property
-    def dec0(self):
-        """
-        Dec. of the center of the input sky field.
-        unit: deg
-        """
-        return self.my_settings.getfloat("dec0", fallback=-27.0)
-
-    @property
-    def use_gpus(self):
-        """
-        Whether to GPUs
-        """
-        return self.my_settings.getboolean("use_gpus", fallback=False)
-
-    @property
-    def start_time(self):
-        """
-        Start time of the simulating observation
-        """
-        # This default time keeps 'EoR0' region above horizon for 12 hours.
-        # SKA EoR0 region: (ra, dec) = (0, -27) [deg]
-        default = "2000-01-01T03:30:00.000"
-        return self.my_settings.get("start_time", fallback=default)
-
-    @property
-    def obs_length(self):
-        """
-        Observation length of time (unit: s).
-        """
-        default = 12.0 * 3600  # 12 hours
-        return self.my_settings.getfloat("obs_length", fallback=default)
-
-    @property
-    def obs_interval(self):
-        """
-        Observation interval providing the number of time steps in the
-        output data (unit: s).
-        """
-        default = 10.0  # [s]
-        return self.my_settings.getfloat("obs_interval", fallback=default)
-
-    @property
-    def time_average(self):
-        """
-        Correlator time-average duration to simulate time-averaging smearing
-        (unit: s).
-        """
-        default = 10.0  # [s]
-        return self.my_settings.getfloat("time_average", fallback=default)
-
-    def init_oskar_settings(self):
-        """
-        Initialize a `ConfigParser` instance with the default settings
-        for 'oskar_sim_interferometer'.
-        """
-        settings = configparser.ConfigParser()
-        settings.read_dict({
-            "General": {
-                "app": "oskar_sim_interferometer",
-            },
-            "simulator": {
-                "use_gpus": self.use_gpus,
-                "max_sources_per_chunk": 65536,
-                "double_precision": "true",
-                "keep_log_file": "true",
-            },
-            "sky": {
-                "advanced/apply_horizon_clip": "false",
-            },
-            "observation": {
-                "phase_centre_ra_deg": self.ra0,
-                "phase_centre_dec_deg": self.dec0,
-                "start_time_utc": self.start_time,
-                "length": self.obs_length,
-                "num_time_steps":
-                    int(np.ceil(self.obs_length/self.obs_interval)),
-                "num_channels": 1,
-            },
-            "telescope": {
-                "input_directory": self.telescope_model,
-                "pol_mode": "Scalar",
-                "normalise_beams_at_phase_centre": "true",
-                "allow_station_beam_duplication": "true",
-                "aperture_array/array_pattern/enable": "true",
-                "aperture_array/element_pattern/functional_type": "Dipole",
-                "aperture_array/element_pattern/dipole_length": 0.5,
-                "aperture_array/element_pattern/dipole_length_units":
-                    "Wavelengths",
-                "station_type": "Aperture array",
-            },
-            "interferometer": {
-                "channel_bandwidth_hz": self.bandwidth * 1e6,
-                "time_average_sec": self.time_average,
-                "uv_filter_min": "min",
-                "uv_filter_max": "max",
-                "uv_filter_units": "Wavelengths",
-            }
-        })
-        self.oskar_settings = settings
-        logger.info("Initialized 'oskar_settings'")
-
-    def update_oskar_settings(self, config):
-        """
-        Update the OSKAR settings with the loaded user configurations.
-        """
-        for section in self.oskar_settings.sections():
-            if section in config:
-                for key, value in config[section].items():
-                    self.oskar_settings[section][key] = value
-                    logger.info("oskar_settings: [%s]%s = %s" % (
-                        section, key, value))
-        logger.info("Updated 'oskar_settings'")
-
-    def write_oskar_settings(self, outfile, clobber=False):
-        """
-        Write the settings file for 'oskar_sim_interferometer'.
-        """
-        if os.path.exists(outfile) and (not clobber):
-            raise OSError("oskar settings file already exists: " % outfile)
-        with open(outfile, "w") as fp:
-            # NOTE: OSKAR do NOT like space around '='
-            self.oskar_settings.write(fp, space_around_delimiters=False)
-        logger.info("Wrote oskar settings file: %s" % outfile)
-
-
-class SpectralCube:
-    """
-    Manipulate the FITS spectral cube.
-
-    NOTE: The FITS data as `numpy.ndarray` has the opposite index
-          ordering, which likes the Fortran style, i.e., fastest
-          changing axis last: data[frequency, y, x]
-    """
-    def __init__(self, infile):
-        self.infile = infile
-        with fits.open(infile) as hdulist:
-            self.header = hdulist[0].header
-            self.cube = hdulist[0].data
-        self.wcs = WCS(self.header)
-        logger.info("Loaded FITS spectral cube: %s" % infile)
-        logger.info("Spectral cube: width=%d, height=%d" %
-                    (self.width, self.height))
-        if not self.is_cube:
-            logger.warning("NOT a spectral cube!")
-        else:
-            logger.info("Number of frequencies: %d" % self.nfreq)
-
-    @property
-    def naxis(self):
-        return self.header["NAXIS"]
-
-    @property
-    def is_cube(self):
-        return self.naxis == 3
-
-    @property
-    def width(self):
-        """
-        Width of the image, i.e., X axis.
-        """
-        return self.header["NAXIS1"]
-
-    @property
-    def height(self):
-        """
-        Height of the image, i.e., Y axis.
-        """
-        return self.header["NAXIS2"]
-
-    @property
-    def nfreq(self):
-        return self.header["NAXIS3"]
-
-    @property
-    def frequencies(self):
-        """
-        Frequencies of this cube. (unit: MHz)
-        If the input file is not a cube, then return 'None'.
-        """
-        if not self.is_cube:
-            logger.warning("Input FITS file is not a spectral cube: %s" %
-                           self.infile)
-            return None
-
-        nfreq = self.nfreq
-        pix = np.zeros(shape=(nfreq, self.naxis), dtype=np.int)
-        pix[:, -1] = np.arange(nfreq)
-        world = self.wcs.wcs_pix2world(pix, 0)
-        freqMHz = world[:, -1] / 1e6  # Hz -> MHz
-        return freqMHz
-
-    def get_slice(self, nth=0):
-        """
-        Extract the specified nth frequency slice from the cube.
-        """
-        if not self.is_cube:
-            logger.warning("Input FITS file is not a spectral cube: %s" %
-                           self.infile)
-            return self.cube
-        else:
-            return self.cube[nth, :, :]
-
-
-class SkyModel:
-    """
-    OSKAR sky model.
-    """
-    def __init__(self, image, freq, pixsize, ra0, dec0):
-        self.image = image  # K (brightness temperature)
-        self.freq = freq  # MHz
-        self.pixsize = pixsize  # arcsec
-        self.ra0 = ra0  # deg
-        self.dec0 = dec0  # deg
-        logger.info("SkyModel: Loaded image @ %.2f [MHz]" % freq)
-
-    @property
-    def wcs(self):
-        """
-        WCS for the given image assuming the 'SIN' projection.
-        """
-        shape = self.image.shape
-        delta = self.pixsize / 3600.0  # deg
-        wcs_ = WCS(naxis=2)
-        wcs_.wcs.ctype = ["RA---SIN", "DEC--SIN"]
-        wcs_.wcs.crval = np.array([self.ra0, self.dec0])
-        wcs_.wcs.crpix = np.array([shape[1], shape[0]]) / 2.0 + 1
-        wcs_.wcs.cdelt = np.array([delta, delta])
-        return wcs_
-
-    @property
-    def fits_header(self):
-        header = self.wcs.to_header()
-        header["BUNIT"] = ("Jy/pixel", "Brightness unit")
-        header["FREQ"] = (self.freq, "Frequency [MHz]")
-        header["RA0"] = (self.ra0, "Center R.A. [deg]")
-        header["DEC0"] = (self.dec0, "Center Dec. [deg]")
-        return header
-
-    @property
-    def factor_K2JyPixel(self):
-        """
-        Conversion factor to convert brightness unit from 'K' to 'Jy/pixel'
-
-        http://www.iram.fr/IRAMFR/IS/IS2002/html_1/node187.html
-        """
-        pixarea = np.deg2rad(self.pixsize/3600.0) ** 2  # [sr]
-        kB = ac.k_B.si.value  # Boltzmann constant [J/K]
-        c0 = ac.c.si.value  # speed of light in vacuum [m/s]
-        freqHz = self.freq * 1e6  # [Hz]
-        factor = 2*kB * 1.0e26 * pixarea * (freqHz/c0)**2
-        return factor
-
-    @property
-    def ra_dec(self):
-        """
-        Calculate the (ra, dec) of each image pixel using the above WCS.
-
-        NOTE: axis ordering difference between numpy array and FITS
-        """
-        shape = self.image.shape
-        wcs = self.wcs
-        x, y = np.meshgrid(np.arange(shape[1]), np.arange(shape[0]))
-        pix = np.column_stack([x.flatten(), y.flatten()])
-        world = wcs.wcs_pix2world(pix, 0)
-        ra = world[:, 0].reshape(shape)
-        dec = world[:, 1].reshape(shape)
-        return (ra, dec)
-
-    @property
-    def sky(self):
-        """
-        OSKAR sky model array converted from the input image.
-
-        Columns
-        -------
-        ra : (J2000) right ascension (deg)
-        dec : (J2000) declination (deg)
-        flux : source (Stokes I) flux density (Jy)
-        """
-        ra, dec = self.ra_dec
-        ra = ra.flatten()
-        dec = dec.flatten()
-        flux = self.image.flatten() * self.factor_K2JyPixel
-        mask = flux > 1e-40
-        sky_ = np.column_stack([ra[mask], dec[mask], flux[mask]])
-        return sky_
-
-    def write_sky_model(self, outfile, clobber=False):
-        """
-        Write the converted sky model for simulation.
-        """
-        if os.path.exists(outfile) and (not clobber):
-            raise OSError("oskar sky model file already exists: " % outfile)
-        sky = self.sky
-        header = ("Frequency = %.3f [MHz]\n" % self.freq +
-                  "Pixel size = %.2f arcsec\n" % self.pixsize +
-                  "RA0 = %.4f [deg]\n" % self.ra0 +
-                  "Dec0 = %.4f [deg]\n" % self.dec0 +
-                  "Number of sources = %d\n\n" % len(sky) +
-                  "R.A.[deg]    Dec.[deg]    flux[Jy]")
-        np.savetxt(outfile, sky, fmt='%.10e, %.10e, %.10e', header=header)
-        logger.info("Wrote oskar sky model file: %s" % outfile)
-
-    def write_fits(self, outfile, oldheader=None, clobber=False):
-        if os.path.exists(outfile) and (not clobber):
-            raise OSError("Sky FITS already exists: " % outfile)
-        if oldheader is not None:
-            header = oldheader
-            header.extend(self.fits_header, update=True)
-        else:
-            header = self.fits_header
-        image = self.image * self.factor_K2JyPixel
-        hdu = fits.PrimaryHDU(data=image, header=header)
-        try:
-            hdu.writeto(outfile, overwrite=True)
-        except TypeError:
-            hdu.writeto(outfile, clobber=True)  # old astropy versions
-        logger.info("Wrote sky FITS to file: %s" % outfile)
-
-
-class Oskar:
-    """
-    Run OSKAR simulations
-    """
-    def __init__(self, settings):
-        self.settings = settings
-
-    def run(self, settingsfile, dryrun=True, shfile=None):
-        cmd = [self.settings.oskar_bin]
-        if self.settings.quiet:
-            cmd += ["--quiet"]
-        cmd += [settingsfile]
-        shellcmd = cmd[0]
-        for arg in cmd[1:]:
-            shellcmd += ' "%s"' % arg
-        logger.info("Running OSKAR simulator:")
-        logger.info("$ %s" % shellcmd)
-        if shfile:
-            open(shfile, "a").write("%s\n" % shellcmd)
-        if dryrun:
-            logger.info("Dry run!")
-        else:
-            subprocess.check_call(cmd)
-
-
-def main():
-    parser = argparse.ArgumentParser(
-        description="Run OSKAR to simulate visibilities")
-    parser.add_argument("-C", "--clobber", dest="clobber",
-                        action="store_true",
-                        help="overwrite existing files")
-    parser.add_argument("-R", "--run-oskar", dest="run",
-                        action="store_true",
-                        help="run OSKAR instead of just print commands")
-    parser.add_argument("-S", "--sh-file", dest="shfile",
-                        help="also write OSKAR commands to this file")
-    parser.add_argument("config", help="Configuration file")
-    args = parser.parse_args()
-
-    settings = Settings(args.config)
-    clobber = args.clobber if args.clobber else settings.clobber
-    dryrun = settings.dryrun
-    if args.run:
-        dryrun = False
-    if args.shfile and os.path.exists(args.shfile):
-        if clobber:
-            os.remove(args.shfile)
-        else:
-            raise OSError("Output file already exists: %s" % args.shfile)
-
-    image_cube = SpectralCube(settings.input_cube)
-    frequencies = image_cube.frequencies  # [MHz]
-    if frequencies is None:
-        frequencies = [settings.frequency]
-    logger.info("Number of image slices/frequencies: %d" % len(frequencies))
-
-    for nth, freq in enumerate(frequencies):
-        logger.info(">>> Processing #%d/%d image slice @ %.2f [MHz] <<<" %
-                    (nth+1, len(frequencies), freq))
-        settingsfile = settings.output_settings_fn.format(freq=freq)
-        skymodelfile = settings.output_skymodel_fn.format(freq=freq)
-        skyfitsfile = settings.output_skyfits_fn.format(freq=freq)
-        msfile = settings.output_ms_fn.format(freq=freq)
-        visfile = settings.output_vis_fn.format(freq=freq)
-        for filepath in [settingsfile, skymodelfile, skyfitsfile,
-                         msfile, visfile]:
-            dname = os.path.dirname(filepath)
-            if not os.path.isdir(dname):
-                os.makedirs(dname)
-
-        newconfig = configparser.ConfigParser()
-        newconfig.read_dict({
-            "sky": {
-                "oskar_sky_model/file": skymodelfile,
-            },
-            "observation": {
-                "start_frequency_hz": freq * 1e6,
-            },
-            "interferometer": {
-                "oskar_vis_filename": "",
-                "ms_filename": msfile,
-            },
-        })
-        settings.update_oskar_settings(newconfig)
-        settings.write_oskar_settings(outfile=settingsfile, clobber=clobber)
-
-        image_slice = image_cube.get_slice(nth)
-        skymodel = SkyModel(image=image_slice, freq=freq,
-                            pixsize=settings.image_pixsize,
-                            ra0=settings.ra0, dec0=settings.dec0)
-        skymodel.write_sky_model(skymodelfile, clobber=clobber)
-        skymodel.write_fits(skyfitsfile, oldheader=image_cube.header,
-                            clobber=clobber)
-
-        oskar = Oskar(settings)
-        oskar.run(settingsfile, dryrun=dryrun, shfile=args.shfile)
-
-
-if __name__ == '__main__':
-    main()