sourcefinder.utils#

This module contain utilities for the source finding routines

Functions#

calculate_beamsize(semimajor, semiminor)

Calculate the beamsize based on the semi-major and minor axes.

calculate_correlation_lengths(semimajor, semiminor)

Calculate the Condon correlation lengths.

circular_mask(xdim, ydim, radius)

Returns a numpy array of shape (xdim, ydim). All points within

complement_gaussian_args(initial_params, fixed_params, ...)

Complements initial parameters for Gaussian fitting, with

flatten(nested_list)

Flatten a nested list

fudge_max_pix(semimajor, semiminor, theta)

Estimate peak spectral brightness correction at pixel of maximum spectral

generate_result_maps(data, sourcelist)

Return an image with Gaussian reconstructions of the sources and the

generate_subthresholds(min_value, max_value, ...)

Generate a series of num_thresholds logarithmically spaced values

get_error_radius(wcs, x_value, x_error, y_value, y_error)

Estimate an absolute angular error on the position (x_value,

interp_per_row(grid_row, y_initial, y_sought, interp_row)

Interpolate one row of the grid along the second dimension

is_valid_beam_tuple(→ bool)

make_subimages(a_data, a_mask, back_size_x, back_size_y)

Make subimages.

nearest_nonzero(some_arr, rms)

Replace values in some_arr based on the nearest non-zero

newton_raphson_root_finder(f, sigma0, min_sigma, max_sigma)

Solve the transcendental equation for sigma using Newton's

two_step_interp(grid, new_xdim, new_ydim)

Perform two-step interpolation.

Module Contents#

sourcefinder.utils.calculate_beamsize(semimajor, semiminor)[source]#

Calculate the beamsize based on the semi-major and minor axes.

Parameters:
semimajorfloat

The semi-major axis length in pixels.

semiminorfloat

The semi-minor axis length in pixels.

Returns:
float

The calculated beamsize.

sourcefinder.utils.calculate_correlation_lengths(semimajor, semiminor)[source]#

Calculate the Condon correlation lengths.

In order to derive the error bars for Gaussian fits from the Condon (1997, PASP 109, 116C) formulae, one needs a quantity called the correlation length. The Condon formulae assume a circular area with diameter theta_N (in pixels) for the correlation: i.e. all noise within that area is assumed completely correlated while outside that area all noise is assumed completely uncorrelated. This was later generalized by Hopkins et al. (2003, AJ 125, 465, paragraph 3) for correlation areas which are not circular, i.e. for anisotropic restoring beams.

Basically one has theta_N^2 = theta_B*theta_b.

Good estimates in general are:

  • theta_B = 2.0 * semimajor

  • theta_b = 2.0 * semiminor

but we have included this function to provide a convenient way of altering these dependencies.

Parameters:
semimajorfloat

The semi-major axis length in pixels.

semiminorfloat

The semi-minor axis length in pixels.

Returns:
tuple[float,float]

A tuple containing the correlation lengths (theta_B, theta_b), in pixels.

sourcefinder.utils.circular_mask(xdim, ydim, radius)[source]#

Returns a numpy array of shape (xdim, ydim). All points within radius from the centre are set to 0; outside that region, they are set to 1.

Parameters:
xdimint

The dimension of the array along the x-axis.

ydimint

The dimension of the array along the y-axis.

radiusfloat

The radius from the center within which points are set to 0.

Returns:
np.ndarray

A 2D ndarray with points within the radius set to 0 and outside set to 1.

sourcefinder.utils.complement_gaussian_args(initial_params, fixed_params, fit_params)[source]#

Complements initial parameters for Gaussian fitting, with fixed values. The end result should be a list of six elements, corresponding to ‘peak’, ‘xbar’, ‘ybar’, ‘semimajor’, ‘semiminor’ and ‘theta’, in that order

Parameters:
initial_paramsnp.ndarray

A Numpy float array containing the initial values for at least one, but a maximum of six Gaussian parameters: ‘peak’, ‘xbar’, ‘ybar’, ‘semimajor’, ‘semiminor’ and ‘theta’ , in that order.

fixed_paramsdict

A dictionary where the keys are zero, one or more of these parameter names: ‘peak’, ‘xbar’, ‘ybar’, ‘semimajor’, ‘semiminor’ and ‘theta’. The values are the fixed values for those parameters. It must complement initial_parms, such that the total number of parameters is six: len(initial_params) + len(fixed_params) == 6.

fit_paramstuple

A tuple of the six Gaussian parameters to be fitted. If a parameter is found in fixed_params, its value is used from there; otherwise, the value is taken from initial_params. Almost always this should be (‘peak’, ‘xbar’, ‘ybar’, ‘semimajor’, ‘semiminor’, ‘theta’). len(fit_params) == 6.

Returns:
gaussian_argslist[float]

A list of Gaussian fitting arguments, where each value corresponds to either a fixed value or an initial value from initial_params. len(gaussian_args) == 6.

Examples

>>> initial_parms = np.array([1.0, 4.0, 5.0, 6.0])
>>> fixed_parms = {'center_x': 2.5, 'center_y': 3.5}
>>> fit_parms = ('peak', 'center_x', 'center_y', 'semi-major axis',
...               'semi-minor axis', 'position angle')
>>> complement_gaussian_args(initial_parms, fixed_parms, fit_parms)
[1.0, 2.5, 3.5, 4.0, 5.0, 6.0]
sourcefinder.utils.flatten(nested_list)[source]#

Flatten a nested list

Nested lists are made in the deblending algorithm. They’re awful. This is a piece of code I grabbed from http://www.daniweb.com/code/snippet216879.html.

The output from this method is a generator, so make sure to turn it into a list, like this:

flattened = list(flatten(nested)).
Parameters:
nested_listlist

A list that may contain other lists or tuples.

Yields:
element

The next element in the flattened list.

Notes

The keyword “yield” is used; i.e. a generator object is returned.

sourcefinder.utils.fudge_max_pix(semimajor, semiminor, theta)[source]#

Estimate peak spectral brightness correction at pixel of maximum spectral brightness.

Previously, we adopted Rengelink’s correction for the underestimate of the peak of the Gaussian by the maximum pixel method: fudge_max_pix = 1.06. See the WENSS paper (1997A&AS..124..259R) or his thesis. The peak of the Gaussian is, of course, never at the exact center of the pixel, that’s why the maximum pixel method will underestimate it, when averaged over an ensemble. This effect is smaller when the clean beam is more densely sampled.

But, instead of just taking 1.06 one can make an estimate of the overall correction by assuming that the true peak is at a random position on the peak pixel and averaging over all possible corrections. This overall correction makes use of the beamshape, so strictly speaking only accurate for unresolved sources. After some investigation, it turns out that this method requires not only unresolved sources, but also a circular beam shape. With the general elliptical beam, the peak pixel may be located more than half a pixel away from the true peak. In that case the integral below will not suffice as a correction. Calculating an overall correction for an ensemble of sources will, in the case of an elliptical beam shape, become much more involved.

Parameters:
semimajorfloat

The semi-major axis length in pixels.

semiminorfloat

The semi-minor axis length in pixels.

thetafloat

The position angle of the major axis in radians.

Returns:
correctionfloat

The estimated peak spectral brightness correction.

sourcefinder.utils.generate_result_maps(data, sourcelist)[source]#

Return an image with Gaussian reconstructions of the sources and the corresponding residual image.

Given a data array (image) and list of sources, return two images, one showing the sources themselves and the other the residual after the sources have been removed from the input data.

Parameters:
datanp.ndarray

The input data array (image).

sourcelistlist

A list of sources to be removed from the input data.

Returns:
gaussian_mapnp.ndarray (2D)

Shows the Gaussian reconstructions of the sources.

residual_mapnp.ndarray (2D)

Shows the residuals from the subtractions of these reconstructions from the image data.

sourcefinder.utils.generate_subthresholds(min_value, max_value, num_thresholds)[source]#

Generate a series of num_thresholds logarithmically spaced values in the range (min_value, max_value) (both exclusive).

First, we calculate a logarithmically spaced sequence between exp(0.0) and (max_value - min_value + 1). That is, the total range is between 1 and one greater than the difference between max_value and min_value. We subtract 1 from this to get the range between 0 and (max_value - min_value). We add min to that to get the range between min and max.

This formula sets the subthreshold levels:

\[t_i = \exp\left( \frac{i \cdot \log(\text{max\_value} + 1 - \text{min\_value})} {\text{num\_thresholds} + 1} \right) + \text{min\_value} - 1\]

for \(i = 1, 2, \ldots, \text{num\_thresholds}\).

Parameters:
min_valuefloat

The minimum value of the range (not included in the output).

max_valuefloat

The maximum value of the range (not included in the output).

num_thresholdsint

The number of threshold values to generate.

Returns:
np.ndarray

An array of logarithmically spaced values between min_value and max_value (exclusive).

sourcefinder.utils.get_error_radius(wcs, x_value, x_error, y_value, y_error)[source]#

Estimate an absolute angular error on the position (x_value, y_value) with the given errors.

Parameters:
wcsobject

The WCS (World Coordinate System) object used for transforming pixel coordinates to sky coordinates.

x_valuefloat

Position along first pixel coordinate (row index of ndarray with image data).

x_errorfloat

The 1-sigma error in x-value.

y_valuefloat

Position along second pixel coordinate (column index of ndarray with image data).

y_errorfloat

The 1-sigma error in y-value.

Returns:
float

The estimated absolute angular error in arcseconds.

Notes

This is a pessimistic estimate, because we take sum of the error along the X and Y axes. Better might be to project them both back on to the major/minor axes of the elliptical fit, but this should do for now.

sourcefinder.utils.interp_per_row(grid_row, y_initial, y_sought, interp_row)[source]#

Interpolate one row of the grid along the second dimension (y-axis).

Parameters:
grid_row

1D array representing a single row of the grid.

y_initial

Original grid coordinates along the y-axis.

y_sought

Target coordinates for interpolation along the y-axis.

interp_row

Output array to store the interpolated row.

sourcefinder.utils.is_valid_beam_tuple(b) → bool[source]#
Return type:

bool

sourcefinder.utils.make_subimages(a_data, a_mask, back_size_x, back_size_y)[source]#

Make subimages.

Reshape the image data such that it is suitable for guvectorized kappa * sigma clipping. The idea is that we have designed a function that will perform kappa * sigma clipping on a single flattened subimage, i.e. on a 1D ndarray. If we decorate that function with Numba’s guvectorize decorator, we can use a 2D grid of subimages as input instead of a single flattened subimage. This means that the input to that guvectorized kappa * sigma clipper should be 3D. This function makes that 3D input, which is essentially a reshape using Numba with parallelization.

Parameters:
a_data: np.ndarray

The data of the masked array (without the mask).

a_mask: np.ndarray

The mask of the masked array (True means the value is masked).

back_size_x: int

The size of the subimage along the row indices.

back_size_y: int

The size of the subimages along the column indices.

Returns:
b: np.ndarray

3D array where each subimage of size (back_size_x * back_size_y) is flattened and padded with NaN.

c: np.ndarray

2D array where each element indicates the number of unmasked values in the corresponding subimage.

sourcefinder.utils.nearest_nonzero(some_arr, rms)[source]#

Replace values in some_arr based on the nearest non-zero values in rms.

Parameters:
some_arrnp.ndarray

A 2D array whose values will be replaced where rms == 0.

rmsnp.ndarray

A 2D array of the same shape as some_arr. Nearest non-zero neighbors in rms determine the replacement indices for some_arr.

Returns:
np.ndarray

A copy of some_arr with values replaced based on nearest non-zero neighbors in rms.

sourcefinder.utils.newton_raphson_root_finder(f, sigma0, min_sigma, max_sigma, tol=1e-08, max_iter=100, *args)[source]#

Solve the transcendental equation for sigma using Newton’s method with interval safeguards.

Parameters:
ffunction

The function to find the root of. It should take three parameters: sigma, sigma_meas, and D.

sigma0float

Initial guess for the value of sigma.

min_sigmafloat

Minimum bound for sigma.

max_sigmafloat

Maximum bound for sigma.

tolfloat, default: 1e-8

The tolerance for convergence.

max_iterint, default: 100

The maximum number of iterations.

*argstuple

Additional arguments for the function f.

Returns:
sigmafloat

The value of sigma that solves the equation, constrained within the bounds [min_sigma, max_sigma].

iint

The number of iterations performed. If convergence was reached, this is the iteration count. If max iterations were reached, this is equal to max_iter.

Raises:
ValueError

If the derivative of the function is near zero, causing a division by zero error.

Notes

The method employs Newton’s method for root-finding, with safeguards to ensure that sigma stays within the bounds [min_sigma, max_sigma]. The method terminates when the absolute change in sigma is smaller than the specified tolerance tol or when the maximum number of iterations is reached. If the derivative is too small (near zero), a ValueError is raised.

sourcefinder.utils.two_step_interp(grid, new_xdim, new_ydim)[source]#

Perform two-step interpolation.

It is done on a grid to upsample it to new dimensions. This function proivdes fast pieceswise bilinear interpolation in two steps:

  1. Interpolation across columns, i.e. per row of the input grid. Each row of the input grid is handled independently.

  2. Transpose the result.

  3. Again, interpolate across columns.

This method was inspired by a comment from “tiago” in this SO discussion: https://stackoverflow.com/q/14530556

Parameters:
gridnumpy.ndarray

The input 2D array to be upsampled.

new_xdimint

The desired number of rows in the upsampled grid.

new_ydimint

The desired number of columns in the upsampled grid.

Returns:
numpy.ndarray

The upsampled grid with dimensions (new_xdim, new_ydim).