from astropy.io import fits from scipy.optimize import curve_fit import glob from photutils import CircularAperture BLACK = 0 def OpenImageHist(filename): """Opens a FIT image and puts pixel ADU values into an array""" image = fits.getdata(filename) A = image.flatten() return A def OpenDirectory(Directory): filenames = glob.glob(Directory+"*.fit") return filenames def Histogram(Array, bins): """Plots histogram from flattened array of data, automatic range for gaussian distribution in array.""" h_test = plt.hist(Array, bins, (min(Array), max(Array))) #Histogram #1 over entire array plt.clf() test_y = h_test[0] test_x = h_test[1] PeakPosition = np.argmax(test_y) #Find peak of histogram Peak = test_x[PeakPosition] #Find position of peak along x axis BottomofRange = Peak-300 #Tailor new range around peak, very vague estimate, could probably use FWHM instead? Or a percentage of peak TopofRange = Peak+300 EstRange = (BottomofRange, TopofRange)#Estimate 1 of range Test2_h = plt.hist(Array, bins, EstRange) #Histogram #2 has slightly better range, still not great plt.clf() Test2_y = Test2_h[0] Test2_x = Test2_h[1] Counter = 0 #loop over the y values (Counts in histogram) to find the first non-zero value. for i in Test2_y: if i == 0: Counter+=1 if i >0: break Test2RangeLower = Test2_x[Counter] #The next test range is between all values that are non-zero, #this should be changed to a number dependent on the values in the histogram, #because they might not all be non-zero. Counter = 0 #loop over reversed list, to cut off any non-zero values. for i in reversed(Test2_y): if i == 0: Counter+=1 if i >0: break Test2RangeUpper = Test2_x[len(Test2_y) - Counter] Test2Range = (Test2RangeLower, Test2RangeUpper) Test3_h = plt.hist(Array, bins, Test2Range) plt.clf() Test3_x = Test3_h[1] Test3_y = Test3_h[0] Counter = 0 for i in Test3_y: if i<2000: #Instead of value here could have a percentage of peak count? Counter +=1 if i>=2000: break NewRangeLower = int(round(Test3_x[Counter])) Counter = 0 for i in reversed(Test3_y): if i<2000: Counter +=1 if i>=2000: break NewRangeHigher = int(round(Test3_x[len(Test3_y) - Counter])) FinalRange = (NewRangeLower, NewRangeHigher) Difference = NewRangeHigher-NewRangeLower bins = (Difference/2)-7 h = plt.hist(Array, bins, FinalRange) xedge = h[1] y = h[0] xcentre = (xedge + np.roll(xedge, -1))/2.0 x = xcentre[0:len(y)] plt.clf() return x, y, NewRangeLower, NewRangeHigher def GaussFunc(x,a,x0,sigma): return a*np.exp(-(x-x0)**2/(2*sigma**2)) def FitAGaussian(Array, bins): #Not Working """Fits a Gaussian Distribution over histogram""" h = Histogram(Array, bins) x = h[0] y = h[1] n = len(x) mean = sum(x)/n popt, pcov = curve_fit(GaussFunc, x, y, (1, mean, 1),absolute_sigma=True) return popt, pcov def FindExposureTime(path): head, tail = ntpath.split(path) B = tail #B = ntpath.basename(head) ExposureTime = float(B[:3]) return ExposureTime def Bins(array): binwidth = 10 bins = np.arange(0, max(array) + binwidth, binwidth) return bins def OpenImage(FileName): """Opens a FIT image and puts pixel ADU values into an array""" image = fits.getdata(FileName) return image def OpenDirectory(Directory): filenames = glob.glob(Directory+"*.fit") return filenames def count(pic): xsize, ysize = np.shape(pic) temp = pic clust = [] for x in range(xsize): for y in range(ysize): if temp[x,y] >= Threshold and temp[x,y] != WHITE: clust.append((fill(temp, xsize, ysize, x, y))) plt.imshow(temp, cmap='gray_r', origin='lower',interpolation='none', vmin=WHITE - 1,vmax=WHITE) return clust def fill(pic, xsize, ysize, x_init, y_init): thisClustList = [] queue = [(x_init, y_init)] while len(queue) != 0: top = queue.pop() x, y = top[0], top[1] pval = pic[x, y] if pval >= Threshold and pval != WHITE: thisClustList.append([x, y]) pic[x, y] = WHITE if x > 0: queue.append((x-1, y)) if x < (xsize - 1): queue.append((x+1, y)) if y > 0: queue.append((x, y-1)) if y < (ysize - 1): queue.append((x, y+1)) return thisClustList def blob_centre(blob, pic): blobxlist = [] blobylist = [] for pixel in blob: blobxlist.append(pixel[1]) blobylist.append(pixel[0]) nsum = 0 nxsum = 0 nysum = 0 for pixel in blob: nsum += pic[pixel[0]][pixel[1]] nxsum += (pic[pixel[0]][pixel[1]] * pixel[1]) nysum += (pic[pixel[0]][pixel[1]] * pixel[0]) xcen = nxsum/nsum ycen = nysum/nsum return xcen, ycen #path = "C:\\Users\\Blake\\Dropbox\\Msci_project\\SimulatedStars\\50StarsWithStotFluctuation.fits" path = "C:\\Users\\Blake\\Dropbox\\Msci_project\\Alpha\\alpha_persei_10.fit" # picHist = OpenImageHist(path) picFill = OpenImage(path) pic = OpenImage(path) WHITE = max(picHist) + 1 bins = Bins(picHist) popt, pcov = FitAGaussian(picHist, bins) mean = popt[1] std = abs(popt[2]) Threshold = mean+std*6 #Threshold = picHist.mean() + picHist.std() plt.figure(1) plt.imshow(pic, cmap='gray', origin='lower',interpolation='none', vmin=np.median(pic),vmax=np.median(pic)*2.2) plt.colorbar() plt.figure(2) clist = count(picFill) print len(clist) bloblist = [] # i = 0 # for blob in clist: # print i # print len(blob) # i+=1 blob = clist[179] for pixel in blob: bloblist.append(pic[pixel[0]][pixel[1]]) plt.figure(3) print len(bloblist) plt.hist(bloblist, bins = np.arange(0, max(bloblist) + 10,10)) plt.axis([200, 4000, 0, 120]) # maxb = max(bloblist) # bloblist.index(maxb) # centreCo = blob[] # apertures = CircularAperture(centreCo, r=50.) # apertures.plot(color='blue', lw=1.5, alpha=0.5) plt.show() <\verbatim> -- Public.LunaBorella - 27 Oct 2017