Math Tools¶
Coordinates Summary

Pixel distances 

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

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

Rotate (x,y) positions to new coords 

Telescope coords converted to Science coords 

Same as det_to_sci 

Same as sci_to_det 

Compass arrows 
Image Manipulation Summary

Histogram indices 

Binned statistic 

Fractional rebin 

Fractional image shift 

Fourier shift image 

Shift and subtract image 

Find best shift value 

Reference image scaling 

Optimize subtraction of ref PSF 

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

Iteratively fix NaNs with surrounding Real data 

Rotate and offset an array. 
Polynomial Fitting Summary

Fast polynomial fitting 

Evaluate polynomial 
Robust Summary

Biweight Mean 

Determine the quality of a fit and biweights. 

Outlier resistance twovariable linear regression function. 

Robust Mean 

Median Absolute Deviation 

Robust estimator of the mode of a data set using the halfsample mode. 

Outlier resistance twovariable polynomial function fitter. 

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 counterclockwise.
 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.
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 skyright coordinates (ie., North/V3 up, and East/V2 to the left).
 Parameters
ax (axis) – matplotlib axis to plot coordiante arrows.
position (tuple) – XYlocation of joined arrows as a fraction (0.01.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) – XYdirection values to point “horizontal” arrow.
dir2 (array like) – XYdirection values to point “vertical” arrow.
angle (float) – Rotate coordinate axis by some angle. Positive values rotate counterclockwise.
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) – Selfexplanatory.
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 equalwidth 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. Dropin 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 equalwidth 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, 1d or 2d 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 noninteger 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 leastsquare 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 subpixel shifting
 Returns
ndarray – 1D array of targetreference 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 subpixel shifting. Options are fourier_imshift or fshift. fshift tends to be 35 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 noiselimited 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 cutoff 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 2element 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 05.
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.0prefilter (bool, optional) – The parameter prefilter determines if the input is prefiltered 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 leastsquares. 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 speedup 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) – Xvalues of the data array (1D).
yvals (ndarray) – Yvalues (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 (yintercept) >>> 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 yvalues 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:
Robust standard deviation analog
Fractional median absolute deviation of the residuals
Number of input points given nonzero weight in the calculation
Bisquare weights of the input points
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 twovariable 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 halfsample mode.

pynrc.maths.robust.
polyfit
(inputX, inputY, order, iterMax=25)[source]¶ Outlier resistance twovariable 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 nonuniform 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.