diff options
| -rwxr-xr-x | astro/oskar/fits2skymodel.py | 67 |
1 files changed, 37 insertions, 30 deletions
diff --git a/astro/oskar/fits2skymodel.py b/astro/oskar/fits2skymodel.py index 8627324..90d5836 100755 --- a/astro/oskar/fits2skymodel.py +++ b/astro/oskar/fits2skymodel.py @@ -9,10 +9,11 @@ Convert a FITS image to OSKAR sky model for simulation usage. NOTE ---- -The OSKAR sky model consists of all the valid (>threshold) pixels -from the given image (slice), and fluxes are given in unit [Jy], -therefore, the input image should be converted from brightness -temperature [K] to unit [Jy/pixel]. +The OSKAR sky model consists of all the valid pixels (with absolute +values within the specified minimum and maximum thresholds) from +the given image (i.e., slice at a frequency channel), and fluxes +are given in unit [Jy], therefore, the input image should be converted +from brightness temperature [K] to unit [Jy/pixel]. References ---------- @@ -57,24 +58,28 @@ class SkyModel: ra0, dec0 : float The coordinate of the image center; unit [deg] minvalue : float, optional - The minimum threshold for the image values + The minimum threshold for the image absolute values + maxvalue : float, optional + The maximum threshold for the image absolute values mask : 2D bool `~numpy.ndarray`, optional Use this mask to select the sources of the output sky model, - instead of the above ``minvalue``. - NOTE: this overwrite the above ``minvalue`` if provided. + instead of the above ``minvalue`` and ``maxvalue``. + NOTE: Will overwrite the above ``minvalue`` and ``maxvalue``. projection : str, optional The WCS projection for the image; Default: "CAR" (Cartesian) TODO: support "SIN" etc. """ def __init__(self, image, freq, pixelsize, ra0, dec0, - minvalue=1e-4, mask=None, projection="CAR"): + minvalue=1e-4, maxvalue=np.inf, mask=None, + projection="CAR"): self.image = image # [K] (brightness temperature) self.freq = freq # [MHz] self.pixelsize = pixelsize # [arcsec] self.ra0 = ra0 # [deg] self.dec0 = dec0 # [deg] self.minvalue = minvalue + self.maxvalue = maxvalue self.mask = mask self.projection = projection logger.info("SkyModel: Loaded image @ %.2f [MHz], " % freq + @@ -125,7 +130,7 @@ class SkyModel: @property def factor_K2JyPixel(self): """ - Conversion factor to convert brightness unit from 'K' to 'Jy/pixel' + Conversion factor from [K] to [Jy/pixel] """ pixarea = (self.pixelsize * au.arcsec) ** 2 freq = self.freq * au.MHz @@ -161,8 +166,10 @@ class SkyModel: flux : source (Stokes I) flux density (Jy) """ if self.mask is None: - self.mask = np.abs(self.image) >= self.minvalue - logger.info("Use minimum threshold to determine output sky") + self.mask = ((np.abs(self.image) >= self.minvalue) & + (np.abs(self.image) <= self.maxvalue)) + logger.info("Use minimum and maximum thresholds: [%.4e, %.4e]" % + (self.minvalue, self.maxvalue)) else: logger.info("Use provided mask to determine output sky") idx = self.mask.flatten() @@ -188,8 +195,10 @@ class SkyModel: "RA0 = %.4f [deg]\n" % self.ra0 + "Dec0 = %.4f [deg]\n" % self.dec0 + "Minimum value = %.4e [K]\n" % self.minvalue + + "Maximum value = %.4e [K]\n" % self.maxvalue + "Source counts = %d (%.1f%%)\n\n" % (counts, percent) + "R.A.[deg] Dec.[deg] flux[Jy]") + logger.info("Writing sky model ...") np.savetxt(outfile, sky, fmt='%.10e, %.10e, %.10e', header=header) logger.info("Wrote OSKAR sky model to file: %s" % outfile) @@ -246,17 +255,17 @@ def main(): parser.add_argument("-f", "--freq", dest="freq", type=float, help="frequency [MHz] the image measured; " + "(default: obtain from the FITS header 'FREQ')") - exgrp = parser.add_mutually_exclusive_group() - exgrp.add_argument("-m", "--min-value", dest="minvalue", type=float, - help="minimum threshold to the output sky model; " + - "unit: [K]; (default: 1e-4, i.e., 0.1 mK)") - exgrp.add_argument("-M", "--min-peak-fraction", dest="minpfrac", - type=float, - help="minimum threshold determined as the fraction " + - "the peak value to the output sky model") - exgrp.add_argument("--mask", dest="mask", - help="use a mask to determine the output sky model") - # + parser.add_argument("-m", "--min-value", dest="minvalue", + type=float, default=1e-4, + help="[K] minimum threshold to the output sky model " + + "(default: 1e-4, i.e., 0.1 mK)") + parser.add_argument("-M", "--max-value", dest="maxvalue", + type=float, default=np.inf, + help="[K] maximum threshold to the output sky model " + + "(default: inf)") + parser.add_argument("--mask", dest="mask", + help="use a mask to determine the output sky model " + + "(NOTE: will override --min-value and --max-value)") parser.add_argument("-F", "--osm-fits", dest="osmfits", action="store_true", help="save a FITS version of the converted sky model") @@ -300,17 +309,15 @@ def main(): pixelsize = header["PixSize"] # [arcsec] logger.info("Frequency: %.2f [MHz]" % freq) logger.info("Pixel size: %.2f [arcsec]" % pixelsize) - minvalue = 1e-4 # i.e., 0.1 [mK] - mask = None - if args.minvalue: - minvalue = args.minvalue - if args.minpfrac: - minvalue = args.minpfrac * image.max() if args.mask: mask = fits.open(args.mask)[0].data.astype(np.bool) - logger.info("Minimum threshold: %g [K]" % minvalue) + logger.info("Loaded sky mask from file: %s" % args.mask) + else: + mask = None + logger.info("Threshold: %g - %g [K]" % (args.minvalue, args.maxvalue)) skymodel = SkyModel(image=image, freq=freq, ra0=args.ra0, dec0=args.dec0, - pixelsize=pixelsize, minvalue=minvalue, mask=mask) + pixelsize=pixelsize, minvalue=args.minvalue, + maxvalue=args.maxvalue, mask=mask) logger.info("Conversion [K] -> [Jy/pixel]: %g" % skymodel.factor_K2JyPixel) skymodel.write_sky_model(outfile, clobber=args.clobber) if args.osmfits: |