Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Computes an iteratively sigma-clipped mean on a data set. Clipping is done about median, but mean is returned.
import numpy
def meanclip(indata, clipsig=3.0, maxiter=5, converge_num=0.02, verbose=0):
"""
Computes an iteratively sigma-clipped mean on a
data set. Clipping is done about median, but mean
is returned.
.. note:: MYMEANCLIP routine from ACS library.
:History:
* 21/10/1998 Written by RSH, RITSS
* 20/01/1999 Added SUBS, fixed misplaced paren on float call, improved doc. RSH
* 24/11/2009 Converted to Python. PLL.
Examples
--------
>>> mean, sigma = meanclip(indata)
Parameters
----------
indata: array_like
Input data.
clipsig: float
Number of sigma at which to clip.
maxiter: int
Ceiling on number of clipping iterations.
converge_num: float
If the proportion of rejected pixels is less than
this fraction, the iterations stop.
verbose: {0, 1}
Print messages to screen?
Returns
-------
mean: float
N-sigma clipped mean.
sigma: float
Standard deviation of remaining pixels.
"""
# Flatten array
skpix = indata.reshape( indata.size, )
ct = indata.size
iter = 0; c1 = 1.0 ; c2 = 0.0
while (c1 >= c2) and (iter < maxiter):
lastct = ct
medval = numpy.median(skpix)
sig = numpy.std(skpix)
wsm = numpy.where( abs(skpix-medval) < clipsig*sig )
ct = len(wsm[0])
if ct > 0:
skpix = skpix[wsm]
c1 = abs(ct - lastct)
c2 = converge_num * lastct
iter += 1
# End of while loop
mean = numpy.mean( skpix )
sigma = robust_sigma( skpix )
if verbose:
prf = 'MEANCLIP:'
print '%s %.1f-sigma clipped mean' % (prf, clipsig)
print '%s Mean computed in %i iterations' % (prf, iter)
print '%s Mean = %.6f, sigma = %.6f' % (prf, mean, sigma)
return mean, sigma
import numpy, pylab, scipy
from scipy import optimize
# Uses MEANCLIP from above
from meanclip import meanclip
def msky(inarray, do_plot=0, verbose=0, ptitle='', func=0):
"""
Find modal sky on an array.
First step is determination of median value and sigma.
Histogram of the data and fit parabola to the
logaritmic histogram. The coefficient of the parabola
are used to get mode and sigma of the sky on the
assumption that it is well fitted by a gaussian or
2nd-degree polynomial.
.. note:: MYSKY5 routine from ACS library.
:History:
* 11/25/2009 Created by PLL based on IDL routine by MS.
* 01/29/2010 PLL fixed broken where statement for histogram.
Parameters
----------
inarray: array_like
Input data.
do_plot: {0, 1}
Do plot?
verbose: {0, 1}
Print info to screen?
ptitle: string
Title of plot. Only used if `do_plot`=1.
func: {0, 1}
Function for fitting:
* 0: 2nd degree polynomial
* 1: Gaussian
Returns
-------
mmean: float
Mode of fitted function.
sigma: float
Sigma of fitted function.
"""
# Constants
func_name = 'MSKY'
nsig = 8.0
c1 = 2.5 # 2.8
c2 = 0.8 # 1.3
# Min/max of input array
arr_min = numpy.min(inarray)
arr_max = numpy.max(inarray)
# Get sigma
mmean, sigma = meanclip(inarray, clipsig=5.0, maxiter=10, verbose=verbose)
if sigma <= 0:
print ' '
print '%s: Weird distribution' % func_name
print '%-6s: %.6f' % ('MEAN', mmean)
print '%-6s: %.6f' % ('STDDEV', sigma)
print '%-6s: %.6f' % ('MIN', arr_min)
print '%-6s: %.6f' % ('MAX', arr_max)
return mmean, sigma
# Print info
if verbose:
print ' '
print '%s input array info' % func_name
print '%-6s: %.6f' % ('MIN', arr_min)
print '%-6s: %.6f' % ('MAX', arr_max)
# Flatten input array
arr_1d = inarray.reshape( inarray.size, )
# Define min and max for the histogram
x = nsig * sigma
mmean = numpy.median(arr_1d)
minhist = mmean - x
maxhist = mmean + x
ufi = inarray[ numpy.where((inarray > minhist) & (inarray < maxhist)) ]
# Calculate 25% and 75% percentile to get the interquartile range
# IRQ = pc75-pc25
# zenman, A. J. 1991.
sixd = numpy.argsort( ufi )
ndata = ufi.size
pc25 = ufi[ sixd[0.25*ndata] ]
pc75 = ufi[ sixd[0.75*ndata] ]
irq = pc75 - pc25
step = 2.0 * irq * ndata**(-1.0/3.0)
# Calculate number of bins to use
nbin = round(2*x / step - 1)
# Histogram
# http://www.scipy.org/Tentative_NumPy_Tutorial
yhist, hbin = numpy.histogram(arr_1d, range=(minhist,maxhist),bins=nbin)
xhist = 0.5 * ( hbin[1:] + hbin[:-1] )
# Define xmin and xmax for the 2-0rder fit
x1 = mmean - c1*sigma
x2 = mmean + c2*sigma
# Select the points beween x1 and x2 for the fit
w = numpy.where((xhist > x1) & (xhist < x2) & (yhist > 0))
count = len(w[0])
xwg = xhist[w]
nywg = yhist[w]
if count < 2:
print ' '
print '%s: Singular matrix' % func_name
print '%-6s: %.6f' % ('X[W]', xwg)
print '%-6s: %.6f' % ('NY[W]', nywg)
print '%-6s: %.6f' % ('MEDIAN', mmean)
return mmean, sigma
# Change to log scale
yhist = numpy.log10(yhist)
iyh = numpy.where( ~numpy.isinf(yhist) )
xhist = xhist[iyh]
yhist = yhist[iyh]
nywg = numpy.log10(nywg)
# Calculate the fit coefficients
ymax = nywg.max()
# Gaussian
# http://www.scipy.org/Cookbook/FittingData
if func == 1:
if verbose: print '%s: Fitting gaussian' % func_name
# Initial guess
ysum = numpy.sum(nywg)
mmean = numpy.sum(xwg * nywg) / ysum
sigma = numpy.sqrt(numpy.abs(numpy.sum((xwg-mmean)**2*nywg)/ysum))
params = (ymax, mmean, sigma)
# Error function to minimize
errorfunction = lambda p: scipy.ravel(gaussian(*p)(xhist) - yhist)
# Linear least square fitting
a_opt, a_success = optimize.leastsq(errorfunction, params)
ymax = a_opt[0]
mmean = a_opt[1]
sigma = a_opt[2]
# Fit to entire range
yall = ymax * numpy.exp(-(xhist-mmean)**2/(2.0*sigma**2))
# 2nd degree polynomial
else:
if verbose: print '%s: Fitting 2nd deg polynomial' % func_name
# Polynomial
a_opt = numpy.polyfit(xwg, nywg, 2)
mmean = -0.5*a_opt[1]/a_opt[0]
sigma = numpy.sqrt(-0.5 / ( a_opt[0]*numpy.log(10) ) )
# Fit to entire range
yall = a_opt[0]*xhist**2 + a_opt[1]*xhist + a_opt[2]
# Print results
if verbose:
print ' '
print '%s: Results' % func_name
print '%-6s: %.6f' % ('MODE', mmean)
print '%-6s: %.6f' % ('SIGMA', sigma)
print ' '
plot_y1 = 0
plot_y2 = ymax + 0.1
# Plot results
if do_plot:
pylab.clf()
# Data points
pylab.plot(xhist, yhist, 'ko')
# Mark fitted region
pylab.plot(xwg, nywg, 'bo', hold=1)
# Draw fitted function
pylab.plot(xhist, yall, 'r--', hold=1)
pylab.axvline(x=mmean, color='r')
# Plot xmin xmax for the fit
pylab.axvline(x=x1, ymin=0, ymax=0.1, color='b')
pylab.axvline(x=x2, ymin=0, ymax=0.1, color='b')
# Axis limits and labels
pylab.axis([minhist, maxhist, plot_y1, plot_y2])
pylab.xlabel('Pix value')
pylab.ylabel('Log num of pix')
pylab.title(ptitle)
#pylab.show()
# Pause for visual examination of plot
x_user = raw_input('ENTER ANY CHAR TO CONTINUE: ')
return mmean, sigma
def gaussian(height, center_x, width_x):
"""
Returns a gaussian function with the given parameters.
This is used for least square fitting optimization.
.. note:: This is used by `msky`.
Parameters
----------
height: float
Peak amplitude.
center_x: float
Peak location.
width_x: float
Sigma of gaussian curve.
Returns
-------
x: lambda function
Function used for optimization.
"""
return lambda x: height * numpy.exp(-(center_x-x)**2/(2.0*width_x**2))
import numpy
# Uses MEANCLIP from above
from meanclip import meanclip
def mytotal(inarray, axis, type='meanclip'):
"""
Collapse 2-D array in one dimension.
.. note:: MYTOTAL routine from ACS library.
:History:
* Obtained from M. Sirianni.
* Modified and converted to Python by P. L. Lim in 2009.
Examples
--------
>>> collapsed_array = mytotal(inarray, 1, type='median')
Parameters
----------
inarray: array_like
Input 2-D array.
axis: {1, 2}
Axis to collapse.
* 1: Return values along Y.
* 2: Return values along X.
type: {'median', 'meanclip', 'stdev'}
Algorithm to use.
Returns
-------
out_arr: array_like
1-D array collapsed along desired axis with desired
algorithm.
"""
func_name = 'MYTOTAL'
out_arr = 0.0
# Check inarray
if inarray.ndim != 2:
print '%s: Input array must be 2D' % func_name
return out_arr
# Check axis
if axis == 1:
n_out = inarray.shape[0]
elif axis == 2:
n_out = inarray.shape[1]
else:
print func_name, ': Axis not supported -', axis
return out_arr
# Initialize output array
out_arr = numpy.zeros(n_out)
out_rng = range(0, n_out)
# Check type
if type == 'meanclip':
for i in out_rng:
if axis == 1:
im_i = inarray[i,:]
else:
im_i = inarray[:,i]
mmean, msigma = meanclip(im_i, maxiter=10, converge_num=0.001)
out_arr[i] = mmean
elif type == 'stdev':
for i in out_rng:
if axis == 1:
im_i = inarray[i,:]
else:
im_i = inarray[:,i]
mmean, msigma = meanclip(im_i, maxiter=10, converge_num=0.001)
out_arr[i] = msigma
elif type == 'median':
for i in out_rng:
if axis == 1:
im_i = inarray[i,:]
else:
im_i = inarray[:,i]
out_arr[i] = numpy.median(im_i)
else:
print func_name, ': Type not supported -', type
out_arr = 0.0
return out_arr
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.