PyBDSF Documentation

PyBDSF (the Python Blob Detector and Source Finder) is a tool designed to decompose radio interferometry images into sources and make available their properties for further use. PyBDSF can decompose an image into a set of Gaussians, shapelets, or wavelets as well as calculate spectral indices and polarization properties of sources and measure the psf variation across an image. PyBDSF uses an interactive environment based on CASA [1] that will be familiar to most radio astronomers. Additionally, PyBDSF may also be used in Python scripts.

Why PyBDSF?

PyBDSF was developed to serve the needs of LOFAR, the LOw Frequency ARray [1], functioning primarily in The Netherlands. LOFAR achieves orders of magnitude better sensitivity and resolution at low frequencies (15-80 MHz, 120-240 MHz) than previous radio interferometers. In addition, given the large primary beam, the field of view is large (2-16 degrees), and the spectral coverage is wider than usual (up to 48 MHz bandwidth and up to 62,464 channels). Lastly, up to 244 independent beams can be electronically generated on the sky. These capabilities make LOFAR an ideal survey instrument, but also create challenges in data processing. In particular, for surveys a good source extraction software is essential. Before PyBDSF was developed, a survey was made of the existing source extraction packages (SAD in AIPS, SFIND in MIRIAD and SExtractor). It was concluded that none of these packages were adequate to the task, and further, it would be difficult to modify any of these to suit the needs of LOFAR. Hence, PyBDSF was written (first called PyBDSM). However, in the recent years, new low-frequency telescope projects have started which are a similar to LOFAR in some parameters (e.g., MEERKAT, Mileura Array, ASKAP and other SKA pathfinders) and there has been considerable effort to develop source extraction software at some these project sites, e.g. DUCHAMP [2].

Traditionally, source extraction software, at least in radio astronomy, has defined the process as fitting (multiple) Gaussians to source pixels. This makes sense since all interferometric images are convolved with a Gaussian (fit to the main lobe of the dirty beam) after deconvolution. This process is adequate also because most radio images have primarily consisted of point (or slightly extended) sources. LOFAR images, however, are quite different. Note that the antenna diameter is 50 m, maximum baselines extend to 100 km or more, and in addition, LOFAR has almost no missing short spacing measurements (unless flagged due to RFI). Hence, the images will have a much wider range of scales of emission than usual - from point sources up to extended 3C sources. Decomposing such sources into Gaussians may not be very effective (as well as highly degenerate and hence not very useful). Hence alternative basis sets which can capture a variety of scales are essential. AIPS (Classic AIPS as well as CASA) has been experimenting with multi-resolution CLEAN methods for many years now and in the same spirit, we have included shapelet and wavelet decomposition as well. With the kind of image morphologies LOFAR routinely produces, complex ways of describing sources are needed, not just to catalog them but also to perform other filtering operations post-extraction for science purposes. Note that although PyBDSF is written for LOFAR, it will obviously work for images from any radio interferometric telescope.

Footnotes

Capabilities of PyBDSF

PyBDSF can be run on FITS images (using Astropy [1]) or CASA images (using python-casacore [2]), including 3-D and 4-D cubes, and can handle blanked image pixels. If a spectral cube is given, then all source extraction as well as other computation (psf variation, wavelet decomposition, etc.) are done on a collapsed 2-D stokes I image. Once sources have been identified, their spectral and polarisation properties are then extracted from the full cubes. If you need a full 3-D Gaussian decomposition, then DUCHAMP [3] is what you need.

PyBDSF performs the following tasks:

Reads in the image, collapses specific frequency channels, with weights, and produces a ‘continuum’ image (the ‘ch0’ image) for all polarisations. The Stokes I ch0 image is used for all further computation.

Preprocessing is done, whereby some basic parameters like image statistics are computed. Also, any input parameters that are left to default are calculated using sensible algorithms.

The background rms and mean images are computed. If the variation in these images is not statistically significant, then a constant value is taken. The parameters for this calculation are computed generically and hence do not have information about, for example, the typical size of the artifacts around bright sources. These parameters (e.g., rms_box) are probably the only ones the user needs to take care to specify.

A constant threshold for separating source and noise pixels is set. This threshold can be either a hard threshold or calculated using the False Detection Rate algorithm.

Using these parameters, islands of contiguous source emission are identified. Islands are the basic units which are operated upon subsequently.

Each island is now fit with multiple Gaussians. Depending on the number of degrees of freedom, etc, the size of the fitted Gaussian could be fixed to be the restoring beam. The fitted Gaussians are then flagged to produce a list of acceptable set of Gaussians.

Each island can also be decomposed into shapelets. Currently only cartesian shapelets are implemented, and only one shapelet set, with the same scale in both dimensions, can be fit to an island. The shapelets parameters can be written out as ASCII or FITS tables.

Residual FITS images are computed, for both Gaussians and shapelets. The Gaussian parameters can be written out in various formats (ASCII, FITS tables, LOFAR BBS, ds9 region files, AIPS star, Kvis, etc). Shapelet parameters can be written out to FITS tables.

Gaussians within a given island are grouped into discrete sources.

If a frequency cube is input, then for each source identified in an island, the spectral index is computed. If possible, a spectral index is calculated for each Gaussian as well. This is done for point as well as extended sources.

If all four Stokes images are present, then the polarisation percentage and angle are calculated for each source.

The residual ch0 image, after subtracting fitted Gaussians, is processed using the à trous wavelet transform to generate images at various scales. Islands are identified in each of these wavelet images and fitted with Gaussians, all of which are then grouped to form pyramidal sources. These can be used further by the user as a starting point for morphological filters.

Since the ionosphere affects low frequencies significantly, PyBDSF can also estimate the spatial variation of the PSF across the image, which can be used to correct various source parameters.

Footnotes

Downloading and installing

The latest version of the code and installation instructions are available on the PyBDSF GitHub page at https://github.com/lofar-astron/PyBDSF. Please go there to obtain PyBDSF.

Note

Most of the documentation assumes that the pybdsf interactive shell is used, so it is recommended that you install it (it is not installed by default; see the installation instructions at the GitHub page above for details). However, if you are only interested in using PyBDSF inside Python scripts (see Scripting PyBDSF) and are already familiar with its usage, the interactive shell is not needed.

Note

If you are working on the LOFAR CEP3 cluster, then PyBDSF is already installed. All that needs to be done is to initialize your environment as follows:

$ module load lofim

What’s New

Version 1.10.3 (2023/05/22):

Fix build issue with Python 3.11 (#205)

Use cibuildwheel to build binary wheels for Linux and MacOS (Intel); drop support for Python 3.6 (#203)

Fix #198. Use the new method call canvas.manager.set_window_title (#199)

Replace Travis CI with GitHub actions (#196)

Version 1.10.2 (2023/02/10):

Fix issues with numpy versions >= 1.24 (#193)
Switch to manylinux2014 for building binary wheels (#191)
Fix ImportError in setuptools (#190)
Add binary wheels for Python 3.10 (#186)
Fix various documentation issues (#185)
Add logfilename option (#181)
Use len() instead of numpy.alen() (#180)

Version 1.10.1 (2022/02/14):

Fix NumPy API compatibility issue

Version 1.10.0 (2022/02/09):

Update some functions as required by scipy version >= 1.8.0 (PR #172)
Fix build issues with Python 3.8, end support for Python < 3.6, add support for Python 3.8 and 3.9, and make installation of the interactive pybdsf shell optional (PR #169)
Improve handling of the beam in the spectral index module (PR #165)
Improve handling of large, complex islands (PR #160)
Allow a file to be supplied for the ch0 image used in the spectral index module (PR #127)

Version 1.9.2 (2019/12/05):

Fix exception behaviour if spline order change does not work
Add check for frequency info in header

Version 1.9.1 (2019/09/25):

Fix various shapelet decomposition issues

Fix crash in Gaussian fitting (#96)

Fix blank_limit check_low error (#100)

Fix various minor bugs

Version 1.9.0 (2019/03/25):

Add support for Python 3

Fix various minor bugs

Version 1.8.15 (2018/10/12):

Fix segfault in Gaussian fitting (#63)

Fix math domain error (#76)

Fix setup.py for boost versions > 1.63 (#58)

Version 1.8.14 (2018/05/18):

Fix an error on total flux density (#50)

Add the possibility to provide an external noise and mean maps (#43)

Append the image FITS header into catalog FITS header (#53)

Make PyBDSF compatible with newer boost libraries, specifically those used in Ubuntu 18.04 (#55)

Version 1.8.13 (2017/11/17):

Remove deprecated boolean operators

Version 1.8.12 (2017/09/01):

Fix crash with tiny regions

Fix very low centroid peak fluxes

Fix compile error with numpy 1.13

Version 1.8.11 (2017/06/01):

Fix for interactive shell problem

Version 1.8.10 (2017/05/31):

Fixes for various installation and runtime issues on modern systems.

Version 1.8.9 (2017/03/23):

Fix to bug that causes an error when grouping Gaussians into sources

Version 1.8.8 (2017/03/17):

Rename to PyBDSF

Move to python setuptools

Move out of the LOFAR tree to github

Fix to issues related to numpy versions >= 1.12.0

Version 1.8.7 (2016/06/10):

Fix to bug that caused incorrect output images when input image was not square.

Version 1.8.6 (2016/01/20):

Fix to bug that caused incorrect island mask when two islands are very close together.

Fix to bug that caused crash when image is masked and the src_ra_dec option is used.

Version 1.8.5 (2015/11/30):

Fix to bug in export_image that resulted in incorrect output image when both trim_box and pad_image were used.

Fix to bug in wavelet module related to merging of islands.

Fix to bug in polarization module related to numbering of new islands.

Added option to use much faster (but also much more memory intensive) SciPy fftconvolve function instead of custom PyBDSM one. The option (use_scipy_fft) defaults to True.

Increased number of digits for values in output text catalogs

Version 1.8.4 (2015/08/06):

Improved speed of wavelet module.

Added option to use PyFFTW in wavelet module if available.

Fix to IPython version check.

Fix to bug that caused a failure to write shapelet models in FITS format.

Fix to bug that caused a crash when both atrous_do = True and output_all = True.

Fixed a bug that caused a crash on machines with only one core.

Version 1.8.3 (2014/09/26):

Fix to bug that caused a crash when using the wavelet module and all Gaussians in an island were flagged.

Version 1.8.2 (2014/05/14):

Island masks generated by the export_image task will now be expanded to match input image shape (as required by the AWimager).

Fix to bug that caused image read failure when image lacks a Stokes axis.

Fix to bug in CASA masks generated with export_image that caused cleaning to fail in CASA 4.2 and above.

Fix to bug that resulted in output file names being converted to lower case inappropriately.

Version 1.8.1 (2014/01/14):

Added option (bbs_patches = 'mask') to allow patches in an output BBS sky model to be defined using a mask image (set with the bbs_patches_mask option).

Fix to bug that caused the incl_empty option to be ignored when format = 'fits' in the write_catalog task.

Enabled output of images in CASA format in the export_image task (img_format = 'casa').

Added an option to export_image to export an island-mask image, with ones where there is emission and zeros elsewhere (image_type = 'island_mask'). Features in the island mask may be optionally dilated by specifying the number of dilation iterations with the mask_dilation parameter. The mask image may be padded with zeros to match the original image when the trim_box option was used to analyze only a portion of the image (pad_image = True).

Added an option to write a CASA region file to the write_catalog task (format = 'casabox').

Added an option to write a CSV catalog to the write_catalog task (format = 'csv').

Added error message when the rms is zero in some part of the rms map.

Version 1.8.0 (2013/10/16):

Improved wavelet fitting. Added option so that wavelet fitting can be done to the sum of images on the remaining wavelet scales, improving the signal for fitting (controlled with the atrous_sum option). Added option so that user can choose whether to include new islands found only in the wavelet images in the final fit or not (controlled with the atrous_orig_isl option).

Fixed a bug that could lead to incomplete fitting of some islands.

Improved overall convergence of fits.

Version 1.7.7 (2013/10/10):

Improved fitting of bright sources under certain circumstances.

Version 1.7.6 (2013/09/27):

Changed caching behavior to ensure that temporary files are always deleted after they are no longer needed or on exit.

Renamed blank_zeros to blank_limit. The blank_limit option now specifies a limit below which pixels are blanked.

Enabled SAGECAL sky-model output.

Version 1.7.5 (2013/09/02):

Fix to bug that caused a crash when images with 2 or 3 axes were used.

Improved rms and mean calculation (following the implementation used in PySE, see http://dare.uva.nl/document/174052 for details). The threshold used to determine the clipped rms and mean values is now determined internally by default (i.e., kappa_clip = None).

Version 1.7.4 (2013/08/29):

Fix to bug in show_fit that caused error when i is pressed in the plot window and shapelet decomposition had not been done.

Tweak to pybdsm startup shell script to avoid problems with the Mac OS X matplotlib backend on non-framework Python installations (such as Anaconda Python).

Fix to bug in process_image that could result in wavelet Gaussians being excluded from model image under certain conditions.

Version 1.7.3 (2013/08/27):

Fix to bug in image reading that caused images to be distorted.

Version 1.7.2 (2013/08/23):

Improved handling of non-standard FITS CUNIT keywords.

Improved loading of FITS images when trim_box is specified.

Version 1.7.1 (2013/08/22):

Fix to bug that caused cached images to be deleted when rerunning an analysis.

Fix to bug in show_fit due to undefined images.

Fix to bug in process_image (and img.process()) that would result in unneeded reprocessing.

Version 1.7.0 (2013/08/20):

PyBDSM will now use Astropy if installed for FITS and WCS modules.

Fix to avoid excessive memory usage in the wavelet module (replaced scipy.signal.fftconvolve with a custom function).

Added option to use disk caching for internally derived images (do_cache). Caching can reduce memory usage and is therefore useful when processing large images.

Implemented a variable operation chain for process_image (and img.process()) that allows unneeded steps to be skipped if the image is being reprocessed.

Fixed a bug that could cause Gaussian fitting to hang during iterative fitting of large islands.

Added option (fix_to_beam) to fix the size and position angle of Gaussians to the restoring beam during fitting.

Fix to bug that caused the position angle used to initialize fitting to be incorrect.

Version 1.6.1 (2013/03/22):

Fix to bug in ds9 and kvis catalog files that resulted in incorrect position angles.

Fix to bug in position-dependent WCS transformations that caused incorrect source parameters in output catalogs.

Added option to output uncorrected source parameters to a BBS sky model file (correct_proj).

Removed sky transformations for flagged Gaussians, as these could sometimes give math domain errors.

Disabled aperture flux measurement on wavelet images as it is not used/needed.

Version 1.6.0 (2013/03/05):

Improved speed and accuracy of aperture flux calculation.

Added option to use the curvature map method of Hancock et al. (2012) for the initial estimation of Gaussian parameters (ini_method = 'curvature') and for grouping of Gaussians into sources (group_method = 'curvature').

Fix to bug in spectral index module that caused sources with multiple Gaussians to be skipped. Minor adjustments to the wavelet module to improve performance.

Implemented position-dependent WCS transformations.

Added option to fit to any arbitrary location in the image within a given radius (src_ra_dec and src_radius_pix).

Fix to bug in wavelet module that caused crash when no Gaussians were fit to the main image.

Fix to bug that resulted in incorrect numbering of wavelet Gaussians. Added 'srl' output in ds9 format when using output_all = True.

Fix to bug in source grouping algorithm. Improved fitting when background mean is nonzero. Fix to allow images with GLAT and GLON WCS coordinates. Fix to bug when equinox is taken from the epoch keyword.

Version 1.5.1 (2012/12/19):

Fix to bug in wavelet module that occurred when the center of the wavelet Gaussian lies outside of the image.

Fix to re-enable srl output catalogs in ds9 region format.

Fix to bug that resulted in the output directory not always being created.

Added an option (aperture_posn), used when aperture fluxes are desired, to specify whether to center the aperture on the source centroid or the source peak.

Changes to reduce memory usage, particularly in the wavelet module.

Fix to bypass bug in matplotlib when display variable is not set.

Fixed bug that caused a crash when a detection image was used.

Fixed a bug with incorrect save directory when “plot_allgaus” is True.

Version 1.5.0 (2012/10/29):

Improved WCS handling. PyBDSM can now read images with a much greater variety of WCS systems (e.g., the VOPT spectral system).

Fixed a bug related to the use of a detection image when a subimage is specified (with trim_box).

Version 1.4.5 (2012/10/12):

Added option (incl_empty) to include empty islands (that have no un-flagged Gaussians) in output catalogs. Any such empty islands are given negative source IDs and have positions given by the location of the peak of the island.

Fixed a bug in Gaussian fitting that could cause a crash when fitting fails.

Fixed a bug in parallelization that could cause a crash due to improper concatenation of result lists.

Version 1.4.4 (2012/10/09):

Fixed a bug related to the parallelization of Gaussian fitting that could cause a crash due to improper mapping of island lists to processes.

Improved logging.

Added a warning when one or more islands are not fit (i.e., no valid, unflagged Gaussians were found).

Added code to handle images with no unblanked pixels.

Improved fitting robustness.

Version 1.4.3 (2012/10/04):

Fixed a bug in the mean map calculation that caused mean maps with constant values (i.e., non-2D maps) to have values of 0.0 Jy/beam unless mean_map = 'const' was explicitly specified.

Fixed a bug in the PSF vary module that resulted in incorrect PSF generators being used. Added an option to smooth the resulting PSF images (psf_smooth). Parallelized the PSF interpolation and smoothing steps. Improved PSF vary documentation.

Version 1.4.2 (2012/09/25):

Dramatically reduced time required to identify valid wavelet islands. Fixed bug that resulted in output FITS gaul tables being improperly sorted.

Version 1.4.1 (2012/09/11):

Added SAMP (Simple Application Messaging Protocol) support to the write_catalog, export_image, and show_fit tasks. These tasks can now use SAMP to communicate with other programs connected to a SAMP hub (e.g., ds9, Topcat, Aladin).

Version 1.4.0 (2012/09/11):

Parallelized Gaussian fitting, shapelet decomposition, validation of wavelet islands, and mean/rms map generation. The number of cores to be used can be specified with the ncores option (default is to use all).

Version 1.3.2 (2012/08/22):

Fixed a bug that could cause the user-specified rms_box value to be ignored. Added an option to enable the Monte Carlo error estimation for ‘M’-type sources (the do_mc_errors option), which is now disabled by default.

Version 1.3.1 (2012/07/11):

Fixed a bug that caused images written when output_all = True to be transposed. Added frequency information to all output images. Improved fitting robustness to prevent rare cases in which no valid Gaussians could be fit to an island. Modified the island-finding routine to handle NaNs properly.

Version 1.3.0 (2012/07/03):

Fixed a bug in the calculation of positional errors for Gaussians.

Adjusted rms_box algorithm to check for negative rms values (due to interpolation with cubic spline). If negative values are found, either the box size is increased or the interpolation is done with order=1 (bilinear) instead.

Output now includes the residual image produced using only wavelet Gaussians (if any) when atrous_do=True and output_all=True.

Improved organization of files when output_all=True.

Added logging of simple statistics (mean, std. dev, skew, and kurtosis) of the residual images.

Included image rotation (if any) in beam definition. Rotation angle can vary across the image (defined by image WCS).

Added Sagecal output format for Gaussian catalogs.

Added check for newer versions of the PyBDSM software tar.gz file available on ftp.strw.leidenuniv.nl.

Added total island flux (from sum of pixels) to gaul and srl catalogs.

Version 1.2 (2012/06/06):

Added option to output flux densities for every channel found by the spectral index module.

Added option to spectral index module to allow use of flux densities that do not meet the desired SNR.

Implemented an adaptive scaling scheme for the rms_box parameter that shrinks the box size near bright sources and expands it far from them (enabled with the adaptive_rms_box option when rms_box is None). This scheme generally results in improved rms and mean maps when both strong artifacts and extended sources are present.

Improved speed of Gaussian fitting to wavelet images.

Added option to calculate fluxes within a specified aperture radius in pixels (set with the aperture option). Aperture fluxes, if measured, are output in the srl format catalogs.

Version 1.1 (2012/03/28):

Tweaked settings that affect fitting of Gaussians to improve fitting in general.

Modified calculation of the rms_box parameter (when the rms_box option is None) to work better with fields composed mainly of point sources when strong artifacts are present.

Modified fitting of large islands to adopt an iterative fitting scheme that limits the number of Gaussians fit simultaneously per iteration to 10. This change speeds up fitting of large islands considerably.

Added the option to use a “detection” image for island detection (the detection_image option); source properties are still measured from the main input image. This option is particularly useful with images made with LOFAR’s AWImager, as the uncorrected, flat-noise image (the *.restored image) is better for source detection than the corrected image (the *.restored.corr image).

Modified the polarization module so that sources that appear only in Stokes Q or U (and hence not in Stokes I) are now identified. This identification is done using the polarized intensity (PI) image.

Improved the plotting speed (by a factor of many) in show_fit when there are a large number of islands present.

Simplified the spectral index module to make it more user friendly and stable.

Altered reading of images to correctly handle 4D cubes.

Extended the psf_vary module to include fitting of stacked PSFs with Gaussians, interpolation of the resulting parameters across the image, and correction of the deconvolved source sizes using the interpolated PSFs.

Added residual rms and mean values to source catalogs. These values can be compared to background rms and mean values as a quick check of fit quality.

Added output of shapelet parameters as FITS tables.

Fixed many minor bugs.

See the changelog (accessible from the interactive shell using help changelog) for details of all changes since the last version.

PyBDSF Basics

PyBDSF has been designed to share many similarities with the CASA interactive environment (known as casa [1] ), which is in turn based on AIPS. Therefore, the commands used in PyBDSF should be familiar to anyone who has used these software packages.

Starting PyBDSF

After installing (see Downloading and installing) you can start the interactive PyBDSF shell by simply opening a terminal and typing:

$ pybdsf

at the terminal prompt.

Note

If the above command does not work, please make sure you have installed the interactive shell (it is not installed by default). See Downloading and installing for details.

The interactive environment will then load, and a welcome screen listing common commands and tasks will be shown. You will then be at the PyBDSF prompt, which looks like this:

BDSF [1]:

Quitting PyBDSF

To quit PyBDSF, type quit or enter CNTL-D at the prompt.

Getting help

PyBDSF has an extensive built-in help system. To get help on a command or task, type:

help <command or task name>

For example, to get help on the process_image task, type:

help process_image

To get help on a input parameter to a task, type:

help '<parameter name>'

Note the quotes, which are necessary (since parameter names are strings). For example, to get help on the 'rms_box' paramter, type:

help 'rms_box'

Simply typing help will start the Python help system.

Logging

Logging of all task output is done automatically to a log file. Logs for subsequent runs on the same image are appended to the end of the log file. The log for each run includes a listing of all the non-default and internally derived parameters, so that a run can be easily reproduced using only information in the log.

Commands

As in CASA, PyBDSF uses a number of commands to list input parameters for tasks, to execute the tasks, etc. The PyBDSF commands are as follows:

inp task ............ : Set current task and list parameters
go .................. : Run the current task
default ............. : Set current task parameters to default values
tput ................ : Save parameter values
tget ................ : Load parameter values

inp

This command sets the current task (e.g., inp process_image) and lists the relevant parameters for that task. If entered without a task name, the parameters of the previously set task will be listed.

Note

At startup, the current task is set to the process_image task.

go

This command executes the current task.

default

This command resets all parameters for a task to their default values.

If a task name is given (e.g.,``default show_fit``), the parameters for that task are reset. If no task name is given, the parameters of the current task are reset.

tput

This command saves the processing parameters to a file.

Note

After the successful completion of a task, the current parameters are saved to the file ‘pybdsf.last’.

A file name may be given (e.g., tput 'savefile.sav'), in which case the parameters are saved to the file specified. If no file name is given, the parameters are saved to the file ‘pybdsf.last’. The saved parameters can be loaded using the tget command.

tget

This command loads the processing parameters from a parameter save file.

A file name may be given (e.g., tget 'savefile.sav'), in which case the parameters are loaded from the file specified. If no file name is given, the parameters are loaded from the file ‘pybdsf.last’ if it exists.

Normally, the save file is created by the tput command.

Tasks

The following tasks are available in PyBDSF:

process_image ....... : Process an image: find sources, etc.
show_fit ............ : Show the results of a fit
write_catalog ....... : Write out list of sources to a file
export_image ........ : Write residual/model/rms/mean image to a file

process_image: This task processes an image to find and measure sources. See process_image: processing an image for details.
show_fit: This task shows the result of a fit. See show_fit: visualizing the fit results for details.
write_catalog: This task writes the source catalog. See write_catalog: writing source catalogs for details.
export_image: This task exports an internally derived image. See export_image: exporting internally derived images for details.

Hierarchy of an astronomical image

The following figure shows the basic hierarchy of an image adopted by PyBDSF. Islands of emission are identified and decomposed into Gaussians. The Gaussians are then grouped into sources.

image hierarchy — The hierarchy of an image.

Quick-start example

Below is an example of using PyBDSF to find and measure sources in an image:

$ pybdsf
PyBDSF version 1.1
========================================================================
PyBDSF commands
  inp task ............ : Set current task and list parameters
  par = val ........... : Set a parameter (par = '' sets it to default)
                          Autocomplete (with TAB) works for par and val
  go .................. : Run the current task
  default ............. : Set current task parameters to default values
  tput ................ : Save parameter values
  tget ................ : Load parameter values
PyBDSF tasks
  process_image ....... : Process an image: find sources, etc.
  show_fit ............ : Show the results of a fit
  write_catalog ....... : Write out list of sources to a file
  export_image ........ : Write residual/model/rms/mean image to a file
PyBDSF help
  help command/task ... : Get help on a command or task
                          (e.g., help process_image)
  help 'par' .......... : Get help on a parameter (e.g., help 'rms_box')
  help changelog ...... : See list of recent changes
________________________________________________________________________

BDSF [1]: inp process_image
--------> inp(process_image)
PROCESS_IMAGE: Find and measure sources in an image.
=================================================================================
filename ................. '': Input image file name
advanced_opts ........ False : Show advanced options
adaptive_rms_box ..... False : Use adaptive rms_box when determining rms and
                               mean maps
atrous_do ............ False : Decompose Gaussian residual image into multiple
                               scales
beam .................. None : FWHM of restoring beam. Specify as (maj, min, pos
                               ang E of N) in degrees. E.g., beam = (0.06, 0.02,
                               13.3). None => get from header
flagging_opts ........ False : Show options for Gaussian flagging
frequency ............. None : Frequency in Hz of input image. E.g., frequency =
                               74e6. None => get from header.
interactive .......... False : Use interactive mode
mean_map .......... 'default': Background mean map: 'default' => calc whether to
                               use or not, 'zero' => 0, 'const' => clipped mean,
                               'map' => use 2-D map
multichan_opts ....... False : Show options for multi-channel images
output_opts .......... False : Show output options
polarisation_do ...... False : Find polarisation properties
psf_vary_do .......... False : Calculate PSF variation across image
rms_box ............... None : Box size, step size for rms/mean map calculation.
                               Specify as (box, step) in pixels. E.g., rms_box =
                               (40, 10) => box of 40x40 pixels, step of 10
                               pixels. None => calculate inside program
rms_map ............... None : Background rms map: True => use 2-D rms map; False
                               => use constant rms; None => calculate inside
                               program
shapelet_do .......... False : Decompose islands into shapelets
spectralindex_do ..... False : Calculate spectral indices (for multi-channel
                               image)
thresh ................ None : Type of thresholding: None => calculate inside
                               program, 'fdr' => use false detection rate
                               algorithm, 'hard' => use sigma clipping
thresh_isl ............. 3.0 : Threshold for the island boundary in number of
                               sigma above the mean. Determines extent of island
                               used for fitting
thresh_pix ............. 5.0 : Source detection threshold: threshold for the
                               island peak in number of sigma above the mean. If
                               false detection rate thresholding is used, this
                               value is ignored and thresh_pix is calculated
                               inside the program

BDSF [2]: filename = 'sb48.fits'
BDSF [3]: go
--------> go()
--> Opened 'sb48.fits'
Image size .............................. : (256, 256) pixels
Number of channels ...................... : 1
Beam shape (major, minor, pos angle) .... : (0.002916, 0.002654, -173.36) degrees
Frequency of averaged image ............. : 146.497 MHz
Blank pixels in the image ............... : 0 (0.0%)
Flux from sum of (non-blank) pixels ..... : 29.565 Jy
Derived rms_box (box size, step size) ... : (61, 20) pixels
--> Variation in rms image significant
--> Using 2D map for background rms
--> Variation in mean image significant
--> Using 2D map for background mean
Min/max values of background rms map .... : (0.05358, 0.25376) Jy/beam
Min/max values of background mean map ... : (-0.03656, 0.06190) Jy/beam
--> Expected 5-sigma-clipped false detection rate < fdr_ratio
--> Using sigma-clipping thresholding
Number of islands found ................. : 4
Fitting islands with Gaussians .......... : [====] 4/4
Total number of Gaussians fit to image .. : 12
Total flux in model ..................... : 27.336 Jy
Number of sources formed from Gaussians   : 6

BDSF [4]: show_fit
--------> show_fit()

The figure made by show_fit is shown in the figure below. In the plot window, one can zoom in, save the plot to a file, etc. The list of best-fit Gaussians found by PyBDSF may be written to a file for use in other programs as follows:

BDSF [5]: write_catalog
--------> write_catalog()
--> Wrote FITS file 'sb48.pybdsf.srl.fits'

The output Gaussian or source list contains source positions, fluxes, etc.

show_fit example output — Output of `show_fit`, showing the original image with and without sources, the model image, and the residual (original minus model) image. Boundaries of the islands of emission found by PyBDSF are shown in light blue. The fitted Gaussians are shown for each island as ellipses (the sizes of which correspond to the FWHMs of the Gaussians). Gaussians that have been grouped together into a source are shown with the same color. For example, the two red Gaussians of island #1 have been grouped together into one source, and the nine Gaussians of island #0 have been grouped into 4 separate sources.

Footnotes

`process_image`: processing an image

A standard analysis is performed using the process_image task. This task reads in the input image, calculates background rms and mean images, finds islands of emission, fits Gaussians to the islands, and groups the Gaussians into sources. Furthermore, the process_image task encompases a number of modules that allow decomposing an image into shapelets, calculating source spectral indices, deriving source polarization properties, and correcting for PSF variations across the image.

When process_image is executed, PyBDSF performs the following steps in order:

Reads in the image and collapses specific frequency channels with weights (see Multichannel options) and produces a ‘continuum’ image (the ch0 image) for all polarisations with which source detection is done.
Calculates basic statistics of the image and sensible values of the processing parameters. First, the number of beams per source is calculated (see Determining whether an image is confused for details), using a sensible estimate of box size and step size (which can be set using the rms_box parameter). Next, the thresholds are set. They can either be hard thresholded (by the user or set as 5-sigma for pixel threshold and 3-sigma for island boundaries by default) or can be calculated using the False Detection Rate (FDR) method using a user defined value for \(\alpha\) (the fdr_alpha parameter). If the user does not specify whether hard thresholding or FDR thresholding should be applied, one or the other is chosen internally based on the ratio of expected false pixels to true pixels.
Calculates the rms and mean images. The 3-sigma clipped rms and mean are calculated inside boxes of defined by the rms_box parameter. Optionally, these images can be calculated using adaptive scaling of this box, so that a smaller box (defined the the rms_box_bright parameter) is used near bright sources (where strong artifacts are more likely). Intermediate values are calculated using bicubic spline interpolation by default (the order of the spline interpolation can be set with the spline_rank parameter). Depending on the resulting statistics (see Determining whether an image is confused for details), we either adopt the rms image or a constant rms in the following analysis.
Identifies islands of contiguous emission. First all pixels greater than the pixel threshold are identified. Next, starting from each of these pixels, all contiguous pixels (defined by 8-connectivity, i.e., the surrounding eight pixels) higher than the island boundary threshold are identified as belonging to one island, accounting properly for overlaps of islands.
Fits multiple Gaussians to each island. The number of multiple Gaussians to be fit can be determined by three different methods (using the ini_gausfit parameter). With initial guesses corresponding to these peaks, Gaussians are simultaneously fit to the island using the Levenberg-Marqhardt algorithm. Sensible criteria for bad solutions are defined (see Flagging options). If multiple Gaussians are fit and one of them is a bad solution then the number of Gaussians is decreased by one and fit again, until all solutions in the island are good (or zero in number, in which case it’s flagged). After the final fit to the island, the deconvolved size is computed assuming the theoretical beam, and the statistics in the source area and in the island are computed and stored. Errors on each of the fitted parameters are computed using the formulae in Condon (1997) [1].
Groups nearby Gaussians within an island into sources. See Grouping of Gaussians into sources for details. Fluxes for the grouped Gaussians are summed to obtain the total flux of the source (the uncertainty is calculated by summing the Gaussian uncertainties in quadrature). The source position is set to be its centroid (the position of the maximum of the source is also calculated and output). The total source size is measured using moment analysis (see http://en.wikipedia.org/wiki/Image_moment for a nice overview of moment analysis). Errors on the source position and size are estimated using Condon (1997), but additional uncertainties due to uncertainties in the constituent Gaussians may be estimated using a Monte Carlo technique.
Continues with further processing, if the user has specified that one or more of the following modules be used:
- Shapelet decomposition (see Shapelet decomposition module for details)
- Wavelet decomposition (see À trous wavelet decomposition module for details)
- Estimation of PSF variation (see PSF variation module for details)
- Calculation of polarization properties (see Polarization module for details)
- Calculation of spectral indices (see Spectral index module for details)

General reduction parameters

Type inp process_image to list the main reduction parameters:

PROCESS_IMAGE: Find and measure sources in an image.
================================================================================
filename ................. '': Input image file name
adaptive_rms_box ..... False : Use adaptive rms_box when determining rms and
                               mean maps
advanced_opts ........ False : Show advanced options
atrous_do ............ False : Decompose Gaussian residual image into multiple
                               scales
beam .................. None : FWHM of restoring beam. Specify as (maj, min, pos
                               ang E of N) in degrees. E.g., beam = (0.06, 0.02,
                               13.3). None => get from header
flagging_opts ........ False : Show options for Gaussian flagging
frequency ............. None : Frequency in Hz of input image. E.g., frequency =
                               74e6. None => get from header. For more than one
                               channel, use the frequency_sp parameter.
interactive .......... False : Use interactive mode
mean_map .......... 'default': Background mean map: 'default' => calc whether to
                               use or not, 'zero' => 0, 'const' => clipped mean,
                               'map' => use 2-D map.
multichan_opts ....... False : Show options for multi-channel images
output_opts .......... False : Show output options
polarisation_do ...... False : Find polarisation properties
psf_vary_do .......... False : Calculate PSF variation across image
rms_box ............... None : Box size, step size for rms/mean map calculation.
                               Specify as (box, step) in pixels. E.g., rms_box =
                               (40, 10) => box of 40x40 pixels, step of 10
                               pixels. None => calculate inside program
rms_map ............... None : Background rms map: True => use 2-D rms map;
                               False => use constant rms; None => calculate
                               inside program
shapelet_do .......... False : Decompose islands into shapelets
spectralindex_do ..... False : Calculate spectral indices (for multi-channel
                               image)
thresh ................ None : Type of thresholding: None => calculate inside
                               program, 'fdr' => use false detection rate
                               algorithm, 'hard' => use sigma clipping
thresh_isl ............. 3.0 : Threshold for the island boundary in number of
                               sigma above the mean. Determines extent of
                               island used for fitting
thresh_pix ............. 5.0 : Source detection threshold: threshold for the
                               island peak in number of sigma above the mean. If
                               false detection rate thresholding is used, this
                               value is ignored and thresh_pix is calculated
                               inside the program

Each of the parameters is described in detail below.

filename

This parameter is a string (no default) that sets the input image file name. The input image can be a FITS or CASA 2-, 3-, or 4-D cube.

adaptive_rms_box

This parameter is a Boolean (default is False). If True, an adaptive box is used when calculating the rms and mean maps. See Adaptive box options for details of the options.

advanced_opts

This parameter is a Boolean (default is False). If True, the advanced options are shown. See Advanced options for details of the advanced options.

atrous_do

This parameter is a Boolean (default is False). If True, wavelet decomposition will be performed. See À trous wavelet decomposition module for details of the options.

beam

This parameter is a tuple (default is None) that defines the FWHM of restoring beam. Specify as (maj, min, pos ang E of N) in degrees. E.g., beam = (0.06, 0.02, 13.3). For more than one channel, use the beam_spectrum parameter. If the beam is not given by the user, then it is looked for in the image header. If not found, then an error is raised. PyBDSF will not work without knowledge of the restoring beam.

flagging_opts

This parameter is a Boolean (default is False). If True, the Gaussian flagging options will be listed. See Flagging options for details of the options.

frequency

This parameter is a float (default is None) that defines the frequency in Hz of the input image. E.g., frequency = 74e6. For more than one channel, use the frequency_sp parameter. If the frequency is not given by the user, then it is looked for in the image header. If not found, then an error is raised. PyBDSF will not work without knowledge of the frequency.

interactive

This parameter is a Boolean (default is False). If True, interactive mode is used. In interactive mode, plots are displayed at various stages of the processing so that the user may check the progress of the fit.

First, plots of the rms and mean background images are displayed along with the islands found, before fitting of Gaussians takes place. The user should verify that the islands and maps are reasonable before preceding.

Next, if atrous_do = True, the fits to each wavelet scale are shown. The wavelet fitting may be truncated at the current scale if desired.

Lastly, the final results are shown.

mean_map

This parameter is a string (default is 'default') that determines how the background mean map is computed and how it is used further.

If 'const', then the value of the clipped mean of the entire image (set by the kappa_clip option) is used as the background mean map.

If 'zero', then a value of zero is used.

If 'map', then the 2-dimensional mean map is computed and used. The resulting mean map is largely determined by the value of the rms_box parameter (see the rms_box parameter for more information).

If 'default', then PyBDSF will attempt to determine automatically whether to use a 2-dimensional map or a constant one as follows. First, the image is assumed to be confused if bmpersrc_th < 25 or the ratio of the clipped mean to rms (clipped mean/clipped rms) is > 0.1, else the image is not confused. Next, the mean map is checked to see if its spatial variation is significant. If so, then a 2-D map is used and, if not, then the mean map is set to either 0.0 or a constant depending on whether the image is thought to be confused or not.

Generally, 'default' works well. However, if there is significant extended emission in the image, it is often necessary to force the use of a constant mean map using either 'const' or 'mean'.

multichan_opts

This parameter is a Boolean (default is False). If True, the multichannel options will be listed. See Multichannel options for details of the options.

output_opts

This parameter is a Boolean (default is False). If True, the output options will be listed. See Output options for details of the options.

polarisation_do

This parameter is a Boolean (default is False). If True, polarization properties will be calculated for the sources. See Polarization module for details of the options.

psf_vary_do

This parameter is a Boolean (default is False). If True, the spatial variation of the PSF will be estimated and its effects corrected. See PSF variation module for details of the options.

rms_box

This parameter is a tuple (default is None) of two integers and is probably the most important input parameter for PyBDSF. The first integer, boxsize, is the size of the 2-D sliding box for calculating the rms and mean over the entire image. The second, stepsize, is the number of pixels by which this box is moved for the next measurement. If None, then suitable values are calculated internally.

In general, it is best to choose a box size that corresponds to the typical scale of artifacts in the image, such as those that are common around bright sources. Too small of a box size will effectively raise the local rms near a source so much that a source may not be fit at all; too large a box size can result in underestimates of the rms due to oversmoothing. A step size of 1/3 to 1/4 of the box size usually works well.

Note

The spline_rank parameter also affects the rms and mean maps. If you find ringing artifacts in the rms or mean maps near bright sources, try adjusting this parameter.

rms_map

This parameter is a Boolean (default is None). If True, then the 2-D background rms image is computed and used. If False, then a constant value is assumed (use rms_value to force the rms to a specific value). If None, then the 2-D rms image is calculated, and if the variation is statistically significant then it is taken, else a constant value is assumed. The rms image used for each channel in computing the spectral index follows what was done for the channel-collapsed image.

Generally, the default value works well. However, if there is significant extended emission in the image, it is often necessary to force the use of a constant rms map by setting rms_map = False.

shapelet_do

This parameter is a Boolean (default is False). If True, shapelet decomposition of the islands will be performed. See Shapelet decomposition module for details of the options.

spectralindex_do

This parameter is a Boolean (default is False). If True, spectral indices will be calculated for the sources. See Spectral index module for details of the options.

thresh

This parameter is a string (default is None). If thresh = 'hard', then a hard threshold is assumed, given by thresh_pix. If thresh = 'fdr', then the False Detection Rate algorithm of Hopkins et al. (2002) is used to calculate the value of thresh_pix. If thresh = None, then the false detection probability is first calculated, and if the number of false source pixels is more than fdr_ratio times the estimated number of true source pixels, then the 'fdr' threshold option is chosen, else the 'hard' threshold option is chosen.

thresh_isl

This parameter is a float (default is 3.0) that determines the region to which fitting is done. A higher value will produce smaller islands, and hence smaller regions that are considered in the fits. A lower value will produce larger islands. Use the thresh_pix parameter to set the detection threshold for sources. Generally, thresh_isl should be lower than thresh_pix.

Only regions above the absolute threshold will be used. The absolute threshold is calculated as abs_thr = mean + thresh_isl * rms. Use the mean_map and rms_map parameters to control the way the mean and rms are determined.

thresh_pix

This parameter is a float (default is 5.0) that sets the source detection threshold in number of sigma above the mean. If false detection rate thresholding is used, this value is ignored and thresh_pix is calculated inside the program

This parameter sets the overall detection threshold for islands (i.e. thresh_pix = 5 will find all sources with peak flux densities per beam of 5-sigma or greater). Use the thresh_isl parameter to control how much of each island is used in fitting. Generally, thresh_pix should be larger than thresh_isl.

Only islands with peaks above the absolute threshold will be used. The absolute threshold is calculated as abs_thr = mean + thresh_pix * rms. Use the mean_map and rms_map parameters to control the way the mean and rms are determined.

Adaptive box options

If adaptive_rms_box = True, the rms_box is reduced in size near bright sources and enlarged far from them. This scaling attempts to account for possible strong artifacts around bright sources while still acheiving accurate background rms and mean values when extended sources are present. This option is generally slower than non-adaptive scaling.

Use the rms_box parameter to set the large-scale box and the rms_box_bright parameter to set the small-scale box. The threshold for bright sources can be set with the adaptive_thresh parameter:

adaptive_rms_box ...... True : Use adaptive rms_box when determining rms and mean maps
  adaptive_thresh ..... None : Sources with pixels above adaptive_thresh*
                               clipped_rms will be considered as bright sources (i.e.,
                               with potential artifacts). None => calculate inside
                               program
  rms_box_bright ...... None : Box size, step size for rms/mean map
                               calculation near bright sources. Specify as (box, step)
                               in pixels. None => calculate inside program

adaptive_thresh

This parameter is a float (default is None) that sets the SNR above which sources may be affected by strong artifacts Sources that meet the SNR threshold will use the small-scale box (set by the rms_box_bright parameter) if their sizes at a threshold of 10.0 is less than 25 beam areas.

If None, the threshold is varied from 500 to 50 to attempt to obtain at least 5 candidate bright sources.

rms_box_bright

This parameter is a tuple (default is None) of two integers that sets the box and step sizes to use near bright sources (determined by the adaptive_thresh parameter). The large-scale box size is set with the rms_box parameter.

Advanced options

If advanced_opts = True, a number of additional options are listed. The advanced options do not usually need to be altered from the default values, but can be useful, for example, for fine tuning a fit or for quickly fitting a small region of a much larger image.

The advanced options are:

advanced_opts ......... True : Show advanced options
  aperture ............ None : Radius of aperture in pixels inside which aperture
                               fluxes are measured for each source. None => no aperture
                               fluxes measured
  aperture_posn .. 'centroid': Position the aperture (if aperture is not None) on: 'centroid' or
                               'peak' of the source.
  blank_limit ......... None : Limit in Jy/beam below which pixels are blanked. None => no such
                               blanking is done
  bmpersrc_th ......... None : Theoretical estimate of number of beams per
                               source. None => calculate inside program
  check_outsideuniv .. False : Check for pixels outside the universe
  detection_image ........ '': Detection image file name used only for
                               detecting islands of emission. Source
                               measurement is still done on the main image
  do_cache ........... False : Cache internally derived images to disk
  do_mc_errors ....... False : Estimate uncertainties for 'M'-type sources
                               using Monte Carlo method
  fdr_alpha ........... 0.05 : Alpha for FDR algorithm for thresholds
  fdr_ratio ............ 0.1 : For thresh = None; if #false_pix / #source_pix <
                               fdr_ratio, thresh = 'hard' else thresh = 'fdr'
  fittedimage_clip ..... 0.1 : Sigma for clipping Gaussians while creating fitted
                               image
  fix_to_beam ........ False : Fix major and minor axes and PA of Gaussians to beam?
  group_by_isl ....... False : Group all Gaussians in each island into a single
                               source
  group_method .. 'intensity': Group Gaussians into sources using 'intensity' map or
                               'curvature' map
  group_tol ............ 1.0 : Tolerance for grouping of Gaussians into sources:
                               larger values will result in larger sources
  ini_gausfit ..... 'default': Initial guess for Gaussian parameters: 'default',
                               'fbdsm', or 'nobeam'
  ini_method .... 'intensity': Method by which inital guess for fitting of Gaussians is chosen:
                               'intensity' or 'curvature'
  kappa_clip ........... 3.0 : Kappa for clipped mean and rms
  maxpix_isl .......... None : Maximum number of pixels with emission per island.
                               None -> no limit
  minpix_isl .......... None : Minimum number of pixels with emission per island.
                               None -> calculate inside program
  ncores .............. None : Number of cores to use during fitting, None => use
                               all
  peak_fit ............ True : Find and fit peaks of large islands before fitting
                               entire island
  peak_maxsize ........ 30.0 : If island size in beam area is more than this,
                               attempt to fit peaks separately (if
                               peak_fit=True). Min value is 30
  rmsmean_map_filename  None : Filenames of FITS files to use as the mean and rms maps,
                               given as a list [<mean_map.fits>, <rms_map.fits>]. If
                               supplied, the internally generated mean and rms maps
                               are not used
  rms_value ........... None : Value of constant rms in Jy/beam to use if rms_map
                               = False. None => calculate inside program
  spline_rank ............ 3 : Rank of the interpolating function for rms/mean
                               map
  split_isl ........... True : Split island if it is too large, has a large
                               convex deficiency and it opens well. If it doesn't
                               open well, then isl.mean = isl.clipped_mean, and
                               is taken for fitting. Splitting, if needed, is
                               always done for wavelet images
  splitisl_maxsize .... 50.0 : If island size in beam area is more than this,
                               consider splitting island. Min value is 50
  src_ra_dec .......... None : List of source positions at which fitting is done.  E.g.,
                               src_ra_dec = [(197.1932, 47.9188), (196.5573, 42.4852)].
  src_radius_pix ...... None : Radius of the island (if src_ra_dec is not None) in pixels. None
                               => radius is set to the FWHM of the beam major axis.
  stop_at ............. None : Stops after: 'isl' = island finding step or 'read'
                               = image reading step
  trim_box ............ None : Do source detection on only a part of the image.
                               Specify as (xmin, xmax, ymin, ymax) in pixels.
                               E.g., trim_box = (120, 840, 15, 895). None => use
                               entire image

aperture

This parameter is a float (default is None) that sets the radius (in pixels) inside which the aperture flux is measured for each source. The aperture is centered on the either the centroid or the peak of the source (depending on the value of the aperture_posn option). Errors are calculated from the mean of the rms map inside the aperture.

aperture_posn

This parameter is a string (default is 'centroid') that sets the how the aperture is positioned relative to the source. If ‘centroid’, the aperture is centered on the source centroid. If ‘peak’, the aperture is centered on the source peak. If aperture=None (i.e., no aperture radius is specified), this parameter is ignored.

blank_limit

This parameter is a float (default is None) that sets the limit in Jy/beam below which pixels are blanked. All pixels in the ch0 image with a value less than the specified limit and with at least 4 neighboring pixels with values also less than this limit are blanked. If None, any such pixels are left unblanked (and hence will affect the rms and mean maps, etc.). Pixels with a value of NaN are always blanked.

bmpersrc_th

This parameter is a float (default is None) that sets the theoretical estimate of number of beams per source. If None, the value is calculated as N/[n*(alpha-1)], where N is the total number of pixels in the image, n is the number of pixels in the image whose value is greater than 5 times the clipped rms, and alpha is the slope of the differential source counts distribution, assumed to be 2.5.

The value of bmpersrc_th is used to estimate the average separation in pixels between two sources, which in turn is used to estimate the boxsize for calculating the background rms and mean images. In addition, if the value is below 25 (or the ratio of clipped mean to clipped rms of the image is greater than 0.1), the image is assumed to be confused and hence the background mean is put to zero.

check_outsideuniv

This parameter is a Boolean (default is False). If True, then the coordinate of each pixel is examined to check if it is outside the universe, which may happen when, e.g., an all sky image is made with SIN projection (commonly done at LOFAR earlier). When found, these pixels are blanked (since imaging software do not do this on their own). Note that this process takes a lot of time, as every pixel is checked in case weird geometries and projections are used.

detection_image

This parameter is a string (default is '') that sets the detection image file name used only for detecting islands of emission. Source measurement is still done on the main image. The detection image can be a FITS or CASA 2-, 3-, or 4-D cube and must have the same size and WCS parameters as the main image.

do_cache

This parameter is a Boolean (default is False) that controls whether internally derived images are stored in memory or are cached to disk. Caching can reduce the amount of memory used, and is therefore useful when analyzing large images.

do_mc_errors

This parameter is a Boolean (default is False). If True, uncertainties on the sizes and positions of ‘M’-type sources due to uncertainties in the constituent Gaussians are estimated using a Monte Carlo technique. These uncertainties are added in quadrature with those calculated using Condon (1997). If False, these uncertainties are ignored, and errors are calculated using Condon (1997) only.

Enabling this option will result in longer run times if many ‘M’-type sources are present, but should give better estimates of the uncertainites, particularly for complex sources composed of many Gaussians.

fdr_alpha

This parameter is a float (default is 0.05) that sets the value of alpha for the FDR algorithm for thresholding. If thresh is 'fdr', then the estimate of fdr_alpha (see Hopkins et al. 2002 [2] for details) is stored in this parameter.

fdr_ratio

This parameter is a float (default is 0.1). When thresh = None, if #false_pix / #source_pix < fdr_ratio, thresh = 'hard' otherwise thresh = 'fdr'.

fittedimage_clip

This parameter is a float (default is 0.1). When the residual image is being made after Gaussian decomposition, the model images for each fitted Gaussian are constructed up to a size 2b, such that the amplitude of the Gaussian falls to a value of fitted_image_clip times the local rms, b pixels from the peak.

fix_to_beam

This parameter is a Boolean (default is False). If True, then during fitting the major and minor axes and PA of the Gaussians are fixed to the beam. Only the amplitude and position are fit. If False, all parameters are fit.

group_by_isl

This parameter is a Boolean (default is False). If True, all Gaussians in the island belong to a single source. If False, grouping is controlled by the group_tol parameter.

group_method

This parameter is a string (default is 'intensity'). Gaussians are deemed to be a part of the same source if: (1) no pixel on the line joining the centers of any pair of Gaussians has a (Gaussian-reconstructed) value less than the island threshold, and (2) the centers are separated by a distance less than half the sum of their FWHMs along the line joining them. If 'curvature', the above comparisons are done on the curature map (see Hancock et al. 2012). If 'intensity', the comparisons are done on the intensity map.

group_tol

This parameter is a float (default is 1.0) that sets the tolerance for grouping of Gaussians into sources: larger values will result in larger sources. Sources are created by grouping nearby Gaussians as follows: (1) If the difference between the minimum value between two Gaussians and the lower of the peak flux densities of the Gaussians in an island is less than group_tol * thresh_isl * rms_clip, and (2) if the centers are separated by a distance less than 0.5 * group_tol of the sum of their FWHMs along the PA of the line joining them, they belong to the same island.

ini_gausfit

This parameter is a string (default is 'default'). These are three different ways of estimating the initial guess for fitting of Gaussians to an island of emission. If 'default', the maximum number of Gaussians is estimated from the number of peaks in the island. An initial guess is made for the parameters of these Gaussians before final fitting is done. This method should produce the best results when there are no large sources present. If 'simple', the maximum number of Gaussians per island is set to 25, and no initial guess for the Gaussian parameters is made. Lastly, the 'nobeam' method is similar to the 'default' method, but no information about the beam is used. This method is best used when source sizes are expected to be very different from the beam and is generally slower than the other methods. For wavelet images, the value used for the original image is used for wavelet order j <= 3 and 'nobeam' for higher orders.

ini_method

This parameter is a string (default is 'intensity'). If 'intensity', the inital guess described in the help for the ini_gausfit parameter is calculated using the intensity (ch0) image. If 'curvature', it is done using the curvature map (see Hancock et al. 2012).

kappa_clip

This parameter is a float (default is 3.0) that is the factor used for Kappa-alpha clipping, as in AIPS. For an image with few source pixels added on to (Gaussian) noise pixels, the dispersion of the underlying noise will need to be determined. This is done iteratively, whereby the actual dispersion is first computed. Then, all pixels whose value exceeds kappa clip times this rms are excluded and the rms is computed again. This process is repeated until no more pixels are excluded. For well behaved noise statistics, this process will converge to the true noise rms with a value for this parameter ~3-5. A large fraction of source pixels, less number of pixels in total, or significant non-Gaussianity of the underlying noise will all lead to non-convergence.

maxpix_isl

This parameter is an integer (default is None) that sets the maximum number of pixels in an island for the island to be included. If None, there is no limit.

minpix_isl

This parameter is an integer (default is None) that sets the minimum number of pixels in an island for the island to be included. If None, the number of pixels is set to 1/3 of the area of an unresolved source using the beam and pixel size information in the image header. It is set to 6 pixels for all wavelet images.

ncores

This parameter is an integer (default is None) that sets the number of cores to use during fitting. If None, all available cores are used (one core is reserved for plotting).

peak_fit

This parameter is a Boolean (default is True). When True, PyBDSF will identify and fit peaks of emission in large islands iteratively (the size of islands for which peak fitting is done is controlled with the peak_maxsize option), using a maximum of 10 Gaussians per iteration. Enabling this option will generally speed up fitting (by factors of many for large islands), but may result in somewhat higher residuals.

peak_maxsize

This parameter is a float (default is 30.0). If island size in beam area is more than this value, attempt to fit peaks iteratively (if peak_fit = True). The minimum value is 30.

rmsmean_map_filename

This parameter is a list (default is None) that sets the filenames of FITS files to use as the mean and rms maps, given as a list [<mean_map.fits>, <rms_map.fits>]. If supplied, the internally generated mean and rms maps are not used.

rms_value

This parameter is a float (default is None) that sets the value of constant rms in Jy/beam to use if rms_map = False. If None, the value is calculated inside the program.

spline_rank

This parameter is an integer (default is 3) that sets the order of the interpolating spline function to interpolate the background rms and mean maps over the entire image.

Note

Bicubic interpolation (the default) can cause ringing artifacts to appear in the rms and mean maps in regions where sharp changes occur. These artifacts can result in regions with negative values. If you find such artifacts, try changing the spline_rank parameter.

split_isl

This parameter is a Boolean (default is True). If True, an island is split if it is too large, has a large convex deficiency and it opens well. If it doesn’t open well, then isl.mean = isl.clipped_mean, and is taken for fitting. Splitting, if needed, is always done for wavelet images

splitisl_maxsize

This parameter is a float (default is 50.0). If island size in beam area is more than this, consider splitting island. The minimum value is 50.

src_ra_dec

This parameter is a list of tuples (default is None) that defines the center positions at which fitting will be done. The size of the region used for the fit is given by the src_radius_pix parameter. Positions should be given as a list of RA and Dec, in degrees, one set per source. These positions will override the normal island finding module.

src_radius_pix

This parameter is a float (default is None) that determines the size of the region used to fit the source positions specified by the src_ra_dec parameter. If None, the radius is set to the FWHM of the beam major axis.

stop_at

This parameter is a string (default is None) that stops an analysis after: ‘isl’ = island finding step or ‘read’ = image reading step.

trim_box

This parameter is a tuple (default is None) that defines a subregion of the image on which to do source detection. It is specified as (xmin, xmax, ymin, ymax) in pixels. E.g., trim_box = (120, 840, 15, 895). If None, the entire image is used.

Flagging options

If flagging_opts = True, a number of options are listed for flagging unwanted Gaussians that occur durring a fit. Flagged Gaussians are not included in any further analysis or catalog. They may be viewed using the show_fit task (see show_fit: visualizing the fit results). A flag value is associated with each flagged Gaussian that allows the user to determine the reason or reasons that it was flagged. If multiple flagging conditions are met by a single Gaussian, the flag values are summed. For example, if a Gaussian is flagged because it is too large (its size exceeds that implied by flag_maxsize_bm, giving a flag value of 64) and because it is too bright (its peak flux density per beam exceeds that implied by flag_maxsnr, giving a flag value of 2) then the final flag value is 64 + 2 = 66.

Note

If a fit did not produce good results, it is often useful to check whether there are flagged Gaussians and adjust the flagging options as necessary. Flagged Gaussians can be viewed by setting ch0_flagged = True in the show_fit task.

The options for flagging of Gaussians are:

flagging_opts ......... True : Show options for Gaussian flagging
  flag_bordersize ........ 0 : Flag Gaussian if centre is outside border -
                               flag_bordersize pixels
  flag_maxsize_bm ..... 25.0 : Flag Gaussian if area greater than flag_maxsize_bm
                               times beam area
  flag_maxsize_isl ..... 1.0 : Flag Gaussian if x, y bounding box around
                               sigma-contour is factor times island bbox
  flag_maxsnr .......... 1.5 : Flag Gaussian if peak is greater than flag_maxsnr
                               times max value in island
  flag_minsize_bm ...... 0.7 : Flag Gaussian if flag_smallsrc = True and area
                               smaller than flag_minsize_bm times beam area
  flag_minsnr .......... 0.9 : Flag Gaussian if peak is less than flag_minsnr
                               times thresh_pix times local rms
  flag_smallsrc ...... False : Flag sources smaller than flag_minsize_bm times
                               beam area

flag_bordersize

This parameter is an integer (default is 0). Any fitted Gaussian whose centre is flag_bordersize pixels outside the island bounding box is flagged. The flag value is increased by 4 (for x) and 8 (for y).

flag_maxsize_bm

This parameter is a float (default is 25.0). Any fitted Gaussian whose size is greater than flag_maxsize_bm times the synthesized beam is flagged. The flag value is increased by 64.

flag_maxsize_fwhm

This parameter is a float (default is 0.3). Any fitted Gaussian whose contour of flag_maxsize_fwhm times the FWHM falls outside the island is flagged. The flag value is increased by 256.

flag_maxsize_isl

This parameter is a float (default is 1.0). Any fitted Gaussian whose maximum x-dimension is larger than flag_maxsize_isl times the x-dimension of the island (and likewise for the y-dimension) is flagged. The flag value is increased by 16 (for x) and 32 (for y).

flag_maxsnr

This parameter is a float (default is 1.5). Any fitted Gaussian whose peak is greater than flag_maxsnr times the value of the image at the peak of the Gaussian is flagged. The flag value is increased by 2.

flag_minsize_bm

This parameter is a float (default is 0.7). If flag_smallsrc is True, then any fitted Gaussian whose size is less than flag_maxsize_bm times the synthesized beam is flagged. The Gaussian flag is increased by 128.

flag_minsnr

This parameter is a float (default is 0.7). Any fitted Gaussian whose peak is less than flag_minsnr times thresh_pix times the local rms is flagged. The flag value is increased by 1.

flag_smallsrc

This parameter is a Boolean (default is False). If True, then fitted Gaussians whose size is less than flag_minsize_bm times the synthesized beam area are flagged. When combining Gaussians into sources, an error is raised if a 2x2 box with the peak of the Gaussian does not have all four pixels belonging to the source. Usually this means that the Gaussian is an artifact or has a very small size.

If False, then if either of the sizes of the fitted Gaussian is zero, then the Gaussian is flagged.

If the image is barely Nyquist sampled, this flag is best set to False. This flag is automatically set to False while decomposing wavelet images into Gaussians.

Output options

If output_opts = True, options to control the output generated by process_image are listed. By default, only a log file is generated and output is controlled with the export_image (see export_image: exporting internally derived images) and write_catalog (see write_catalog: writing source catalogs) tasks. However, the user can specify that a number of optional output files be made automatically whenever process_image is run. These options are most useful for debugging or when running PyBDSF non-interactively in a Python script (see Scripting PyBDSF).

The output options are:

output_opts ........... True : Show output options
  bbs_patches ......... None : For BBS format, type of patch to use: None => no
                               patches. 'single' => all Gaussians in one patch.
                               'gaussian' => each Gaussian gets its own patch.
                               'source' => all Gaussians belonging to a single
                               source are grouped into one patch. 'mask' => use mask
                               file specified by bbs_patches_mask
  bbs_patches_mask .... None : Name of the mask file (of same size as input image)
                               that defines the patches if bbs_patches = 'mask'
  indir ............... None : Directory of input FITS files. None => get from
                               filename
  opdir_overwrite .. 'overwrite': 'overwrite'/'append': If output_all=True,
                               delete existing files or append a new directory
  outdir .............. None : Directory to use for all output files
                               (including log files). None => parent directory of the
                               input filename
  output_all ......... False : Write out all files automatically to directory
                               'outdir/filename_pybdsf'
  plot_allgaus ....... False : Make a plot of all Gaussians at the end
  plot_islands ....... False : Make separate plots of each island during fitting
                               (for large images, this may take a long time and a
                               lot of memory)
  print_timing ....... False : Print basic timing information
  quiet .............. False : Suppress text output to screen. Output is still
                               sent to the log file as usual
  savefits_meanim .... False : Save background mean image as fits file
  savefits_normim .... False : Save norm image as fits file
  savefits_rankim .... False : Save island rank image as fits file
  savefits_residim ... False : Save residual image as fits file
  savefits_rmsim ..... False : Save background rms image as fits file
  solnname ............ None : Name of the run, to be appended to the name of the
                               output directory
  verbose_fitting .... False : Print out extra information during fitting

bbs_patches

This parameter is a string (default is None) that sets the type of patch to use in BBS-formatted catalogs. When the Gaussian catalogue is written as a BBS-readable sky file, this option determines whether all Gaussians are in a single patch ('single'), there are no patches (None), all Gaussians for a given source are in a separate patch ('source'), each Gaussian gets its own patch ('gaussian'), or a mask image is used to define the patches ('mask').

If you wish to have patches defined by island, then set group_by_isl = True before fitting to force all Gaussians in an island to be in a single source. Then set bbs_patches = 'source' when writing the catalog.

bbs_patches_mask

This parameter is a string (default is None) that sets the file name of the mask file to use to define patches in BBS-formatted catalogs. The mask image should be 1 inside the patches and 0 elsewhere and should be the same size as the input image (before any trim_box is applied). Any Gaussians that fall outside of the patches will be ignored and will not appear in the output sky model.

indir

This parameter is a string (default is None) that sets the directory of input FITS files. If None, the directory is defined by the input filename.

opdir_overwrite

This parameter is a string (default is 'overwrite') that determines whether existing output files are overwritten or not.

outdir

This parameter is a string (default is None) that sets the directory to use for all output files (including log files). If None, the parent directory of the input image filename is used.

output_all

This parameter is a Boolean (default is False). If True, all output products are written automatically to the directory 'outdir/filename_pybdsf'.

plot_allgaus

This parameter is a Boolean (default is False). If True, make a plot of all Gaussians at the end.

plot_islands

This parameter is a Boolean (default is False). If True, make separate plots of each island during fitting (for large images, this may take a long time and a lot of memory).

print_timing

This parameter is a Boolean (default is False). If True, print basic timing information.

quiet

This parameter is a Boolean (default is False). If True, suppress text output to screen. Output is still sent to the log file as usual.

savefits_meanim

This parameter is a Boolean (default is False). If True, save background mean image as a FITS file.

savefits_normim

This parameter is a Boolean (default is False). If True, save norm image as a FITS file.

savefits_rankim

This parameter is a Boolean (default is False). If True, save island rank image as a FITS file.

savefits_residim

This parameter is a Boolean (default is False). If True, save residual image as a FITS file.

savefits_rmsim

This parameter is a Boolean (default is False). If True, save background rms image as a FITS file.

solnname

This parameter is a string (default is None) that sets the name of the run, to be appended to the name of the output directory.

verbose_fitting

This parameter is a Boolean (default is False). If True, print out extra information during fitting.

Multichannel options

If multichan_opts = True, the options used to control the way PyBDSF handles images with more than one frequency channel are listed. In particular, these options control how the multichannel image is collapsed to form the ch0 image on which source detection is done.

The options concerning multichannel images are:

multichan_opts ........ True : Show options for multi-channel images
  beam_sp_derive ...... True : If True and beam_spectrum is None, then assume
                               header beam is for lowest frequency and scales
                               with frequency for channels
  beam_spectrum ....... None : FWHM of synthesized beam per channel. Specify as
                               [(bmaj_ch1, bmin_ch1, bpa_ch1), (bmaj_ch2,
                               bmin_ch2, bpa_ch2), etc.] in degrees. E.g.,
                               beam_spectrum = [(0.01, 0.01, 45.0), (0.02, 0.01,
                               34.0)] for two channels. None => all equal to beam
  collapse_av ........... [] : List of channels to average if collapse_mode =
                               'average'; None => all
  collapse_ch0 ........... 0 : Number of the channel for source extraction, if
                               collapse_mode = 'single'
  collapse_mode ... 'average': Collapse method: 'average' or 'single'. Average
                               channels or take single channel to perform source
                               detection on
  collapse_wt ....... 'unity': Weighting: 'unity' or 'rms'. Average channels with
                               weights = 1 or 1/rms_clip^2 if collapse_mode =
                               'average'
  frequency_sp ........ None : Frequency in Hz of channels in input image when
                               more than one channel is present. E.g., frequency
                               = [74e6, 153e6]. None => get from header

beam_sp_derive

This parameter is a Boolean (default is True). If True and the parameter beam_spectrum is None and no channel-dependent beam parameters are found in the header, then we assume that the beam in the header is for the lowest frequency of the image cube and scale accordingly to calculate the beam per channel. If False, then a constant value of the beam (that given in the header or specified with the beam parameter) is taken instead.

beam_spectrum

This parameter is a list of tuples (default is None) that sets the FWHM of synthesized beam per channel. Specify as [(bmaj_ch1, bmin_ch1, bpa_ch1), (bmaj_ch2, bmin_ch2, bpa_ch2), etc.] in degrees. E.g., beam_spectrum = [(0.01, 0.01, 45.0), (0.02, 0.01, 34.0)] for two channels.

If None, then the channel-dependent restoring beam is searched for in the header. If not found, the beam is either assumed to be a constant or to scale with frequency, depending on whether the parameter beam_sp_derive is False or True.

collapse_av

This parameter is a list of integers (default is []) that specifies the channels to be averaged to produce the continuum image for performing source extraction, if collapse_mode is 'average'. If the value is [], then all channels are used. Otherwise, the value is a Python list of channel numbers.

collapse_ch0

This parameter is an integer (default is 0) that specifies the number of the channel for source extraction, if collapse_mode = 'single'.

collapse_mode

This parameter is a string (default is 'average') that determines whether, when multiple channels are present, the source extraction is done on a single channel ('single') or an average of many channels ('average').

collapse_wt

This parameter is a string (default is 'unity'). When collapse_mode is 'average', then if this value is 'unity', the channels given by collapse_av are averaged with unit weights and if 'rms', then they are averaged with weights which are inverse square of the clipped rms of each channel image.

frequency_sp

This parameter is a list of floats (default is None) that sets the frequency in Hz of channels in input image when more than one channel is present. E.g., frequency_sp = [74e6, 153e6].

If the frequency is not given by the user, then it is looked for in the image header. If not found, then an error is raised. PyBDSF will not work without the knowledge of the frequency.

À trous wavelet decomposition module

If atrous_do = True, this module decomposes the residual image that results from the normal fitting of Gaussians into wavelet images of various scales. Such a decomposition is useful if there is extended emission that is not well fit during normal fitting. Such emission therefore remains in the Gaussian residual image and can be further fit by Gaussians whose size is tuned to the various wavelet scales. Therefore, wavelet decomposition should be used when there is significant residual emission that remains after normal Gaussian fitting.

The wavelet module performs the following steps:

The number of wavelet scales to be considered is set by the atrous_jmax parameter. By default, this number is determined automatically from the size of the largest island in the image. Wavelet images are then made for scales of order (j) ranging from 1 to jmax.
For each scale (j), the appropriate à trous wavelet transformation is made (see Holschneider et al. 1989 for details). Additionally, the “remainder” image (called the c_J image) is also made. This image includes all emission not included in the other wavelet images.
Depending on the value of the atrous_sum option, fitting is done to either an image that is a sum over all scales equal to or larger than the scale under consideration (atrous_sum = True) or to an image of a single scale (atrous_sum = False). Fitting to the sum over all larger scales will generally result in increased signal to noise.
If atrous_bdsm = True, an rms map and a mean map are made for each wavelet image and Gaussians are fit in the normal way. Gaussians can be optionally restricted to lie within islands found from the initial image. If a wavelet island overlaps spatially with an existing island, the two islands are merged together to form a single island. The wavelet Gaussians can then be included in source catalogs (see write_catalog: writing source catalogs).

The options for this module are as follows:

atrous_do ............. True : Decompose Gaussian residual image into multiple
                               scales
  atrous_bdsm_do ...... True : Perform source extraction on each wavelet scale
  atrous_jmax ............ 0 : Max allowed wavelength order, 0 => calculate
                               inside program
  atrous_lpf ........... 'b3': Low pass filter, either 'b3' or 'tr', for B3
                               spline or Triangle
  atrous_orig_isl .... False : Restrict wavelet Gaussians to islands found in
                               original image
  atrous_sum .......... True : Fit to the sum of images of the remaining wavelet
                               scales
  use_scipy_fft ....... True : Use fast SciPy FFT for convolution

atrous_bdsm_do

This parameter is a Boolean (default is False). If True, PyBDSF performs source extraction on each wavelet scale.

atrous_jmax

This parameter is an integer (default is 0) which is the maximum order of the à trous wavelet decomposition. If 0 (or <0 or >15), then the value is determined within the program. The value of this parameter is then estimated as the (lower) rounded off value of ln[(nm-l)/(l-1) + 1]/ln2 + 1 where nm is the minimum of the residual image size (n, m) in pixels and l is the length of the filter à trous lpf (see the atrous_lpf parameter for more info).

A sensible value is such that the size of the kernel is not more than 3-4 times smaller than the smallest image dimension.

atrous_lpf

This parameter is a string (default is 'b3') that sets the low pass filter, which can be either the B3 spline or the triangle function, which is used to generate the à trous wavelets. The B3 spline is [1, 4, 6, 4, 1] and the triangle is [1, 2, 1], normalised so that the sum is unity. The lengths of the filters are hence 5 and 3 respectively.

atrous_orig_isl

This parameter is a Boolean (default is False). If True, all wavelet Gaussians must lie within the boundaries of islands found in the original image. If False, new islands that are found only in the wavelet images are included in the final fit.

atrous_sum

This parameter is a Boolean (default is True). If True, fitting is done on an image that is the sum of the remaining wavelet scales. Using the sum will generally result in improved signal. If False, fitting is done on only the wavelet scale under consideration.

use_scipy_fft

This parameter is a Boolean (default is True). If True, the SciPy FFT function will be used instead of the custom version. The SciPy version is much faster but also uses much more memory.

PSF variation module

If psf_vary_do = True, then the spatial variations in the PSF are estimated and their effects corrected for. To this end, PyBDSF performs the following steps:

A list of sources that are likely to be unresolved is constructed. This is done by first selecting only type ‘S’ sources by default (see Definition of output columns for details of source types), but this restriction can be overridden using the psf_stype_only option) and sources with SNRs that exceed psf_snrcut. Next, a function is fit to determine how the size of sources (normalized by the median size) varies with the SNR. The function used is defined as \(\sigma / median = \sqrt(c_1^2 + c_2^2/SNR^2)\), where \(\sigma\) is the size of the Gaussian and \(c_1\) and \(c_2\) are free parameters. Clipping of outliers is done during this fitting, controlled by the psf_nsig parameter. Lastly, unresolved sources are selected by choosing sources that lie within psf_kappa2 times the rms of this best-fit sigma-SNR relation. As this last step can be unreliable for high-SNR sources, an additional selection can be made for the highest SNR sources using the psf_high_snr parameter. All sources with SNRs above psf_high_snr will be taken as unresolved.
Next the image is tessellated using Voronoi tessellation to produce tiles within which the PSF shape is calculated (and assumed to be constant). The list of probable unresolved sources is filtered to select “calibrator” sources to use to determine the tessellation tiles. These sources are the brightest sources (known as the primary generators), defined as those sources that have SNRs in the top fraction of sources defined by psf_snrtop and that also have SNRs greater than psf_snrcutstack. These sources are then grouped by their proximity, if they are within 50% of the distance to third closest source.
The unresolved sources within each tile that have SNRs greater than psf_snrcutstack are then stacked to form a high-SNR PSF. For each tile, this PSF is fit with a Gaussian to recover its size. The significance of the variation in the sizes across the image is quantified.
If the variation is significant, the major axis, minor axis, and position angle are then interpolated across the image. Smoothing can be applied to these images to smooth out artifacts due to noise and the interpolation. Additionally, images are made of the ratio of peak-to-total flux and peak-to-aperture flux (if an aperture is specified). These ratio images provide conversions from total flux to peak flux for point sources. In the absence of smearing effects, these ratios should be around unity. However, if ionospheric effects are present, significant smearing can be present. In this case, these ratio images can be useful, for example, in determining the sensitivity at a particular location in the image to a point source with a given total flux.
Lastly, the deconvolved source sizes are adjusted to include the PSF variation as a function of position.

The options for this module are as follows:

psf_vary_do ........... True : Calculate PSF variation across image
  psf_high_snr ........ None : SNR above which all sources are taken to be
                               unresolved. E.g., psf_high_snr = 20.0. None => no
                               such selection is made
  psf_itess_method ....... 0 : 0 = normal, 1 = 0 + round, 2 = LogSNR, 3 =
                               SqrtLogSNR
  psf_kappa2 ........... 2.0 : Kappa for clipping for analytic fit
  psf_nsig ............. 3.0 : Kappa for clipping within each bin
  psf_over ............... 2 : Factor of nyquist sample for binning bmaj, etc. vs
                               SNR
  psf_smooth .......... None : Size of Gaussian to use for smoothing of
                               interpolated images in arcsec. None => no smoothing
  psf_snrcut .......... 10.0 : Minimum SNR for statistics
  psf_snrcutstack ..... 15.0 : Unresolved sources with higher SNR taken for
                               stacked psfs
  psf_snrtop .......... 0.15 : Fraction of SNR > snrcut as primary generators
  psf_stype_only ...... True : Restrict sources used in PSF variation
                               estimating to be only of type 'S'

psf_high_snr: This parameter is a float (default is None). Gaussians with SNR greater than this are used to determine the PSF variation, even if they are deemed to be resolved. This corrects for the unreliability at high SNRs in the algorithm used to find unresolved sources. The minimum value is 20.0. If None, then no such selection is made.
psf_itess_method: This parameter is an integer (default is 0) which can be 0, 1, 2 or 3, which corresponds to a tessellation method. If 0, 2 or 3, then the weights used for Voronoi tessellation are unity, log(SNR) and sqrt[log(SNR)] where SNR is the signal to noise ratio of the generator in a tile. If 1, then the image is tessellated such that each tile has smooth boundaries instead of straight lines, using pixel-dependent weights.
psf_kappa2: This parameter is a float (default is 2.0). When iteratively arriving at a statistically probable set of ‘unresolved’ sources, the fitted major and minor axis sizes versus SNR are binned and fitted with analytical functions. Those Gaussians which are within psf_kappa2 times the fitted rms from the fitted median are then considered ‘unresolved’ and are used further to estimate the PSFs.
psf_nsig: This parameter is a float (default is 3.0). When constructing a set of ‘unresolved’ sources for psf estimation, the (clipped) median, rms and mean of major and minor axis sizes of Gaussians versus SNR within each bin is calculated using kappa = psf_nsig.
psf_over: This parameter is an integer (default is 2). When constructing a set of ‘unresolved’ sources for psf estimation, this parameter controls the factor of nyquist sample for binning bmaj, etc. vs SNR.
psf_smooth: This parameter is a float (default is None) that sets the smoothing scale (in arcsec) used to smooth the interpolated images. Generally, artifacts due to noise and the interpolation can be significantly reduced if the smoothing scale is similar to the typical source separation scale.
psf_snrcut: This parameter is a float (default is 10.0). Only Gaussians with SNR greater than this are considered for processing. The minimum value is 5.0
psf_snrcutstack: This parameter is a float (default is 15.0). Only Gaussians with SNR greater than this are used for estimating PSF images in each tile.
psf_snrtop: This parameter is a float (default is 0.15). If psf_generators is ‘calibrator’, then the peak pixels of Gaussians which are the psf_snrtop fraction of the SNR distribution are taken as Voronoi generators.
psf_stype_only: This parameter is a Boolean (default is False). If True, sources are restricted to be only of type ‘S’.

Spectral index module

If spectralindex_do = True (and the input image has more than one frequency), then spectral indices are calculated for the sources in the following way:

The rms maps for the remaining channels are determined.
Neighboring channels are averaged to attempt to obtain the target SNR per channel for a given source, set by the specind_snr parameter.

Note

No color corrections are applied during averaging. However, unless the source spectral index is very steep or the channels are very wide, the correction is minimal. See Effect of neglecting color corrections for details.
Flux densities are measured for both individual Gaussians and for total sources. Once source flux densities have been measured in each channel, the SEDs are fit with the following power-law model: \(S_\nu \propto \nu^\alpha\), where \(\alpha\) is the spectral index. The resulting best-fit spectral index is then included in any catalogs that are written out (see write_catalog: writing source catalogs). In addition, plots of the fits can be viewed with the show_fit task (see show_fit: visualizing the fit results).

The options for this module are as follows:

spectralindex_do ...... True : Calculate spectral indices (for multi-channel
                               image)
  flagchan_rms ........ True : Flag channels before (averaging and) extracting
                               spectral index, if their rms if more than 5
                               (clipped) sigma outside the median rms over all
                               channels, but only if <= 10% of channels
  flagchan_snr ........ True : Flag channels that do not meet SNR criterion set
                               by specind_snr
  specind_maxchan ........ 0 : Maximum number of channels to average for a
                               given source when when attempting to meet target
                               SNR. 1 => no averaging; 0 => no maximum
  specind_snr .......... 3.0 : Target SNR to use when fitting power law. If
                               there is insufficient SNR, neighboring channels
                               are averaged to obtain the target SNR

flagchan_rms: This parameter is a Boolean (default is True). If True, then the clipped rms and median (r and m) of the clipped rms of each channel is calculated. Those channels whose clipped rms is greater than 4r away from m are flagged prior to averaging and calculating spectral indices from the image cube. However, these channels are flagged only if the total number of these bad channels does not exceed 10% of the total number of channels themselves.
flagchan_snr: This parameter is a Boolean (default is True). If True, then flux densities in channels that do not meet the target SNR are not used in fitting.
specind_maxchan: This parameter is an integer (default is 0) that sets the maximum number of channels that can be averaged together to attempt to reach the target SNR set by the specind_snr parameter. If 0, there is no limit to the number of channels that can be averaged. If 1, no averaging will be done.
specind_snr: This parameter is a float (default is 3.0) that sets the target SNR to use when fitting for the spectral index. If there is insufficient SNR, neighboring channels are averaged to obtain the target SNR. The maximum allowable number of channels to average is determined by the specind_maxchan parameter. Channels (after averaging) that fail to meet the target SNR are not used in fitting.

Polarization module

If polarisation_do = True, then the polarization properties of the sources are calculated. First, if pi_fit = True, source detection is performed on the polarized intensity (PI) image [3] to detect sources without Stokes I counterparts. The polarization module then calculates the I, Q, U, and V flux densities, the total, linear, and circular polarisation fractions and the linear polarisation angle of each Gaussian and source. The linear polarisation angle is defined from North, with positive angles towards East. Flux densities are calculated by fitting the normalization of the Gaussians found using the Stokes I or PI images.

For linearly polarised emission, the signal and noise add vectorially, giving a Rice distribution instead of a Gaussian one. To correct for this, a bias is estimated and removed from the polarisation fraction using the same method used for the NVSS catalog (see ftp://ftp.cv.nrao.edu/pub/nvss/catalog.ps). Errors on the linear and total polarisation fractions and polarisation angle are estimated using the debiased polarised flux density and standard error propagation. See Sparks & Axon (1999) [4] for a more detailed treatment.

The options for this module are as follows:

polarisation_do ....... True : Find polarisation properties
  pi_fit .............. True : Check the polarized intesity (PI) image for
                               sources not found in Stokes I
  pi_thresh_isl ....... None : Threshold for PI island boundary in number
                               of sigma above the mean. None => use thresh_isl
  pi_thresh_pix ....... None : Source detection threshold for PI image:
                               threshold for the island peak in number of sigma
                               above the mean. None => use thresh_pix

pi_fit: This parameter is a Boolean (default is True). If True, the polarized intensity image is searched for sources not present in the Stokes I image. If any such sources are found, they are added to the the Stokes I source lists. Use the pi_thresh_pix and pi_thresh_isl parameters to control island detection in the PI image.
pi_thresh_isl: This parameter is a float (default is None) that determines the region to which fitting is done in the polarized intensity (PI) image. If None, the value is set to that of the thresh_isl parameter. A higher value will produce smaller islands, and hence smaller regions that are considered in the fits. A lower value will produce larger islands. Use the pi_thresh_pix parameter to set the detection threshold for sources. Generally, pi_thresh_isl should be lower than pi_thresh_pix.
pi_thresh_pix: This parameter is a float (default is None) that sets the overall detection threshold for islands in the polarized intensity (PI) image (i.e. pi_thresh_pix = 5 will find all sources with peak flux densities per beam of 5-sigma or greater). If None, the value is set to that of the thresh_pix parameter. Use the pi_thresh_isl parameter to control how much of each island is used in fitting. Generally, pi_thresh_pix should be larger than pi_thresh_isl.

Shapelet decomposition module

If shapelet_do = True, then islands are decomposed into shapelets. Shapelets are a set of 2-D basis functions (for details, see Refregier 2003 [5]) that can be used to completely model any source, typically with far fewer parameters than pixels in the source. Shapelets are useful in particular for modeling complex islands that are not well modeled by Gaussians alone. PyBDSF can currently fit cartesian shapelets to an image. The shapelet parameters can be written to a catalog using write_catalog (see write_catalog: writing source catalogs).

For each island of emission, a shapelet decomposition is done after estimating the best values of the center, the scale \(\beta\), and nmax in the following way. First, an initial guess of \(\beta\) is taken as \(2\sqrt{[m2(x)m2(y)]}\), where \(m2\) is the second moment over the island, based on shapeelt analysis of simulated images of resolved sources. Similarly, a guess for nmax is taken as the minimum of 14, and maximum of 10 and \(2n + 2\) where \(n=\sqrt{(n^2 + m^2)}/n_p^n - 1\), where (n, m) is the size of the island and \(n^m_p\) is the synthesized beam minor axis FWHM in pixels. This guess for nmax is based partly on simulations and partly on the requirememts of computing time, number of constraints, etc, for shapelet decomposition.

These initial values are then used to calculate the optimal central position around which to decompose the island. First, for every pixel in the island, the coefficients c12 and c21 are computed assuming that pixel as the centre of expansion. Next, the zero crossings for every vertical line of the c12 image and horizontal line of the c21 image are computed. The intersection point of these two zero-crossing vectors is then taken as the proper centre of the expansion for the image. If this procedure does not work, then the first moment is taken as the center.

This updated center position is used to compute the optimal \(\beta\), which is taken as the value of \(\beta\) that minimises the residual rms in the island area. Using this \(\beta\), the center is computed once more and the final shapelet deocmposition is then made.

The options for this module are as follows:

shapelet_do ........... True : Decompose islands into shapelets
  shapelet_basis .. 'cartesian': Basis set for shapelet decomposition:
                               'cartesian' or 'polar'
  shapelet_fitmode .... 'fit': Calculate shapelet coeff's by fitting ('fit') or
                               integrating (None)
  shapelet_gresid ..... 'False': Use Gaussian residual image for shapelet
                               decomposition?

shapelet_basis: This parameter is a string (default is 'cartesian') that determines the type of shapelet basis used. Currently however, only cartesian is supported.
shapelet_fitmode: This parameter is a string (default is 'fit') that determines the method of calculating shapelet coefficients. If None, then these are calculated by integrating (actually, by summing over pixels, which introduces errors due to discretisation). If ‘fit’, then the coefficients are found by least-squares fitting of the shapelet basis functions to the image.
shapelet_gresid: This parameter is a Boolean (default is True). If True, then the shapelet decomposition is done on the Gaussian residual image rather that the ch0 image.

Footnotes

`show_fit`: visualizing the fit results

PyBDSF includes a task named show_fit that allows the user to quickly check the results of the process_image task. Use inp show_fit to list the parameters:

SHOW_FIT: Show results of fit.
================================================================================
broadcast ............ False : Broadcast Gaussian and source IDs and coordinates
                               to SAMP hub when a Gaussian is clicked?
ch0_flagged .......... False : Show the ch0 image with flagged Gaussians (if
                               any) overplotted
ch0_image ............. True : Show the ch0 image. This is the image used for
                               source detection
ch0_islands ........... True : Show the ch0 image with islands and Gaussians (if
                               any) overplotted
gmodel_image .......... True : Show the Gaussian model image
gresid_image .......... True : Show the Gaussian residual image
mean_image ............ True : Show the background mean image
pi_image ............. False : Show the polarized intensity image
psf_major ............ False : Show the PSF major axis variation
psf_minor ............ False : Show the PSF minor axis variation
psf_pa ............... False : Show the PSF position angle variation
rms_image ............. True : Show the background rms image
smodel_image ......... False : Show the shapelet model image
source_seds .......... False : Plot the source SEDs and best-fit spectral
                               indices (if image was processed with
                               spectralindex_do = True). Sources may be chosen
                               by ID with the 'c' key or, if ch0_islands = True,
                               by picking a source with the mouse
sresid_image ......... False : Show the shapelet residual image

Each of the parameters is described in detail below.

broadcast: This parameter is a Boolean (default is False) that determines whether the Gaussian and source IDs and coordinates are sent to a running SAMP Hub when a Gaussian is clicked on. Note that for the IDs to be useful, a catalog must have been sent to the SAMP hub previously using the write_catalog task (with outfile = 'SAMP').
ch0_flagged: This parameter is a Boolean (default is False) that determines whether to plot the ch0 image (the image used for source detection) with any flagged Gaussians overplotted.
ch0_image: This parameter is a Boolean (default is True) that determines whether to plot the ch0 image (the image used for source detection).
ch0_islands: This parameter is a Boolean (default is True) that determines whether to plot the ch0 image (the image used for source detection) with islands and Gaussians overplotted.
gmodel_image: This parameter is a Boolean (default is True) that determines whether to plot the Gaussian model image.
gresid_image: This parameter is a Boolean (default is True) that determines whether to plot the Gaussian residual image.
mean_image: This parameter is a Boolean (default is True) that determines whether to plot the background mean image.
pi_image: This parameter is a Boolean (default is False) that determines whether to plot the polarized intensity image.
psf_major: This parameter is a Boolean (default is False) that determines whether to plot the variation of the major axis of the PSF.
psf_minor: This parameter is a Boolean (default is False) that determines whether to plot the variation of the minor axis of the PSF.
psf_pa: This parameter is a Boolean (default is False) that determines whether to plot the variation of the position angle of the PSF.
rms_image: This parameter is a Boolean (default is True) that determines whether to plot the background rms image.
smodel_image: This parameter is a Boolean (default is False) that determines whether to plot the shapelet model image.
source_seds: This parameter is a Boolean (default is False) that determines whether to plot the source SEDs and best-fit spectral indices.
sresid_image: This parameter is a Boolean (default is False) that determines whether to plot the shapelet residual image.

`export_image`: exporting internally derived images

Internally derived images (e.g, the Gaussian model image) can be exported to FITS or CASA files using the export_image task:

EXPORT_IMAGE: Write one or more images to a file.
================================================================================
outfile ............... None : Output file name. None => file is named
                               automatically; 'SAMP' => send to SAMP hub (e.g., to
                               TOPCAT, ds9, or Aladin)
clobber .............. False : Overwrite existing file?
img_format ........... 'fits': Format of output image: 'fits' or 'casa'
img_type ....... 'gaus_resid': Type of image to export: 'gaus_resid',
                               'shap_resid', 'rms', 'mean', 'gaus_model',
                               'shap_model', 'ch0', 'pi', 'psf_major', 'psf_minor',
                               'psf_pa', 'psf_ratio', 'psf_ratio_aper', 'island_mask'
mask_dilation ............ 0 : Number of iterations to use for island-mask dilation.
                               0 => no dilation
pad_image ............ False : Pad image (with zeros) to original size

Each of the parameters is described in detail below.

outfile

This parameter is a string (default is None) that sets the name of the output file. If None, the file is named automatically. If ‘SAMP’ the image is sent to a running SAMP Hub (e.g., to ds9 or Aladin).

clobber

This parameter is a Boolean (default is False) that determines whether existing files are overwritten or not.

img_format

This parameter is a string (default is 'fits') that sets the output file format: 'fits' - FITS format, 'casa' - CASA format (requires casacore).

img_type

This parameter is a string (default is 'gaus_resid') that sets the type of image to export. The following images can be exported:

'ch0' - image used for source detection
'rms' - rms map image
'mean' - mean map image
'pi' - polarized intensity image
'gaus_resid' - Gaussian model residual image
'gaus_model' - Gaussian model image
'shap_resid' - Shapelet model residual image
'shap_model' - Shapelet model image
'psf_major' - image of major axis FWHM variation (arcsec)
'psf_minor' - image of minor axis FWHM variation (arcsec)
'psf_pa' - image of position angle variation (degrees east of north)
'psf_ratio' - image of peak-to-total flux variation (1/beam)
'psf_ratio_aper' - image of peak-to-aperture flux variation (1/beam)
'island_mask' - mask of islands (0 = outside island, 1 = inside island)

mask_dilation

This parameter is an integer (default is 0) that sets the number of dilation iterations to use when making the island mask. More iterations implies larger masked regions (one iteration expands the size of features in the mask by one pixel in all directions).

pad_image

This parameter is a Boolean (default is False) that determines whether the output image is padded to be the same size as the original image (without any trimming defined by the trim_box parameter). If False, the output image will have the size specified by the trim_box parameter.

`write_catalog`: writing source catalogs

The properties of the fitted Gaussians, sources, and shapelets may be written in a number of formats to a file using the write_catalog task. See below (Definition of output columns) for a detailed description of the output columns.

Note

For BBS and SAGECAL formats, the output catalogs always use the J2000 equinox. If the input image does not have an equinox of J2000, the coordinates of sources will be precessed to J2000. Catalogs in other formats will have the equinox of the image.

The task parameters are as follows:

WRITE_CATALOG: Write the Gaussian, source, or shapelet list to a file.
================================================================================
outfile ............... None : Output file name. None => file is named
                               automatically; 'SAMP' => send to SAMP hub (e.g., to
                               TOPCAT, ds9, or Aladin)
bbs_patches ........... None : For BBS format, type of patch to use: None => no
                               patches. 'single' => all Gaussians in one patch.
                               'gaussian' => each Gaussian gets its own patch.
                               'source' => all Gaussians belonging to a single
                               source are grouped into one patch. 'mask' => use mask
                               file specified by bbs_patches_mask
bbs_patches_mask ...... None : Name of the mask file (of same size as input image)
                               that defines the patches if bbs_patches = 'mask'
catalog_type .......... 'srl': Type of catalog to write:  'gaul' - Gaussian
                               list, 'srl' - source list (formed by grouping
                               Gaussians), 'shap' - shapelet list (FITS
                               format only)
clobber .............. False : Overwrite existing file?
correct_proj .......... True : Correct source parameters for image projection
                               (BBS format only)?
format ............... 'fits': Format of output Gaussian list: 'bbs', 'ds9',
                               'fits', 'star', 'kvis', 'ascii', 'csv', 'casabox', or
                               'sagecal'
incl_chan ............ False : Include fluxes from each channel (if any)?
incl_empty ........... False : Include islands without any valid Gaussians
                               (source list only)?
srcroot ............... None : Root name for entries in the output catalog. None
                               => use image file name

Each of the parameters is described in detail below.

outfile

This parameter is a string (default is None) that sets the name of the output file. If None, the file is named automatically. If ‘SAMP’ the full catalog (i.e., format = 'fits') is sent to a running SAMP Hub (e.g., to TOPCAT or Aladin).

bbs_patches

This parameter is a string (default is None) that sets the type of patch to use in BBS-formatted catalogs. When the Gaussian catalogue is written as a BBS-readable sky file, this option determines whether all Gaussians are in a single patch ('single'), there are no patches (None), all Gaussians for a given source are in a separate patch ('source'), each Gaussian gets its own patch ('gaussian'), or a mask image is used to define the patches ('mask').

If you wish to have patches defined by island, then set group_by_isl = True before fitting to force all Gaussians in an island to be in a single source. Then set bbs_patches = 'source' when writing the catalog.

bbs_patches_mask

This parameter is a string (default is None) that sets the file name of the mask file to use to define patches in BBS-formatted catalogs. The mask image should be 1 inside the patches and 0 elsewhere and should be the same size as the input image (before any trim_box is applied). Any Gaussians that fall outside of the patches will be ignored and will not appear in the output sky model.

catalog_type

This parameter is a string (default is 'srl') that sets the type of catalog to write: 'gaul' - Gaussian list, 'srl' - source list (formed by grouping Gaussians), 'shap' - shapelet list ('fits' format only)

Note

The choice of 'srl' or 'gaul' depends on whether you want all the source structure in your catalog or not. For example, if you are making a sky model for use as a model in calibration, you want to include all the source structure in your model, so you would use a Gaussian list ('gaul'), which writes each Gaussian. On the other hand, if you want to compare to other source catalogs, you want instead the total source flux densities, so use source lists ('srl'). For example, say you have a source that is unresolved in WENSS, but is resolved in your image into two nearby Gaussians that are grouped into a single source. In this case, you want to compare the sum of the Gaussians to the WENSS flux density, and hence should use a source list.

clobber

This parameter is a Boolean (default is False) that determines whether existing files are overwritten or not.

correct_proj

This parameter is a Boolean (default is True) that determines whether the source parameters in the output catalog will be corrected for first-order projection effects. If False, no correction is done. In this case, the position angle is relative to the +y axis, NOT true north, and source sizes are calculated assuming a constant pixel scale (equal to the scale at the image center).

If True, the position angle and source size are corrected using the average pixel size and angle offset (between the +y axis and north) at the location of the source center.

format

This parameter is a string (default is 'fits') that sets the format of the output catalog. The following formats are supported:

'bbs' - BlackBoard Selfcal sky model format (Gaussian list only)
'ds9' - ds9 region format
'fits' - FITS catalog format, readable by many software packages, including IDL, TOPCAT, Python, fv, Aladin, etc.
'star' - AIPS STAR format (Gaussian list only)
'kvis' - kvis format (Gaussian list only)
'ascii' - simple text file with spaces separating the values
'csv' - Comma-separated Values (CSV) text file
'casabox' - CASA region file (boxes only)
'sagecal' - SAGECAL sky model format (Gaussian list only)

Catalogues with the 'fits', 'ascii', and 'csv' formats include all available information (see Definition of output columns for column definitions). The other formats include only a subset of the full information.

incl_chan

This parameter is a Boolean (default is False) that determines whether the total flux densities of each source measured in each channel by the spectral index module are included in the output.

incl_empty

This parameter is a Boolean (default is False) that determines whether islands without any valid Gaussians are included in the output catalog. This option is only available for source lists. If True, islands for which Gaussian fitting failed will be included in the output catalog. In these cases, the source IDs are negative and only a subset of the standard columns will be populated (columns requiring information from Gaussian fits are left blank).

srcroot

This parameter is a string (default is None) that sets the root for source names in the output catalog.

Definition of output columns

The information included in the Gaussian and source catalogs varies by format and can include the following quantities.

Note

For ACSII, CSV, and FITS formats, the reference frequency (in Hz) and equinox are stored in the header of the catalog. The header in ASCII and CSV catalogs is the first few lines of the catalog. For FITS catalogs, this information is stored in the comments as well as in the FREQ0 and EQUINOX keywords in the primary header.

The shapelet catalog contains the following additional columns:

shapelet_basis: the basis coordinate system: ‘c’ for cartesian, ‘s’ for spherical
shapelet_beta: the \(\beta\) parameter of the shapelet decomposition
shapelet_nmax: the maximum order of the shapelet
shapelet_cf: a (flattened) array of the shapelet coefficients

Scripting PyBDSF

Because PyBDSF is written in Python, it is straightforward to use PyBDSF non-interactively within Python scripts (for example, to automate source detection in a large number of images for which the optimal analysis parameters are known). To use PyBDSF in a Python script, import it by calling:

import bdsf

inside your script.

Processing may then be done using process_image() as follows:

img = bdsf.process_image(filename, <args>)

where filename is the name of the image (in FITS or CASA format) or PyBDSF parameter save file and <args> is a comma-separated list of arguments defined as in the interactive environment (e.g., beam = (0.033, 0.033, 0.0), rms_map=False). See process_image: processing an image for details.

Note

The filename of the input image is also stored in the parameter save file. If you wish to override this filename (e.g., to use the saved parameters on a different image), give the save file as the first parameter and then explicitly set the filename. For example: img = bdsf.process_image('image1_savefile.sav', filename='image2.fits').

If the fit is successful, PyBDSF will return an Image object (in this example named img) which contains the results of the fit (among many other things).

When run in a Python script, it may be desirable to set output_all = True to write all output, including source lists, residual images, etc. to a directory named filename_pybdsf. Optionally, the same tasks used in the interactive PyBDSF shell are available for examining the fit and writing out the source list, residual image, etc. These tasks are methods of the Image object returned by bdsf.process_image() and are described below. The input parameters to each of these tasks are the same as those available in the interactive shell (see the relevant task section for details).

img.show_fit(): This method shows a quick summary of the fit by plotting the input image with the islands and Gaussians found, along with the model and residual images. See show_fit: visualizing the fit results for details.
img.export_image(): Write an internally derived image (e.g., the model image) to a FITS file. See export_image: exporting internally derived images for details.
img.write_catalog(): This method writes the Gaussian or source list to a file. See write_catalog: writing source catalogs for details.

An example of using PyBDSF within a Python script is given in Scripting example.

Alphabetical listing of all parameters

For a listing of all parameters, please see the Index (Index).

Simple image with point sources

Below is an example of running PyBDSF on an image composed primarily of point sources (a VLSS image).

$ pybdsf

PyBDSF version 1.7.0
========================================================================
PyBDSF commands
  inp task ............ : Set current task and list parameters
  par = val ........... : Set a parameter (par = '' sets it to default)
                          Autocomplete (with TAB) works for par and val
  go .................. : Run the current task
  default ............. : Set current task parameters to default values
  tput ................ : Save parameter values
  tget ................ : Load parameter values
PyBDSF tasks
  process_image ....... : Process an image: find sources, etc.
  show_fit ............ : Show the results of a fit
  write_catalog ....... : Write out list of sources to a file
  export_image ........ : Write residual/model/rms/mean image to a file
PyBDSF help
  help command/task ... : Get help on a command or task
                          (e.g., help process_image)
  help 'par' .......... : Get help on a parameter (e.g., help 'rms_box')
  help changelog ...... : See list of recent changes
________________________________________________________________________


BDSF [1]: filename='VLSS.fits'

Note

When PyBDSF starts up, the process_image task is automatically set to be the current task, so one does not need to set it with inp process_image.

BDSF [2]: frequency=74e6

Note

For this image, no frequency information was present in the image header, so the frequency must be specified manually.

BDSF [3]: interactive=T

Note

It is often advisable to use the interactive mode when processing an image for the first time. This mode will display the islands that PyBDSF has found before proceeding to fitting, allowing the user to check that they are reasonable.

BDSF [4]: go
---------> go()
--> Opened 'VLSS.fits'
Image size .............................. : (1024, 1024) pixels
Number of channels ...................... : 1
Number of Stokes parameters ............. : 1
Beam shape (major, minor, pos angle) .... : (0.02222, 0.02222, 0.0) degrees
Frequency of image ...................... : 74.000 MHz
Number of blank pixels .................. : 0 (0.0%)
Flux from sum of (non-blank) pixels ..... : 177.465 Jy
Derived rms_box (box size, step size) ... : (196, 65) pixels
--> Variation in rms image significant
--> Using 2D map for background rms
--> Variation in mean image significant
--> Using 2D map for background mean
Min/max values of background rms map .... : (0.06305, 0.16508) Jy/beam
Min/max values of background mean map ... : (-0.01967, 0.01714) Jy/beam
--> Expected 5-sigma-clipped false detection rate < fdr_ratio
--> Using sigma-clipping thresholding
Minimum number of pixels per island ..... : 5
Number of islands found ................. : 115
--> Displaying islands and rms image...
========================================================================
NOTE -- With the mouse pointer in plot window:
  Press "i" ........ : Get integrated fluxes and mean rms values
                       for the visible portion of the image
  Press "m" ........ : Change min and max scaling values
  Press "n" ........ : Show / hide island IDs
  Press "0" ........ : Reset scaling to default
  Click Gaussian ... : Print Gaussian and source IDs (zoom_rect mode,
                       toggled with the "zoom" button and indicated in
                       the lower right corner, must be off)
________________________________________________________________________

Note

At this point, because interactive=True, PyBDSF plots the islands. Once the plot window is closed, PyBDSF prompts the user to continue or to quit fitting:

Press enter to continue or 'q' to quit .. :
Fitting islands with Gaussians .......... : [==========================================] 115/115
Total number of Gaussians fit to image .. : 147
Total flux in model ..................... : 211.800 Jy
Number of sources formed from Gaussians   : 117

The process_image task has now finished. PyBDSF estimated a reasonable value for the rms_box parameter and determined that 2-D rms and mean maps were required to model the background of the image. Straightforward island thresholding at the 5-sigma level was used, and the minimum island size was set at 5 pixels. In total 115 islands were found, and 147 Gaussians were fit to these islands. These 147 Gaussians were then grouped into 117 sources. To check the fit, call the show_fit task:

BDSF [5]: show_fit
---------> show_fit()
========================================================================
NOTE -- With the mouse pointer in plot window:
  Press "i" ........ : Get integrated fluxes and mean rms values
                       for the visible portion of the image
  Press "m" ........ : Change min and max scaling values
  Press "n" ........ : Show / hide island IDs
  Press "0" ........ : Reset scaling to default
  Click Gaussian ... : Print Gaussian and source IDs (zoom_rect mode,
                       toggled with the "zoom" button and indicated in
                       the lower right corner, must be off)
________________________________________________________________________

The show_fit task produces the figure below. It is clear that the fit worked well and all significant sources were identified and modeled successfully.

example output — Example fit with default parameters of an image with mostly point sources.

Lastly, the plot window is closed, and the source catalog is written out to an ASCII file with the write_catalog task:

BDSF [6]: inp write_catalog
--------> inp(write_catalog)
WRITE_CATALOG: Write the Gaussian, source, or shapelet list to a file.
================================================================================
outfile ............... None : Output file name. None => file is named
                               automatically; 'SAMP' => send to SAMP hub (e.g.,
                               to TOPCAT, ds9, or Aladin)
bbs_patches ........... None : For BBS format, type of patch to use: None => no
                               patches. 'single' => all Gaussians in one patch.
                               'gaussian' => each Gaussian gets its own patch.
                               'source' => all Gaussians belonging to a single
                               source are grouped into one patch
bbs_patches_mask ...... None : Name of the mask file (of same size as input image)
                               that defines the patches if bbs_patches = 'mask'
catalog_type .......... 'srl': Type of catalog to write:  'gaul' - Gaussian
                               list, 'srl' - source list (formed by grouping
                               Gaussians), 'shap' - shapelet list
clobber .............. False : Overwrite existing file?
correct_proj .......... True : Correct source parameters for image projection
                               (BBS format only)?
format ............... 'fits': Format of output catalog: 'bbs', 'ds9', 'fits',
                               'star', 'kvis', or 'ascii', 'csv', 'casabox',
                               or 'sagecal'
incl_chan ............ False : Include flux densities from each channel (if any)?
incl_empty ........... False : Include islands without any valid Gaussians (source
                               list only)?
srcroot ............... None : Root name for entries in the output catalog. None
                               => use image file name

BDSF [7]: format='ascii'

BDSF [8]: go
---------> go()
--> Wrote ASCII file 'VLSS.fits.pybdsf.srl'

Image with artifacts

Occasionally, an analysis run with the default parameters does not produce good results. For example, if there are significant deconvolution artifacts in the image, the thresh_isl, thresh_pix, or rms_box parameters might need to be changed to prevent PyBDSF from fitting Gaussians to such artifacts. An example of running PyBDSF with the default parameters on such an image is shown in the figures below.

It is clear that a number of spurious sources are being detected. Simply raising the threshold for island detection (using the thresh_pix parameter) would remove these sources but would also remove many real but faint sources in regions of low rms. Instead, by setting the rms_box parameter to better match the typical scale over which the artifacts vary significantly, one obtains much better results. In this example, the scale of the regions affected by artifacts is approximately 20 pixels, whereas PyBDSF used a rms_box of 63 pixels when run with the default parameters, resulting in an rms map that is over-smoothed. Therefore, one should set rms_box=(20,10) so that the rms map is computed using a box of 20 pixels in size with a step size of 10 pixels (i.e., the box is moved across the image in 10-pixel steps). See the figures below for a summary of the results of this call.

Image with extended emission

If there is extended emission that fills a significant portion of the image, the background rms map will likely be biased high in regions where extended emission is present, affecting the island determination (this can be checked during a run by setting interactive=True). Setting rms_map=False and mean_map='const' or 'zero' will force PyBDSF to use a constant mean and rms value across the whole image. Additionally, setting flag_maxsize_bm to a large value (50 to 100) will allow large Gaussians to be fit, and setting atrous_do=True will fit Gaussians of various scales to the residual image to recover extended emission missed in the standard fitting. Depending on the source structure, the thresh_isl and thresh_pix parameters may also have to be adjusted as well to ensure that PyBDSF finds and fits islands of emission properly. An example analysis of an image with significant extended emission is shown below. Note that large, complex sources can require a long time to fit (on the order of hours).

Scripting example

You can use the complete functionality of PyBDSF within Python scripts (see Scripting PyBDSF for details). Scripting can be useful, for example, if you have a large number of images or if PyBDSF needs to be called as part of an automated reduction. Below is a short example of using PyBDSF to find sources in a number of images automatically. In this example, the best reduction parameters were determined beforehand for a representative image and saved to a PyBDSF save file using the tput command (see Commands for details).

# pybdsf_example.py
#
# This script fits a number of images automatically, writing out source
# catalogs and residual and model images for each input image. Call it
# with "python pybdsf_example.py"

import bdsf

# Define the list of images to process and the parameter save file
input_images = ['a2597.fits', 'a2256_1.fits', 'a2256_2.fits',
                 'a2256_3.fits', 'a2256_4.fits', 'a2256_5.fits']
save_file = 'a2256.sav'

# Now loop over the input images and process them
for input_image in input_images:

    if input_image == 'a2597.fits':
        # For this one image, run with different parameters.
        # Note that the image name is the first argument to
        # process_image:
        img = bdsf.process_image(input_image, rms_box=(100,20))

    else:
        # For the other images, use the 'a2256.sav` parameter save file.
        # The quiet argument is used to supress output to the terminal
        # (it still goes to the log file).
        # Note: when a save file is used, it must be given first in the
        # call to process_image:
        img = bdsf.process_image(save_file, filename=input_image, quiet=True)

    # Write the source list catalog. File is named automatically.
    img.write_catalog(format='fits', catalog_type='srl')

    # Write the residual image. File is named automatically.
    img.export_image(img_type='gaus_resid')

    # Write the model image. Filename is specified explicitly.
    img.export_image(img_type='gaus_model', outfile=input_image+'.model')

Using SAMP interoperability

PyBDSF supports SAMP (Simple Application Messaging Protocol) to provide interoperability to other applications, such as TOPCAT [1], ds9 [2], and Aladin [3]. To use this functionality, a SAMP hub must be running (both TOPCAT and Aladin come with SAMP hubs). Below is an example of using PyBDSF with TOPCAT. In this example, it is assumed that an image has already been processed with process_image.

BDSF [1]: process_image('VLSS.fits')
...

At this point, make sure that TOPCAT is started and its SAMP hub is running (activated by clicking the “Attempt to connect to SAMP hub” icon in the lower right-hand corner and selecting “Start internal hub”). Next, we send the PyBDSF source list to TOPCAT with write_catalog:

BDSF [2]: inp write_catalog

BDSF [3]: outfile='SAMP'

BDSF [4]: go
---------> go()
--> Table sent to SAMP hub.

TOPCAT should automatically load the table. Double-click on the table name in TOPCAT to open the table viewer. We can use now the show_fit task to highlight the table row that corresponds to a source of interest. To do this, we start show_fit with broadcast = True:

BDSF [6]: show_fit(broadcast=T)
========================================================================
NOTE -- With the mouse pointer in plot window:
  Press "i" ........ : Get integrated flux densities and mean rms
                       values for the visible portion of the image
  Press "m" ........ : Change min and max scaling values
  Press "n" ........ : Show / hide island IDs
  Press "0" ........ : Reset scaling to default
  Click Gaussian ... : Print Gaussian and source IDs (zoom_rect mode,
                       toggled with the "zoom" button and indicated in
                       the lower right corner, must be off)
________________________________________________________________________

Now, clicking on a Gaussian will highlight the row corresponding to the source to which the Gaussian belongs. Gaussian catalogs (i.e., made with catalog_type='gaul' in write_catalog) are also supported (and may be used simultaneously in TOPCAT with source catalogs).

Images can be sent to ds9 or Aladin using the export_image task in the same way (with outfile = 'SAMP'). Furthermore, if an image was sent, clicking on a Gaussian in the show_fit window will tell ds9 or Aladin to center their view on the coordinates of the Gaussian’s center.

Footnotes

Determining whether an image is confused

The number of beams per source (if not set with the bmpersrc_th parameter) is calculated by assuming the number of sources in the image, \(N_s\), as:

\[N_s = (\text{No. pixels} > 5\sigma)/(<\text{pix/src}>),\]

where the average number of pixels per source, \(<pix/src>\), is given by:

\[2\pi \sigma_{\text{major}} \sigma_{\text{minor}} \times (\ln(S_{\text{min}}/5\sigma) - 1/(\alpha - 1)),\]

where \(\alpha\) is the slope of the differential source counts taken from Katgert et al. (1988) [1]. Assuming a minimum of one pixel to define a source and ignoring the effect of noise for sources close to the threshold, we can ignore the logarithmic term and hence \(\text{bmpersrc_th} = (n\times m)/(\text{No. pixels} > 5\sigma)/(\alpha-1))\). The value of bmpersrc_th is used to decide whether the image is expected to be confused (and if so, the mean image is taken to be zero) and also to estimate the box size for calculating the rms image (see below).

Calculation of mean and rms maps

The box size and step size for calculating the rms image are estimated as follows (if not set by the rms_box parameter). Typical intersource seperation, \(s_1\), is \(2\sqrt{\text{bmpersrc_th}} \times B_{\text{major}}\). The size of brightest source, \(s_{\text{max}}\), is \(2 B_{\text{major}} \times \sqrt{[2\ln(\text{Max}_{\text{flux}}/\text{threshold})]}\). Lastly, the maximum dimension of the largest island, \(s_{\text{isl}}\), defined at 10–20 sigma above the clipped rms is also found. The box size is estimated as the larger of the quantities \(s_1\), \(s_{\text{max}}\), and \(s_{\text{isl}}\). The step size is then calculated as the minimum of a third of the box size and a tenth of the smallest image dimension. These prescriptions yield reasonable numbers for the images tested.

Either the calculated rms image or a constant rms is used for subsequent analysis based on whether the dispersion in the rms image is consistent with, or is higher than, the expected statistics. Hence if the dispersion of the rms image is higher than 1.1 times the (clipped) rms of the image times the inverse of \(\sqrt{2} \times \text{Boxsize}_{\text{pixels}}\) then the rms image is taken. Otherwise, the constant value of the clipped rms is used.

Gaussian fitting

The current procedure for calculating the number of Gaussians to be fit simultaneously to an island is as follows. First, the number of Gaussians is identified with the number of distinct peaks (higher than the pixel threshold) of emission inside the island (negative gradient in all 8 directions). These peaks are CLEANed from the subimage assuming the theoretical beam. If the (unclipped) rms of the residual subimage is greater than the (clipped) rms in the region, the maximum pixel in the residue is greater than the threshold for this former rms, and is located at least 0.5 beams (and \(\sqrt{5}\) pixels) away from all previous peaks, then this residual peak is identified as a new one.

Grouping of Gaussians into sources

Inside each island, groups of Gaussians are deemed to be a part of the same source if:

the difference between the minimum value along the line joining the centers of any pair of Gaussians and the peak value of the lower Gaussian is less than the product of the island threshold and the island rms, and

the centers are separated by a distance less than half the sum of their FWHMs along the line joining them.

Once the Gaussians that belong to a source are identified, fluxes for the grouped Gaussians are summed to obtain the total flux of the source. The uncertainty on the total flux is calculated by summing the uncertainties on the total fluxes of the individual Gaussians in quadrature. The source RA and Dec position is set to the source centroid determined from moment analysis (the position of the maximum of the source is also calculated). The total source size is also measured using moment analysis (see https://en.wikipedia.org/wiki/Image_moment for an overview of moment analysis).

Effect of neglecting color corrections

No color correction is performed when averaging channels. However, as is shown below, errors in the derived parameters are generally small unless the averaged bandwidth is large.

color correction errors — The correction in frequency in kHz as a function of the frequency resulting from averaging n channels each of bandwidth bw, for sources with various spectral indices \(-1.3 < \alpha < -0.3\).

Footnotes

PyBDSF Documentation

Why PyBDSF?

Capabilities of PyBDSF

Downloading and installing

What’s New

PyBDSF Basics

Starting PyBDSF

Quitting PyBDSF

Getting help

Logging

Commands

Tasks

Hierarchy of an astronomical image

Quick-start example

process_image: processing an image

General reduction parameters

Adaptive box options

Advanced options

Flagging options

Output options

Multichannel options

À trous wavelet decomposition module

PSF variation module

Spectral index module

Polarization module

Shapelet decomposition module

show_fit: visualizing the fit results

export_image: exporting internally derived images

write_catalog: writing source catalogs

Definition of output columns

Scripting PyBDSF

Alphabetical listing of all parameters

Simple image with point sources

Image with artifacts

Image with extended emission

Scripting example

Using SAMP interoperability

Determining whether an image is confused

Calculation of mean and rms maps

Gaussian fitting

Grouping of Gaussians into sources

Effect of neglecting color corrections

`process_image`: processing an image

`show_fit`: visualizing the fit results

`export_image`: exporting internally derived images

`write_catalog`: writing source catalogs