notebook.community

Edit and run 





In [2]:



    


import numpy as np
from astropy.io import fits as pf
import matplotlib.pyplot as plt



    











In [97]:



    


#Import SpIES / SHELA data
data = '/Users/johntimlin/Complete_Clustering_Analysis/Clustering/Data_sets/HZLZ_combined_all_hzclassifiers_wphotoz_zspecflg.fits'
obs = pf.open(data)[1].data
Z = obs.zbest
imag = -2.5/np.log(10) * (np.arcsinh((obs.iflux/1e9)/(2*1.8e-10))+np.log(1.8e-10)) 

dx = ((Z >= 2.9) & (obs.dec>=-1.2) & (obs.dec<=1.2)) & ((obs.Au <0.3) | (obs.zspec>2.9))
gdx = ((Z >= 2.9)&(imag >= 20.2) & (obs.dec>=-1.2) & (obs.dec<=1.2))

rand = '/Users/johntimlin/Complete_Clustering_Analysis/Clustering/Data_sets/Short_randoms_spsh.fits'
rnd = pf.open(rand)[1].data
rdx = ((rnd.DEC>=-1.2) & (rnd.DEC<=1.2)) #& ((rnd.RA>344) | (rnd.RA<60))

RAr = rnd.RA[rdx]
RAd = obs.ra[dx]


rRA = []
for i in range(len(RAr)):
    if RAr[i]>300:
        rRA.append(RAr[i]-360)
    else:
        rRA.append(RAr[i])

        
dRA = []
for j in range(len(RAd)):
    if RAd[j]>300:
        dRA.append(RAd[j]-360)
    else:
        dRA.append(RAd[j])
        

print len(dRA)



    


















    






1592

















In [ ]:



    






    











In [ ]:



    






    











In [98]:



    


#Generate the redshift histogram for the faint objects

rbins = np.arange(-30,37,5)
dbins = np.arange(-1.2,1.2,0.25)

tdnum, rabins, decbins = np.histogram2d(dRA,obs.dec[dx],bins = [rbins,dbins])
rdnum, rabins, decbins = np.histogram2d(rRA,rnd.DEC[rdx],bins = [rbins,dbins])

ra_bin_midpoints = rabins[:-1] + np.diff(rabins)/2
dec_bin_midpoints = decbins[:-1] + np.diff(decbins)/2

rdnum = rdnum*(float(len(dRA))/len(rRA)) / (0.25*5)
tdnum = tdnum / (0.25*5)


denratio = tdnum/rdnum

print denratio

plt.figure(1)
plt.hist(denratio.flatten(),bins = 20)
plt.show()



    


















    






[[ 0.59460947  1.09986327  0.31933791  0.48678416  0.49818631  0.91402778
   1.02061324  1.71678911  2.11416701]
 [ 0.85734389  0.70354556  0.55264141  0.7359776   0.92759078  0.55264141
   0.6739632   0.95342553  0.73243286]
 [ 0.          0.48629493  0.97471235  0.40886271  0.61709625  0.7164729
   0.69209863  0.48416606  0.99311917]
 [ 0.37363974  1.35609844  1.5517487   1.18696664  1.14994199  1.01032913
   0.6731428   1.05645241  1.45934313]
 [ 2.14752159  1.1852363   1.3611983   1.21982384  0.7825285   1.00543056
   0.74590495  0.83011877  1.16607113]
 [ 1.75551368  1.21833079  0.74912006  0.97414     0.9770085   0.96591582
   0.54015806  1.11926489  1.33003303]
 [ 0.36025199  0.76513898  1.03908267  0.95893666  0.93106618  0.69514369
   1.2139514   1.205056    1.82805557]
 [ 1.70149787  1.45522844  1.30114543  0.9569611   0.78068726  1.28154533
   1.08491353  1.10208391  1.30421417]
 [ 1.09864929  1.05771536  0.9172767   1.0336202   1.15272294  1.70885911
   1.35535983  1.00258232  1.0564039 ]
 [ 0.85939373  0.82440694  1.18720332  1.02523626  0.78484231  0.77026601
   0.89449051  1.08994501  1.28228825]
 [ 0.75469664  1.16418276  0.60070987  1.32493687  1.30762282  0.84640838
   0.97873772  0.53012516  1.03121083]
 [ 1.12854451  1.1863841   0.93387398  0.92384431  0.31671639  1.20651667
   1.36299566  0.47206191  1.16828199]
 [ 0.89191421  1.16486396  1.63343796  0.79757713  1.25528666  1.15138954
   1.23900346  0.76493853  1.14270337]]










    
























In [96]:



    


extent = [ra_bin_midpoints[0], ra_bin_midpoints[-1], dec_bin_midpoints[0], dec_bin_midpoints[-1]]

plt.figure(1,figsize = (50,5))
#plt.hist2d(dRA, obs.dec[dx], bins=(rabins, decbins))
plt.imshow(tdnum,extent = extent)
#plt.contourf(ra_bin_midpoints,dec_bin_midpoints,tdnum.T)
#plt.colorbar()#ticks=np.arange(0,10,0.25))
#plt.savefig('density.png')

plt.figure(2)
#plt.hist2d(rRA, rnd.DEC[rdx], bins=(rabins, decbins))
plt.imshow(rdnum,extent = extent)
#plt.contourf(ra_bin_midpoints,dec_bin_midpoints,rdnum.T)
#plt.colorbar()#ticks=np.arange(0,10,0.25))
#plt.savefig('density.png')

plt.show()



    


















    

















    
























In [99]:



    


# Examine the densities of the MW plane and the non-MW plane



    











In [111]:



    


#data
nmw = ((obs.ra > 344) | (obs.ra<60)) & ((obs.zbest >= 2.9) & (obs.dec>=-1.2) & (obs.dec<=1.2))
mw = ((obs.ra < 344) & (obs.ra>60)) & ((obs.zbest >= 2.9) & (obs.dec>=-1.2) & (obs.dec<=1.2))
allx = ((obs.zbest >= 2.9) & (obs.dec>=-1.2) & (obs.dec<=1.2))
# randoms
rnmw = ((rnd.ra > 344) | (rnd.ra<60)) & ((rnd.zbest >= 2.9) & (rnd.dec>=-1.2) & (rnd.dec<=1.2))
rmw = ((rnd.ra < 344) & (rnd.ra>60)) & ((rnd.zbest >= 2.9) & (rnd.dec>=-1.2) & (rnd.dec<=1.2))
rallx = ((rnd.zbest >= 2.9) & (rnd.dec>=-1.2) & (rnd.dec<=1.2))



    


















    






1802
1385
417

















In [ ]:



    






    











In [ ]:



    






    











In [ ]:



    






    











In [ ]:



    






    











In [ ]:



    






    











In [ ]:



    






    











In [ ]:



    


## Do with the KDE instead



    











In [4]:



    


import numpy as np
from astropy.io import fits as pf
import matplotlib.pyplot as plt
from sklearn.neighbors import KernelDensity as kde



    











In [121]:



    


#Import SpIES / SHELA data
data = '/Users/johntimlin/Complete_Clustering_Analysis/Clustering/Data_sets/HZLZ_combined_all_hzclassifiers_wphotoz_zspecflg.fits'
obs = pf.open(data)[1].data
Z = obs.zbest
imag = -2.5/np.log(10) * (np.arcsinh((obs.iflux/1e9)/(2*1.8e-10))+np.log(1.8e-10)) 

dx = ((Z >= 2.9) & (obs.dec>=-1.2) & (obs.dec<=1.2))
gdx = ((Z >= 2.9)&(imag >= 20.2) & (obs.dec>=-1.2) & (obs.dec<=1.2))

rand = '/Users/johntimlin/Complete_Clustering_Analysis/Clustering/Data_sets/Short_randoms_spsh.fits'
rnd = pf.open(rand)[1].data
rdx = ((rnd.DEC>=-1.2) & (rnd.DEC<=1.2))

RAr = rnd.RA[rdx]
RAd = obs.ra[dx]


rRA = []
for i in range(len(RAr)):
    if RAr[i]>300:
        rRA.append(RAr[i]-360)
    else:
        rRA.append(RAr[i])

        
dRA = []
for j in range(len(RAd)):
    if RAd[j]>300:
        dRA.append(RAd[j]-360)
    else:
        dRA.append(RAd[j])



    











In [122]:



    


#Generate the redshift histogram for the faint objects

rbins = np.arange(-30,37,2)
dbins = np.arange(-1.2,1.2,0.1)


dranum, bins = np.histogram(dRA,rbins)
ddecnum, bins = np.histogram(obs.dec[dx],dbins)

rranum, bins = np.histogram(rRA,rbins)
rdecnum, bins = np.histogram(rnd.DEC[rdx],dbins)

rranum = rranum * float(len(obs.ra[dx]))/len(rnd.RA[rdx])
rdecnum = rdecnum * float(len(obs.ra[dx]))/len(rnd.RA[rdx])

plt.figure(1)
plt.plot(rbins[:-1],dranum,linestyle = 'steps-mid', label = 'data')
plt.plot(rbins[:-1],rranum,linestyle = 'steps-mid', label = 'rand')
plt.legend()

plt.figure(2)
plt.plot(dbins[:-1],ddecnum,linestyle = 'steps-mid', label = 'data')
plt.plot(dbins[:-1],rdecnum,linestyle = 'steps-mid', label = 'rand')

plt.legend()
plt.show()



    


















    

















    
























In [139]:



    


#Set up a KDE for dNdz
#tmpz = Z[gdx][:, np.newaxis] #change the array from row shape (1) to column shape (1,)
tmpz = np.hstack([np.asarray(dRA)[:, np.newaxis],obs.dec[dx][:, np.newaxis]])

est = kde(bandwidth=0.25,kernel='gaussian') #Set up the Kernel
histkde = est.fit(tmpz) #fit the kernel to the data and find the density of the grid

val = histkde.sample(n_samples = 100*len(dRA))

new_ra = val[:,0]
new_dec = val[:,1]

print new_ra[0], new_dec[0]



    


















    






24.4794236946 0.22450500458

















In [ ]:



    






    











In [140]:



    


plt.scatter(new_ra,new_dec,s=1)
plt.scatter(dRA,obs.dec[dx],s=1,color = 'c')

plt.show()



    


















    
























In [141]:



    


rbins = np.arange(-30,37,2)
dbins = np.arange(-1.2,1.2,0.11)


dranum, bins = np.histogram(dRA,rbins)
ddecnum, bins = np.histogram(obs.dec[dx],dbins)

rranum, bins = np.histogram(new_ra,rbins)
rdecnum, bins = np.histogram(new_dec,dbins)

rranum = rranum * float(len(obs.ra[dx]))/len(rnd.RA[rdx])
rdecnum = rdecnum * float(len(obs.ra[dx]))/len(rnd.RA[rdx])

plt.figure(1)
plt.plot(rbins[:-1],dranum,linestyle = 'steps-mid', label = 'data')
plt.plot(rbins[:-1],rranum,linestyle = 'steps-mid', label = 'rand')
plt.legend()

plt.figure(2)
plt.plot(dbins[:-1],ddecnum,linestyle = 'steps-mid', label = 'data')
plt.plot(dbins[:-1],rdecnum,linestyle = 'steps-mid', label = 'rand')

#plt.legend()
plt.show()



    


















    

















    
























In [142]:



    


rabins = np.arange(-30,37,2)
decbins = np.arange(-1.2,1.2,0.1)



plt.subplot(121)#, aspect='equal')
plt.hist2d(dRA, obs.dec[dx], bins=(rabins, decbins))
plt.subplot(122)#, aspect='equal')
plt.hist2d(new_ra, new_dec, bins=(rabins, decbins))
plt.show()



    


















    
























In [ ]:



    






    











In [143]:



    


tbhdu=pf.BinTableHDU.from_columns([pf.Column(name='ra',format='D',array=np.asarray(new_ra)), pf.Column(name='dec',format='D',array=np.asarray(new_dec))])
		

prihdr=pf.Header()
prihdr['COMMENT']="position matched sample"
prihdu=pf.PrimaryHDU(header=prihdr)
	
hdulist=pf.HDUList([prihdu,tbhdu])
hdulist.writeto('Position_match_random2.fits')



    











In [ ]:

Content source: JDTimlin/QSO_Clustering

Similar notebooks:

notebook.community | gallery | about