Difference: StewartBoogertPhotometryStarFinder (1 vs. 3)

Revision 323 Nov 2017 - LunaMarieBorella1

Line: 1 to 1
 
META TOPICPARENT name="StewartBoogertPhotometry2017"
Added:
>
>
from astropy.io import fits

from scipy.optimize import curve_fit import glob from photutils import CircularAperture

 BLACK = 0
Changed:
<
<
WHITE=5000
>
>

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"""
Line: 28 to 147
  return clust
Changed:
<
<
def fill(pic, xsize, ysize, x, y):

if pic[x, y] < Threshold: return

>
>
def fill(pic, xsize, ysize, x_init, y_init):
  thisClustList = []
Added:
>
>
queue = [(x_init, y_init)]
 
Deleted:
<
<
queue = [(x, y)]
  while len(queue) = 0:
Changed:
<
<
next = queue.pop() if pic[next[0], next[1]] >= Threshold and pic[next[0], next[1]] = WHITE: pic[next[0], next[1]] = WHITE thisClustList.append((next[0], next[1]))
>
>
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))
Deleted:
<
<
return thisClustList
 
Deleted:
<
<
pic = OpenImage("C:\\Users\\Blake\\Dropbox\\Msci_project\\Alpha\\alpha_persei_10.fit")
 
Changed:
<
<
Threshold = pic.mean()+(np.std(pic))
>
>
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)

Changed:
<
<
stars = count(pic) print stars print len(stars) count(pic) plt.colorbar()
>
>
clist = count(picFill)

print len(clist)

bloblist = [] # i = 0 # for blob in clist: # print i # print len(blob) # i+=1

 
Added:
>
>
blob = clist[179]
 
Added:
>
>
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()

Revision 209 Nov 2017 - LunaMarieBorella1

Line: 1 to 1
 
META TOPICPARENT name="StewartBoogertPhotometry2017"
Deleted:
<
<
from astropy.io import fits import numpy as np import matplotlib.pyplot as plt from scipy.optimize import curve_fit import glob import ntpath import math as m from scipy import stats import operator
 BLACK = 0
Changed:
<
<
WHITE=2**16
>
>
WHITE=5000
 def OpenImage(FileName): """Opens a FIT image and puts pixel ADU values into an array""" image = fits.getdata(FileName)
Line: 26 to 17
 def count(pic): xsize, ysize = np.shape(pic) temp = pic
Changed:
<
<
result = 0
>
>
clust = []
  for x in range(xsize): for y in range(ysize):
Changed:
<
<
if temp[x,y] >= Threshold: result += 1 fill(temp, xsize, ysize, x, y) plt.imshow(temp, cmap=plt.cm.binary) return result
>
>
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, y):
Changed:
<
<
if pic[x, y] = Threshold: pic[x,y] = WHITE
>
>
if pic[x, y] < Threshold:
  return
Deleted:
<
<
pic[x,y] = BLACK if x>0: fill(pic, xsize, ysize, x-1, y) if x < (xsize-1): fill(pic, xsize, ysize, x+1, y) if y > 0: fill(pic, xsize, ysize, x, y-1) if y < (ysize-1): fill(pic, xsize, ysize, x, y+1)
 
Added:
>
>
thisClustList = []

queue = [(x, y)] while len(queue) = 0: next = queue.pop() if pic[next[0], next[1]] >= Threshold and pic[next[0], next[1]] = WHITE: pic[next[0], next[1]] = WHITE thisClustList.append((next[0], next[1])) 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

 
Added:
>
>
pic = OpenImage("C:\\Users\\Blake\\Dropbox\\Msci_project\\Alpha\\alpha_persei_10.fit")
 
Added:
>
>
Threshold = pic.mean()+(np.std(pic))
 
Deleted:
<
<
pic = OpenImage("C:\\Users\\Blake\\Dropbox\\Msci_project\\Alpha\\alpha_persei.fit") Threshold = pic.mean()+(np.std(pic)/20)
 plt.figure(1)
Changed:
<
<
plt.imshow(pic)
>
>
plt.imshow(pic, cmap='gray', origin='lower',interpolation='none', vmin=np.median(pic),vmax=np.median(pic)*2.2) plt.colorbar()
 plt.figure(2)
Changed:
<
<
print count(pic)
>
>
stars = count(pic) print stars print len(stars) count(pic) plt.colorbar()
 
Deleted:
<
<
plt.show()
 
Added:
>
>
plt.show()
 

<\verbatim>

Revision 127 Oct 2017 - LunaMarieBorella1

Line: 1 to 1
Added:
>
>
META TOPICPARENT name="StewartBoogertPhotometry2017"
from astropy.io import fits
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import glob
import ntpath
import math as m
from scipy import stats
import operator

BLACK = 0
WHITE=2**16
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
   result = 0
   for x in range(xsize):
      for y in range(ysize):
         if temp[x,y] >= Threshold:
            result += 1
            fill(temp, xsize, ysize, x, y)
   plt.imshow(temp, cmap=plt.cm.binary)         
   return result

def fill(pic, xsize, ysize, x, y):
   if pic[x, y] != Threshold:
      pic[x,y] = WHITE
      return 
   pic[x,y] = BLACK
   if x>0: fill(pic, xsize, ysize, x-1, y)
   if x < (xsize-1): fill(pic, xsize, ysize, x+1, y)
   if y > 0: fill(pic, xsize, ysize, x, y-1)
   if y < (ysize-1): fill(pic, xsize, ysize, x, y+1)




pic = OpenImage("C:\\Users\\Blake\\Dropbox\\Msci_project\\Alpha\\alpha_persei.fit")
Threshold = pic.mean()+(np.std(pic)/20)
plt.figure(1)
plt.imshow(pic)
plt.figure(2)
print count(pic)

plt.show()



<\verbatim>



-- Public.LunaBorella - 27 Oct 2017
 
This site is powered by the TWiki collaboration platform Powered by PerlCopyright © 2008-2021 by the contributing authors. All material on this collaboration platform is the property of the contributing authors.
Ideas, requests, problems regarding RHUL Physics Department TWiki? Send feedback