aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorAaron LI <aly@aaronly.me>2017-09-02 10:45:57 +0800
committerAaron LI <aly@aaronly.me>2017-09-02 10:46:32 +0800
commit04aef0c64a49f94f0e5db9c2ccc4bed63f7b4f5d (patch)
treed9faa3d191469fbd97a61c9e0220f17e060e372d
parentd31c8bd4901ca4400727f8ca6e68c94edc747d0c (diff)
downloadatoolbox-04aef0c64a49f94f0e5db9c2ccc4bed63f7b4f5d.tar.bz2
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
-rwxr-xr-xastro/oskar/correct-pb.py54
-rwxr-xr-xastro/oskar/runOSKAR.py598
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()