Math Tools

Coordinates Summary

dist_image(image[, pixscale, center, …])

Pixel distances

xy_to_rtheta(x, y)

Convert (x,y) to (r,theta)

rtheta_to_xy(r, theta)

Convert (r,theta) to (x,y)

xy_rot(x, y, ang)

Rotate (x,y) positions to new coords

Tel2Sci_info(channel, coords[, output])

Telescope coords converted to Science coords

det_to_V2V3(image, detid)

Same as det_to_sci

V2V3_to_det(image, detid)

Same as sci_to_det

plotAxes(ax[, position, label1, label2, …])

Compass arrows

Image Manipulation Summary

hist_indices(values[, bins, return_more])

Histogram indices

binned_statistic(x, values[, func, bins])

Binned statistic

frebin(image[, dimensions, scale, total])

Fractional rebin

fshift(image[, delx, dely, pad, cval])

Fractional image shift

fourier_imshift(image, xshift, yshift[, …])

Fourier shift image

shift_subtract(params, reference, target[, …])

Shift and subtract image

align_LSQ(reference, target[, mask, pad, …])

Find best shift value

scale_ref_image(im1, im2[, mask, …])

Reference image scaling

optimal_difference(im_sci, im_ref, scale[, …])

Optimize subtraction of ref PSF

pad_or_cut_to_size(array, new_shape[, fill_val])

Resize an array to a new shape by either padding with zeros or trimming off rows and/or columns.

fix_nans_with_med(im[, niter_max, verbose])

Iteratively fix NaNs with surrounding Real data

rotate_offset(data, angle[, cen, cval, …])

Rotate and offset an array.

Polynomial Fitting Summary

jl_poly_fit(x, yvals[, deg, QR, robust_fit, …])

Fast polynomial fitting

jl_poly(xvals, coeff[, dim_reorder])

Evaluate polynomial

Robust Summary

biweightMean(inputData[, axis, dtype, iterMax])

Biweight Mean

checkfit(inputData, inputFit, epsilon, delta)

Determine the quality of a fit and biweights.

linefit(inputX, inputY[, iterMax, Bisector, …])

Outlier resistance two-variable linear regression function.

mean(inputData[, Cut, axis, dtype, …])

Robust Mean

medabsdev(data[, axis, keepdims, nan])

Median Absolute Deviation

mode(inputData[, axis, dtype])

Robust estimator of the mode of a data set using the half-sample mode.

polyfit(inputX, inputY, order[, iterMax])

Outlier resistance two-variable polynomial function fitter.

std(inputData[, Zero, axis, dtype, …])

Robust Sigma


pynrc.maths.coords

pynrc.maths.coords.dist_image(image, pixscale=None, center=None, return_theta=False)[source]

Pixel distances

Returns radial distance in units of pixels, unless pixscale is specified. Use the center keyword to specify the position (in pixels) to measure from. If not set, then the center of the image is used.

return_theta will also return the angular position of each pixel relative to the specified center

Parameters
  • image (ndarray) – Input image to find pixel distances (and theta).

  • pixscale (int, None) – Pixel scale (such as arcsec/pixel or AU/pixel) that dictates the units of the output distances. If None, then values are in units of pixels.

  • center (tuple) – Location (x,y) in the array calculate distance. If set to None, then the default is the array center pixel.

  • return_theta (bool) – Also return the angular positions as a 2nd element.

pynrc.maths.coords.xy_to_rtheta(x, y)[source]

Convert (x,y) to (r,theta)

Input (x,y) coordinates and return polar cooridnates that use the WebbPSF convention (theta is CCW of +Y)

Input can either be a single value or numpy array.

Parameters
  • x (float or array) – X location values

  • y (float or array) – Y location values

pynrc.maths.coords.rtheta_to_xy(r, theta)[source]

Convert (r,theta) to (x,y)

Input polar cooridnates (WebbPSF convention) and return Carteesian coords in the imaging coordinate system (as opposed to RA/DEC)

Input can either be a single value or numpy array.

Parameters
  • r (float or array) – Radial offset from the center in pixels

  • theta (float or array) – Position angle for offset in degrees CCW (+Y).

pynrc.maths.coords.xy_rot(x, y, ang)[source]

Rotate (x,y) positions to new coords

Rotate (x,y) values by some angle. Positive ang values rotate counter-clockwise.

Parameters
  • x (float or array) – X location values

  • y (float or array) – Y location values

  • ang (float or array) – Rotation angle in degrees CCW

pynrc.maths.coords.Tel2Sci_info(channel, coords, output='Sci')[source]

Telescope coords converted to Science coords

Returns the detector name and position associated with input coordinates.

Parameters
  • channel (str) – ‘SW’ or ‘LW’

  • coords (tuple) – Telescope coordinates (V2,V3) in arcsec.

  • output (str) –

    Type of desired output coordinates.

    • Det: pixels, in raw detector read out axes orientation

    • Sci: pixels, in conventional DMS axes orientation

    • Idl: arcsecs relative to aperture reference location.

    • Tel: arcsecs V2,V3

pynrc.maths.coords.det_to_V2V3(image, detid)[source]

Same as det_to_sci

pynrc.maths.coords.V2V3_to_det(image, detid)[source]

Same as sci_to_det

pynrc.maths.coords.plotAxes(ax, position=(0.9, 0.1), label1='V2', label2='V3', dir1=[-1, 0], dir2=[0, 1], angle=0, alength=0.12, width=2, headwidth=8, color='w')[source]

Compass arrows

Show V2/V3 coordinate axis on a plot. By default, this function will plot the compass arrows in the lower right position in sky-right coordinates (ie., North/V3 up, and East/V2 to the left).

Parameters
  • ax (axis) – matplotlib axis to plot coordiante arrows.

  • position (tuple) – XY-location of joined arrows as a fraction (0.0-1.0).

  • label1 (str) – Label string for horizontal axis (ie., ‘E’ or ‘V2’).

  • label2 (str) – Label string for vertical axis (ie, ‘N’ or ‘V3’).

  • dir1 (array like) – XY-direction values to point “horizontal” arrow.

  • dir2 (array like) – XY-direction values to point “vertical” arrow.

  • angle (float) – Rotate coordinate axis by some angle. Positive values rotate counter-clockwise.

  • alength (float) – Length of arrow vectors as fraction of plot axis.

  • width (float) – Width of the arrow in points.

  • headwidth (float) – Width of the base of the arrow head in points.

  • color (color) – Self-explanatory.


pynrc.maths.image_manip

pynrc.maths.image_manip.hist_indices(values, bins=10, return_more=False)[source]

Histogram indices

This function bins an input of values and returns the indices for each bin. This is similar to the reverse indices functionality of the IDL histogram routine. It’s also much faster than doing a for loop and creating masks/indices at each iteration, because we utilize a sparse matrix constructor.

Returns a list of indices grouped together according to the bin. Only works for evenly spaced bins.

Parameters
  • values (ndarray) – Input numpy array. Should be a single dimension.

  • bins (int or ndarray) – If bins is an int, it defines the number of equal-width bins in the given range (10, by default). If bins is a sequence, it defines the bin edges, including the rightmost edge. In the latter case, the bins must encompass all values.

  • return_more (bool) – Option to also return the values organized by bin and the value of the centers (igroups, vgroups, center_vals).

Example

Find the standard deviation at each radius of an image

>>> rho = dist_image(image)
>>> binsize = 1
>>> bins = np.arange(rho.min(), rho.max() + binsize, binsize)
>>> igroups, vgroups, center_vals = hist_indices(rho, bins, True)
>>> # Get the standard deviation of each bin in image
>>> std = binned_statistic(igroups, image, func=np.std)
pynrc.maths.image_manip.binned_statistic(x, values, func=<function mean>, bins=10)[source]

Binned statistic

Compute a binned statistic for a set of data. Drop-in replacement for scipy.stats.binned_statistic.

Parameters
  • x (ndarray) – A sequence of values to be binned. Or a list of binned indices from hist_indices().

  • values (ndarray) – The values on which the statistic will be computed.

  • func (func) – The function to use for calculating the statistic.

  • bins (int or ndarray) – If bins is an int, it defines the number of equal-width bins in the given range (10, by default). If bins is a sequence, it defines the bin edges, including the rightmost edge. This doesn’t do anything if x is a list of indices.

Example

Find the standard deviation at each radius of an image

>>> rho = dist_image(image)
>>> binsize = 1
>>> radial_bins = np.arange(rho.min(), rho.max() + binsize, binsize)
>>> radial_stds = binned_statistic(rho, image, func=np.std, bins=radial_bins)
pynrc.maths.image_manip.frebin(image, dimensions=None, scale=None, total=True)[source]

Fractional rebin

Python port from the IDL frebin.pro Shrink or expand the size of a 1D or 2D array by an arbitary amount using bilinear interpolation. Conserves flux by ensuring that each input pixel is equally represented in the output array.

Parameters
  • image (ndarray) – Input image, 1-d or 2-d ndarray.

  • dimensions (tuple or None) – Desired size of output array (take priority over scale).

  • scale (tuple or None) – Factor to scale output array size. A scale of 2 will increase the number of pixels by 2 (ie., finer pixel scale).

  • total (bool) – Conserves the surface flux. If True, the output pixels will be the sum of pixels within the appropriate box of the input image. Otherwise, they will be the average.

Returns

ndarray – The binned ndarray

pynrc.maths.image_manip.fshift(image, delx=0, dely=0, pad=False, cval=0.0)[source]

Fractional image shift

Ported from IDL function fshift.pro. Routine to shift an image by non-integer values.

Parameters
  • image (ndarray) – 1D or 2D array to be shifted

  • delx (float) – shift in x (same direction as IDL SHIFT function)

  • dely (float) – shift in y

  • pad (bool) – Should we pad the array before shifting, then truncate? Otherwise, the image is wrapped.

  • cval (sequence or float, optional) – The values to set the padded values for each axis. Default is 0. ((before_1, after_1), … (before_N, after_N)) unique pad constants for each axis. ((before, after),) yields same before and after constants for each axis. (constant,) or int is a shortcut for before = after = constant for all axes.

Returns

ndarray – Shifted image

pynrc.maths.image_manip.fourier_imshift(image, xshift, yshift, pad=False, cval=0.0)[source]

Fourier shift image

Shift an image by use of Fourier shift theorem

Parameters
  • image (nd array) – N x K image

  • xshift (float) – Number of pixels to shift image in the x direction

  • yshift (float) – Number of pixels to shift image in the y direction

  • pad (bool) – Should we pad the array before shifting, then truncate? Otherwise, the image is wrapped.

  • cval (sequence or float, optional) – The values to set the padded values for each axis. Default is 0. ((before_1, after_1), … (before_N, after_N)) unique pad constants for each axis. ((before, after),) yields same before and after constants for each axis. (constant,) or int is a shortcut for before = after = constant for all axes.

Returns

ndarray – Shifted image

pynrc.maths.image_manip.shift_subtract(params, reference, target, mask=None, pad=False, shift_function=<function fshift>)[source]

Shift and subtract image

Use Fourier Shift theorem for subpixel shifts for input into least-square optimizer.

Parameters
  • params (tuple) – xshift, yshift, beta

  • reference (ndarray) – See align_fourierLSQ

  • target (ndarray) – See align_fourierLSQ

  • mask (ndarray, optional) – See align_fourierLSQ

  • pad (bool) – Should we pad the array before shifting, then truncate? Otherwise, the image is wrapped.

  • shift_function (func) – which function to use for sub-pixel shifting

Returns

ndarray – 1D array of target-reference residual after applying shift and intensity fraction.

pynrc.maths.image_manip.align_LSQ(reference, target, mask=None, pad=False, shift_function=<function fshift>)[source]

Find best shift value

LSQ optimization with option of shift alignment algorithm

Parameters
  • reference (ndarray) – N x K image to be aligned to

  • target (ndarray) – N x K image to align to reference

  • mask (ndarray, optional) – N x K image indicating pixels to ignore when performing the minimization. The masks acts as a weighting function in performing the fit.

  • pad (bool) – Should we pad the array before shifting, then truncate? Otherwise, the image is wrapped.

  • shift_function (func) – which function to use for sub-pixel shifting. Options are fourier_imshift or fshift. fshift tends to be 3-5 times faster for similar results.

Returns

list – (x, y, beta) values from LSQ optimization, where (x, y) are the misalignment of target from reference and beta is the fraction by which the target intensity must be reduced to match the intensity of the reference.

pynrc.maths.image_manip.scale_ref_image(im1, im2, mask=None, smooth_imgs=False, return_shift_values=False)[source]

Reference image scaling

Find value to scale a reference image by minimizing residuals. Assumes everything is already aligned if return_shift_values=False.

Or simply turn on return_shift_values to return (dx,dy,scl). Then fshift(im2,dx,dy) to shift the reference image.

Parameters
  • im1 (ndarray) – Science star observation.

  • im2 (ndarray) – Reference star observation.

  • mask (bool array or None) – Use this mask to exclude pixels for performing standard deviation. Boolean mask where True is included and False is excluded.

  • smooth_imgs (bool) – Smooth the images with nearest neighbors to remove bad pixels?

  • return_shift_values (bool) – Option to return x and y shift values

pynrc.maths.image_manip.optimal_difference(im_sci, im_ref, scale, binsize=1, center=None, mask_good=None, sub_mean=True, std_func=<function std>)[source]

Optimize subtraction of ref PSF

Scale factors from scale_ref_image work great for subtracting a reference PSF from a science image where there are plenty of photons, but perform poorly in the noise-limited regime. If we simply perform a difference by scaling the reference image, then we also amplify the noise. In the background, it’s better to simply subtract the unscaled reference pixels. This routine finds the radial cut-off of the dominant noise source.

Parameters
  • im_sci (ndarray) – Science star observation.

  • im_ref (ndarray) – Reference star observation.

  • scale (float) – Scale factor from scale_ref_image()

  • binsize (int) – Radial binsize (in pixels) to perform calculations

  • center (tuple or None) – Location (x,y) to calculate radial distances. Default is center of image.

  • mask_good (bool array) – Only perform operations on pixels where mask_good=True.

  • sub_mean (bool) – Subtract mean (median, actually) of pixels in each radial bin? Basically a background subtraction.

  • std_func (func) – What function do we want to use for calculating the standard deviation in each radial bin? After comparing the standard deviation between the two scaled differences in each radial bin, we only keep the better of the two.

pynrc.maths.image_manip.pad_or_cut_to_size(array, new_shape, fill_val=0.0)[source]

Resize an array to a new shape by either padding with zeros or trimming off rows and/or columns. The ouput shape can be of any arbitrary amount.

Parameters
  • array (ndarray) – A 1D or 2D array representing some image

  • padded_shape (tuple of 2 elements) – Desired size for the output array. For 2D case, if a single value, then will create a 2-element tuple of the same value.

  • fill_val (scalar, optional) – Value to pad borders. Default is 0.0

Returns

output (ndarray) – An array of size new_shape that preserves the central information of the input array.

pynrc.maths.image_manip.fix_nans_with_med(im, niter_max=5, verbose=False, **kwargs)[source]

Iteratively fix NaNs with surrounding Real data

pynrc.maths.image_manip.rotate_offset(data, angle, cen=None, cval=0.0, order=1, reshape=True, recenter=True, **kwargs)[source]

Rotate and offset an array.

Same as rotate in scipy.ndimage.interpolation except that it rotates around a center point given by cen keyword. The array is rotated in the plane defined by the two axes given by the axes parameter using spline interpolation of the requested order.

Parameters
  • data (ndarray) – The input array.

  • angle (float) – The rotation angle in degrees (rotates in CW direction).

  • cen (tuple) – Center location around which to rotate image. Values are expected to be (xcen, ycen).

  • recenter (bool) – Do we want to reposition so that cen is the image center?

Keyword Arguments
  • axes (tuple of 2 ints, optional) – The two axes that define the plane of rotation. Default is the first two axes.

  • reshape (bool, optional) – If reshape is true, the output shape is adapted so that the input array is contained completely in the output. Default is True.

  • order (int, optional) – The order of the spline interpolation, default is 1. The order has to be in the range 0-5.

  • mode (str, optional) – Points outside the boundaries of the input are filled according to the given mode (‘constant’, ‘nearest’, ‘reflect’, ‘mirror’ or ‘wrap’). Default is ‘constant’.

  • cval (scalar, optional) – Value used for points outside the boundaries of the input if mode='constant'. Default is 0.0

  • prefilter (bool, optional) – The parameter prefilter determines if the input is pre-filtered with spline_filter before interpolation (necessary for spline interpolation of order > 1). If False, it is assumed that the input is already filtered. Default is True.

Returns

rotate (ndarray or None) – The rotated data.


pynrc.maths.fast_poly

pynrc.maths.fast_poly.jl_poly_fit(x, yvals, deg=1, QR=True, robust_fit=False, niter=25)[source]

Fast polynomial fitting

Fit a polynomial to a function using linear least-squares. This function is particularly useful if you have a data cube and want to simultaneously fit a slope to all pixels in order to produce a slope image.

Gives the option of performing QR decomposition, which provides a considerable speed-up compared to simply using np.linalg.lstsq(). In addition to being fast, it has better numerical stability than linear regressions that involve matrix inversions (ie., dot(x.T,x)).

Returns the coefficients of the fit for each pixel.

Parameters
  • x (ndarray) – X-values of the data array (1D).

  • yvals (ndarray) – Y-values (1D, 2D, or 3D) where the first dimension must have equal length of x. For instance, if x is a time series of a data cube with size NZ, then the data cube must follow the Python convention (NZ,NY,NZ).

  • deg (int) – Degree of polynomial to fit to the data.

  • QR (bool) – Perform QR decomposition? Default=True.

  • robust_fit (bool) – Perform robust fitting, iteratively kicking out outliers until convergence.

  • niter (int) – Maximum number of iterations for robust fitting. If convergence is attained first, iterations will stop.

Example

Fit all pixels in a data cube to get slope image in terms of ADU/sec

>>> nz, ny, nx = cube.shape
>>> tvals = (np.arange(nz) + 1) * 10.737
>>> coeff = jl_poly_fit(tvals, cube, deg=1)
>>> bias = coeff[0]  # Bias image (y-intercept)
>>> slope = coeff[1] # Slope image (DN/sec)
pynrc.maths.fast_poly.jl_poly(xvals, coeff, dim_reorder=False)[source]

Evaluate polynomial

Replacement for np.polynomial.polynomial.polyval(wgood, coeff) to evaluate y-values given a set of xvals and coefficients. Uses matrix multiplication, which is much faster. Beware, the default output array shapes organization may differ from the polyval routine for 2D and 3D results.

Parameters
  • xvals (ndarray) – 1D array (time, for instance)

  • coeff (ndarray) – 1D, 2D, or 3D array of coefficients from a polynomial fit. The first dimension should have a number of elements equal to the polynomial degree + 1. Order such that lower degrees are first, and higher degrees are last.

  • dim_reorder (bool) – If true, then result to be ordered (nx,ny,nz), otherwise we use the Python preferred ordering (nz,ny,nx)

Returns

float array – An array of values where each xval has been evaluated at each set of supplied coefficients. The output shape has the first dimension equal to the number of xvals, and the final dimensions correspond to coeff’s latter dimensions. The result is flattened if there is either only one xval or one set of coeff (or both).


pynrc.maths.robust

Small collection of robust statistical estimators based on functions from Henry Freudenriech (Hughes STX) statistics library (called ROBLIB) that have been incorporated into the AstroIDL User’s Library.

pynrc.maths.robust.biweightMean(inputData, axis=None, dtype=None, iterMax=25)[source]

Biweight Mean

Calculate the mean of a data set using bisquare weighting.

Based on the biweight_mean routine from the AstroIDL User’s Library.

Changed in version 1.0.3: Added the ‘axis’ and ‘dtype’ keywords to make this function more compatible with np.mean()

pynrc.maths.robust.checkfit(inputData, inputFit, epsilon, delta, BisquareLimit=6.0)[source]

Determine the quality of a fit and biweights. Returns a tuple with elements:

  1. Robust standard deviation analog

  2. Fractional median absolute deviation of the residuals

  3. Number of input points given non-zero weight in the calculation

  4. Bisquare weights of the input points

  5. Residual values scaled by sigma

This function is based on the rob_checkfit routine from the AstroIDL User’s Library.

pynrc.maths.robust.linefit(inputX, inputY, iterMax=25, Bisector=False, BisquareLimit=6.0, CloseFactor=0.03)[source]

Outlier resistance two-variable linear regression function.

Based on the robust_linefit routine in the AstroIDL User’s Library.

pynrc.maths.robust.mean(inputData, Cut=3.0, axis=None, dtype=None, keepdims=False, return_std=False, return_mask=False)[source]

Robust Mean

Robust estimator of the mean of a data set. Based on the resistant_mean function from the AstroIDL User’s Library. NaN values are excluded.

This function trims away outliers using the median and the median absolute deviation. An approximation formula is used to correct for the truncation caused by trimming away outliers.

Parameters

inputData (ndarray) – The input data.

Keyword Arguments
  • Cut (float) – Sigma for rejection; default is 3.0.

  • axis (None or int or tuple of ints, optional) –

    Axis or axes along which the deviation is computed. The default is to compute the deviation of the flattened array.

    If this is a tuple of ints, a standard deviation is performed over multiple axes, instead of a single axis or all the axes as before. This is the equivalent of reshaping the input data and then taking the standard devation.

  • keepdims (bool, optional) – If this is set to True, the axes which are reduced are left in the result as dimensions with size one. With this option, the result will broadcast correctly against the original arr.

  • return_std (bool) – Also return the std dev calculated using only the “good” data?

  • return_mask (bool) – If set to True, then return only boolean array of good (1) and rejected (0) values.

pynrc.maths.robust.medabsdev(data, axis=None, keepdims=False, nan=True)[source]

Median Absolute Deviation

A “robust” version of standard deviation.

Parameters
  • data (ndarray) – The input data.

  • axis (None or int or tuple of ints, optional) – Axis or axes along which the deviation is computed. The default is to compute the deviation of the flattened array.

    If this is a tuple of ints, a standard deviation is performed over multiple axes, instead of a single axis or all the axes as before. This is the equivalent of reshaping the input data and then taking the standard devation.

  • keepdims (bool, optional) – If this is set to True, the axes which are reduced are left in the result as dimensions with size one. With this option, the result will broadcast correctly against the original arr.

  • nan (bool, optional) – Ignore NaNs? Default is True.

pynrc.maths.robust.mode(inputData, axis=None, dtype=None)[source]

Robust estimator of the mode of a data set using the half-sample mode.

pynrc.maths.robust.polyfit(inputX, inputY, order, iterMax=25)[source]

Outlier resistance two-variable polynomial function fitter.

Based on the robust_poly_fit routine in the AstroIDL User’s Library.

Unlike robust_poly_fit, two different polynomial fitters are used because np.polyfit does not support non-uniform weighting of the data. For the weighted fitting, the SciPy Orthogonal Distance Regression module (scipy.odr) is used.

pynrc.maths.robust.std(inputData, Zero=False, axis=None, dtype=None, keepdims=False, return_mask=False)[source]

Robust Sigma

Based on the robust_sigma function from the AstroIDL User’s Library.

Calculate a resistant estimate of the dispersion of a distribution.

Use the median absolute deviation as the initial estimate, then weight points using Tukey’s Biweight. See, for example, “Understanding Robust and Exploratory Data Analysis,” by Hoaglin, Mosteller and Tukey, John Wiley & Sons, 1983, or equation 9 in Beers et al. (1990, AJ, 100, 32).

Parameters

inputData (ndarray) – The input data.

Keyword Arguments
  • axis (None or int or tuple of ints, optional) –

    Axis or axes along which the deviation is computed. The default is to compute the deviation of the flattened array.

    If this is a tuple of ints, a standard deviation is performed over multiple axes, instead of a single axis or all the axes as before. This is the equivalent of reshaping the input data and then taking the standard devation.

  • keepdims (bool, optional) – If this is set to True, the axes which are reduced are left in the result as dimensions with size one. With this option, the result will broadcast correctly against the original arr.

  • return_mask (bool) – If set to True, then only return boolean array of good (1) and rejected (0) values.