notebook.community

Edit and run

Bootcalib code for r camera





In [2]:



    


%matplotlib inline



    











In [3]:



    


# imports
try:
    import seaborn as sns; sns.set(context="notebook",font_scale=2)
except:
    pass

from desispec import bootcalib as desiboot
from desiutil import funcfits as dufits
from astropy.io import fits
from astropy.stats import sigma_clip

import numpy as np
from astropy.modeling import models, fitting



    


















    






/Users/xavier/anaconda/lib/python2.7/site-packages/matplotlib/__init__.py:872: UserWarning: axes.color_cycle is deprecated and replaced with axes.prop_cycle; please use the latter.
  warnings.warn(self.msg_depr % (key, alt_key))

Flat tracing

Read flat





In [4]:



    


fiberflat = '/Users/xavier/DESI/Wavelengths/pix-r0-00000001.fits'
flat_hdu = fits.open(fiberflat)
header = flat_hdu[0].header
flat = flat_hdu[0].data
ny = flat.shape[0]

Find fibers





In [5]:



    


xpk, ypos, cut = desiboot.find_fiber_peaks(flat)



    


















    






INFO:bootcalib.py:684:find_fiber_peaks: starting










    






INFO:DESI:starting










    






INFO:bootcalib.py:721:find_fiber_peaks: Found 500 fibers










    






INFO:DESI:Found 500 fibers










    






INFO:bootcalib.py:729:find_fiber_peaks: Found 20 bundles










    






INFO:DESI:Found 20 bundles

















In [6]:



    


tst = True
if tst:
    xpk = xpk[0:5]



    











In [7]:



    


reload(desiboot)
desiboot.qa_fiber_peaks(xpk, cut)

Trace the fiber flat spectra





In [8]:



    


# Crude first
xset, xerr = desiboot.trace_crude_init(flat,xpk,ypos)
# Polynomial fits
xfit, fdicts = desiboot.fit_traces(xset,xerr)



    











In [9]:



    


# QA
desiboot.qa_fiber_Dx(xfit, fdicts)

Model the PSF with Gaussian





In [10]:



    


gauss = desiboot.fiber_gauss(flat,xfit,xerr)



    


















    






WARNING:bootcalib.py:587:fiber_gauss: fiber_gauss uses astropy.modeling.  Consider an alternative










    






WARNING:DESI:fiber_gauss uses astropy.modeling.  Consider an alternative










    






INFO:bootcalib.py:606:fiber_gauss: Working on fiber 0 of 5










    






INFO:DESI:Working on fiber 0 of 5
WARNING: AstropyDeprecationWarning: The "sig" keyword is now deprecated, use the "sigma" keyword instead. [astropy.stats.sigma_clipping]
WARNING:astropy:AstropyDeprecationWarning: The "sig" keyword is now deprecated, use the "sigma" keyword instead.

















In [11]:



    


reload(desiboot)
desiboot.qa_fiber_gauss(gauss)

Arc

Wavelength info





In [12]:



    


import desimodel.io



    











In [13]:



    


desi_psf = desimodel.io.load_psf('r')



    











In [14]:



    


wave4 = desi_psf.wavelength(4,np.arange(desi_psf.npix_y))



    











In [15]:



    


wave4[1692]



    


















    
Out[15]:








6506.600473176205

















In [16]:



    


wave4[2451]



    


















    
Out[16]:








6907.4548001103185

















In [17]:



    


wave4[2493]



    


















    
Out[17]:








6929.2715774131439

















In [18]:



    


wave4[2692]



    


















    
Out[18]:








7032.0935157913

















In [19]:



    


wave4[2969]



    


















    
Out[19]:








7173.6943327592153

















In [22]:



    


np.min(wave4), np.max(wave4)



    


















    
Out[22]:








(5571.4062126915705, 7747.3883981124163)

















In [24]:



    


np.abs(wave4[0]-wave4[-1])/desi_psf.npix_y



    


















    
Out[24]:








0.52712746739846073

Read arc





In [25]:



    


arcfile = '/Users/xavier/DESI/Wavelengths/pix-r0-00000000.fits'
arc_hdu = fits.open(arcfile)
arc = arc_hdu[0].data

Extract arc





In [26]:



    


all_spec = desiboot.extract_sngfibers_gaussianpsf(arc,xfit,gauss,verbose=False)

Linelist





In [27]:



    


reload(desiboot)
camera = header['CAMERA']
llist = desiboot.load_arcline_list(camera)
dlamb, wmark, gd_lines, line_guess = desiboot.load_gdarc_lines(camera)



    


















    






ERROR:bootcalib.py:453:load_arcline_list: Not ready for this camera










    






ERROR:DESI:Not ready for this camera










    






INFO:bootcalib.py:469:load_arcline_list: Rejecting select HgI lines










    






INFO:DESI:Rejecting select HgI lines










    






INFO:bootcalib.py:469:load_arcline_list: Rejecting select NeI lines










    






INFO:DESI:Rejecting select NeI lines

















In [28]:



    


gd_lines



    


















    
Out[28]:








array([ 5769.598 ,  5852.4878,  5944.834 ,  6143.062 ,  6402.246 ,
        6506.528 ,  6678.2766,  6717.043 ,  6929.4672,  7032.4128,
        7173.938 ,  7245.1665,  7438.898 ])

Solve





In [29]:



    


ii=4

Find Lines





In [30]:



    


spec = all_spec[:,ii]
# Find Lines
pixpk = desiboot.find_arc_lines(spec)



    











In [31]:



    


pixpk



    


















    
Out[31]:








array([  340.20175748,   349.39823329,   386.90565798,   497.12220262,
         549.7341551 ,   662.71309491,   717.95235434,   816.36200462,
         896.70339584,   936.28627496,   995.33931158,  1021.67117306,
        1059.11275615,  1157.25565836,  1247.54257909,  1317.92608131,
        1372.58728901,  1462.29741387,  1497.97477993,  1691.87102916,
        1741.08729611,  1864.8769706 ,  2014.12979678,  2087.61827821,
        2451.0499158 ,  2493.34600895,  2676.34885371,  2692.64293513,
        2744.56022266,  2969.55198639,  3110.1388375 ,  3497.25654999,
        3564.97170851,  3598.20800884,  3693.26997004,  3710.11410812,
        3895.84128468])

















In [32]:



    


#xdb.xplot(spec)

Match a set of 5 gd_lines to detected lines





In [33]:



    


reload(desiboot)
id_dict = desiboot.id_arc_lines(pixpk,gd_lines,dlamb,wmark,line_guess=line_guess)#line_guess)



    











In [34]:



    


id_dict



    


















    
Out[34]:








{u'dlamb': 0.527,
 u'first_id_idx': [19, 22, 23, 25, 27],
 u'first_id_pix': array([ 1691.87102916,  2014.12979678,  2087.61827821,  2493.34600895,
         2692.64293513]),
 u'first_id_wave': array([ 6506.528 ,  6678.2766,  6717.043 ,  6929.4672,  7032.4128]),
 u'fit': {'coeff': array([ -7.27380983e+03,   4.46236920e-01,   1.86590813e-05]),
  'func': u'polynomial',
  'order': 2,
  'xmax': 1.0,
  'xmin': 0.0},
 'guess': 23,
 u'icen': 7,
 'im1': 22,
 'im2': 19,
 'ip1': 25,
 'ip2': 27,
 'rms': 0.069454506158127594,
 u'wmark': 6717.043}

















In [35]:



    


# IDs
plt.clf()
plt.figure(figsize=(18, 11))

yspec = np.log10(np.maximum(spec,1))
xpk = id_dict['first_id_pix']
xplt = np.arange(spec.size)
plt.plot(xplt,yspec,'b-')
plt.scatter(xpk,yspec[np.round(xpk).astype(int)],color='red')
# Guesses
for jj,xpixpk in enumerate(id_dict['first_id_pix']):
    plt.text(xpixpk, yspec[int(np.round(xpixpk))], '{:g}'.format(id_dict['first_id_wave'][jj]),
         ha='center',color='red',size='x-large')
#
plt.ylim(0.,np.max(yspec)*1.05)
plt.xlim(0,spec.size)
#plt.xlim(2650., 2750)
plt.xlabel('pixel')
plt.ylabel('log Flux')
plt.show()
plt.close()



    


















    








<matplotlib.figure.Figure at 0x12405a310>

Find the other good ones





In [36]:



    


desiboot.add_gdarc_lines(id_dict, pixpk, gd_lines)



    











In [37]:



    


id_dict



    


















    
Out[37]:








{u'dlamb': 0.527,
 u'first_fit': {'coeff': array([ -7.27380983e+03,   4.46236920e-01,   1.86590813e-05]),
  'func': u'polynomial',
  'order': 2,
  'xmax': 1.0,
  'xmin': 0.0},
 u'first_id_idx': [19, 22, 23, 25, 27],
 u'first_id_pix': array([ 1691.87102916,  2014.12979678,  2087.61827821,  2493.34600895,
         2692.64293513]),
 u'first_id_wave': array([ 6506.528 ,  6678.2766,  6717.043 ,  6929.4672,  7032.4128]),
 u'fit': {'coeff': array([ -1.06429692e+04,   1.19215407e+00,  -3.63755931e-05,
           1.35319491e-09]),
  'func': u'polynomial',
  'order': 3,
  'xmax': 1.0,
  'xmin': 0.0},
 'guess': 23,
 u'icen': 7,
 u'id_idx': [5, 11, 18, 19, 22, 23, 25, 27, 29, 30, 31],
 u'id_pix': [662.71309491414252,
  1021.6711730591247,
  1497.9747799297438,
  1691.8710291610926,
  2014.1297967788587,
  2087.6182782053656,
  2493.3460089453315,
  2692.6429351338438,
  2969.5519863947907,
  3110.1388375032575,
  3497.2565499931102],
 u'id_wave': [5944.8339999999998,
  6143.0619999999999,
  6402.2460000000001,
  6506.5280000000002,
  6678.2766000000001,
  6717.0429999999997,
  6929.4672,
  7032.4128000000001,
  7173.9380000000001,
  7245.1665000000003,
  7438.8980000000001],
 'im1': 22,
 'im2': 19,
 'ip1': 25,
 'ip2': 27,
 'rms': 0.07717510516330725,
 u'wmark': 6717.043}

















In [38]:



    


# IDs
plt.clf()
plt.figure(figsize=(18, 11))

yspec = np.log10(np.maximum(spec,1))
xpk = id_dict['id_pix']
xplt = np.arange(spec.size)
plt.plot(xplt,yspec,'b-')
plt.scatter(xpk,yspec[np.round(xpk).astype(int)],color='red')
# Guesses
for jj,xpixpk in enumerate(id_dict['id_pix']):
    plt.text(xpixpk, yspec[int(np.round(xpixpk))], '{:g}'.format(id_dict['id_wave'][jj]),
         ha='center',color='red')
#
plt.ylim(0.,np.max(yspec)*1.05)
plt.xlim(0,spec.size)
plt.xlabel('pixel')
plt.ylabel('log Flux')
plt.show()
plt.close()



    


















    








<matplotlib.figure.Figure at 0x123144c90>

Now the rest





In [39]:



    


desiboot.id_remainder(id_dict, pixpk, llist)



    











In [40]:



    


id_dict['rms']



    


















    
Out[40]:








0.07717510516330725

















In [41]:



    


# IDs
plt.clf()
plt.figure(figsize=(18, 11))

yspec = np.log10(np.maximum(spec,1))
xpk = id_dict['id_pix']
xplt = np.arange(spec.size)
plt.plot(xplt,yspec,'b-')
plt.scatter(xpk,yspec[np.round(xpk).astype(int)],color='red')
# Guesses
for jj,xpixpk in enumerate(id_dict['id_pix']):
    plt.text(xpixpk, yspec[int(np.round(xpixpk))], '{:g}'.format(id_dict['id_wave'][jj]),
         ha='center',color='red')
#
plt.ylim(0.,np.max(yspec)*1.05)
plt.xlim(0,spec.size)
plt.xlabel('pixel')
plt.ylabel('log Flux')
plt.show()
plt.close()



    


















    








<matplotlib.figure.Figure at 0x127e68250>

Final fit wave vs. pix too





In [42]:



    


final_fit, mask = dufits.iter_fit(np.array(id_dict['id_wave']), np.array(id_dict['id_pix']), 'polynomial', 3, xmin=0., xmax=1.)
rms = np.sqrt(np.mean((dufits.func_val(np.array(id_dict['id_wave'])[mask==0], final_fit)-np.array(id_dict['id_pix'])[mask==0])**2))
final_fit_pix,mask2 = dufits.iter_fit(np.array(id_dict['id_pix']), np.array(id_dict['id_wave']),'legendre',4, niter=5)



    











In [43]:



    


print("RMS = {:g}".format(rms))



    


















    






RMS = 0.156162

















In [ ]:

Content source: timahutchinson/desispec

Similar notebooks:

notebook.community | gallery | about