doppler mods for apogee work

.gitignore

@@ -0,0 +1 @@
+*.pyc

doppler/__init__.py

@@ -13,8 +13,8 @@
 models = cannon.load_models()
 cannon.models = models
 
-def read(filename=None,format=None):
-    return reader.read(filename=filename,format=None)
+def read(filename=None,format=None,badval=None):
+    return reader.read(filename=filename,format=None,badval=None)
 
 def fit(*args,**kwargs):
     return rv.fit(*args,**kwargs)

doppler/cannon.py

@@ -1146,7 +1146,10 @@ def prepare_cannon_model(model,spec,dointerp=False):
             #  using 3 affects the final profile shape
             nbin = np.round(np.min(fwhmpix)//4).astype(int)
             if nbin==0:
-                raise Exception('Model has lower resolution than the observed spectrum')
+                import pdb
+                nbin=1
+                print(spec.filename, 'Model has lower resolution than the observed spectrum',fwhmpix.min())
+                #raise Exception('Model has lower resolution than the observed spectrum',spec.filename,fwhmpix.min())
             if nbin>1:
                 rmodel = rebin_cannon_model(tmodel,nbin)
                 rmodel.rebin = True

doppler/cannongrid_cubic_norm_training.py

@@ -0,0 +1,88 @@
+#!/usr/bin/env python
+
+# Imports
+import os
+import thecannon as tc
+from astropy.io import fits
+from astropy.table import Table
+import numpy as np
+import matplotlib.pyplot as plt
+import pickle
+
+#tag = '3000_18000_allstars'
+#tag = '3000_18000_rgb'
+#tag = '3000_18000_cooldwarfs'
+tag = '3000_18000_hotstars'
+#tag = '3000_18000_hotdwarfs'
+#tag = '3000_18000_coolstars'
+#tag = '3000_18000_whitedwarfs'
+
+
+print(tag)
+
+# Import the spectra and labels
+normalized_flux = fits.getdata('cannongrid_'+tag+'_synth_data_norm.fits.gz')
+labelled_set=Table.read('cannongrid_'+tag+'_synth_pars.fits')
+
+# Add wavelengths
+nspec, npix = normalized_flux.shape
+wave = np.arange(npix)*0.10+3000.0
+#model3.dispersion = wave
+normalized_ivar = normalized_flux.copy()*0 + 1e4
+vec3 = tc.vectorizer.PolynomialVectorizer(labelled_set.colnames, 3)
+model3 = tc.CannonModel(labelled_set, normalized_flux, normalized_ivar, vec3, wave)
+model3.regularization = 0   # no regularization for now
+
+# Train the model
+nr_theta, nr_s2, nr_metadata = model3.train()  
+if os.path.exists('cannongrid_'+tag+'_norm_cubic_model.pkl'): os.remove('cannongrid_'+tag+'_norm_cubic_model.pkl')
+model3.write('cannongrid_'+tag+'_norm_cubic_model.pkl')
+
+# Check if there is a continuum model to add
+contfile = 'cannongrid_'+tag+'_cont_cubic_model_logflux.pkl'
+if os.path.exists(contfile):
+    print('Continuum model found.  Adding it')
+    f = open(contfile,'rb')
+    cont = pickle.load(f)
+    f.close()
+else:
+    cont = None
+
+# This adds my custom attributes
+#   model.write() doesn't same them
+infile = open('cannongrid_'+tag+'_norm_cubic_model.pkl','rb')
+temp = pickle.load(infile)
+infile.close()
+temp['fwhm'] = 0.001
+temp['wavevac'] = False
+if cont is not None:
+    temp['continuum'] = cont
+md = temp['metadata']
+ta = md['trained_attributes']
+if cont is not None:
+    ta = ta+('fwhm','wavevac','continuum')
+else:
+    ta = ta+('fwhm','wavevac')
+md['trained_attributes'] = ta
+temp['metadata'] = md
+os.remove('cannongrid_'+tag+'_norm_cubic_model.pkl')
+outfile = open('cannongrid_'+tag+'_norm_cubic_model.pkl','wb')
+pickle.dump(temp,outfile)
+outfile.close()
+
+# Test/fit the spectra to get the labels
+labels3, cov, meta = model3.test(normalized_flux, normalized_ivar)
+
+
+# one-to-one comparisons
+names = model3.vectorizer.label_names
+fig = plt.figure(figsize=(15,10))
+for i in range(0, 3):
+    plt.subplot(2, 2, i+1)
+    plt.plot(labelled_set[names[i]],labels3[:,i],'+')
+    plt.plot(labelled_set[names[i]],labelled_set[names[i]])
+    plt.xlabel('Input '+names[i])
+    plt.ylabel('Recovered '+names[i])
+fig.savefig('cannongrid_'+tag+'_comparison.png')    
+
+

doppler/lsf.py

@@ -118,7 +118,7 @@ def ghlsf(x,xcenter,params,nowings=False):
     wingparams2 = np.empty((nlsf*npix,params['nWpar']+1))
     for i in range(params['nWpar']+1):
         wingparams2[:,i] = np.repeat(wingparams1[i,:],nlsf)
-    out += ghwingsbin(xlsf,wingparams2,params['binsize'],params['Wproftype'])
+    if not nowings : out += ghwingsbin(xlsf,wingparams2,params['binsize'],params['Wproftype'])
 
     # Reshape it to [Npix,Nlsf]
     out = out.reshape(npix,nlsf)
@@ -528,7 +528,7 @@ def ghlsf_bovy(x,xcenter,params,nowings=False):
     # Calculate the GH part of the LSF
     out = gausshermitebin_bovy(x,ghparams,params['binsize'])
     # Calculate the Wing part of the LSF
-    out += ghwingsbin_bovy(x,wingparams,params['binsize'],params['Wproftype'])
+    if not nowings: out += ghwingsbin_bovy(x,wingparams,params['binsize'],params['Wproftype'])
 
     return out
 

doppler/reader.py

@@ -17,14 +17,16 @@
 from scipy.ndimage.filters import median_filter
 from dlnpyutils import utils as dln, bindata
 from .spec1d import Spec1D
+import matplotlib.pyplot as plt
+import pdb
 
 # Ignore these warnings, it's a bug
 warnings.filterwarnings("ignore", message="numpy.dtype size changed")
 warnings.filterwarnings("ignore", message="numpy.ufunc size changed")
 
 
 # Load a spectrum
-def read(filename=None,format=None):
+def read(filename=None,format=None,badval=None):
     '''
     This reads in a SDSS-IV MWM training set spectrum and returns an
     object that is guaranteed to have certain information.
@@ -80,7 +82,7 @@ def read(filename=None,format=None):
         if format.lower() not in _readers.keys():
             raise ValueError('reader '+format+' not found')
         # Use the requested reader/format
-        out = _readers[format](filename)
+        out = _readers[format](filename,badval=badval)
         if out is not None: return out
 
     # Loop over all readers until we get a spectrum out
@@ -97,7 +99,7 @@ def read(filename=None,format=None):
 
 
 # Load APOGEE apVisit/asVisit spectra
-def apvisit(filename):
+def apvisit(filename,badval=20735):
     """
     Read a SDSS APOGEE apVisit spectrum.
 
@@ -119,7 +121,7 @@ def apvisit(filename):
     """
 
     base, ext = os.path.splitext(os.path.basename(filename))
-    
+
     # APOGEE apVisit, visit-level spectrum
     if (base.find("apVisit") > -1) | (base.find("asVisit") > -1):
         # HISTORY AP1DVISIT:  HDU0 = Header only                                          
@@ -165,17 +167,25 @@ def apvisit(filename):
         #           'LITTROW_GHOST','PERSIST_HIGH','PERSIST_MED','PERSIST_LOW','SIG_SKYLINE','SIG_TELLURIC','NOT_ENOUGH_PSF','']
         #   badflag = [1,1,1,1,1,1,1,1,
         #              0,0,0,0,0,0,1,0]
-        mask = (np.bitwise_and(spec.bitmask,16639)!=0) | (np.isfinite(spec.flux)==False)
+        #mask = (np.bitwise_and(spec.bitmask,16639)!=0) | (np.isfinite(spec.flux)==False)
+        mask = (np.bitwise_and(spec.bitmask,badval)!=0) | (np.isfinite(spec.flux)==False)
         # Extra masking for bright skylines
-        x = np.arange(spec.npix)
-        nsky = 4
-        for i in range(spec.norder):
-            sky = spec.sky[:,i]
-            medsky = median_filter(sky,201,mode='reflect')
-            medcoef = dln.poly_fit(x,medsky/np.median(medsky),2)
-            medsky2 = dln.poly(x,medcoef)*np.median(medsky)
-            skymask1 = (sky>nsky*medsky2)    # pixels Nsig above median sky
-            mask[:,i] = np.logical_or(mask[:,i],skymask1)    # OR combine
+        # Commented out in favor of using SIG_SKYLINE in bitmask
+        # This can also mask too many pixels
+        #x = np.arange(spec.npix)
+        #nsky = 4
+        ##plt.clf()
+        #for i in range(spec.norder):
+        #    sky = spec.sky[:,i]
+        #    medsky = median_filter(sky,201,mode='reflect')
+        #    medcoef = dln.poly_fit(x,medsky/np.nanmedian(medsky),2)
+        #    medsky2 = dln.poly(x,medcoef)*np.nanmedian(medsky)
+        #    skymask1 = (sky>nsky*medsky2)    # pixels Nsig above median sky
+        #    #mask[:,i] = np.logical_or(mask[:,i],skymask1)    # OR combine
+        #    #plt.plot(spec.wave[:,i],sky)
+        #    #plt.plot(spec.wave[:,i],nsky*medsky2)
+        #    #plt.plot(spec.wave[:,i],spec.flux[:,i])
+        ##plt.draw()
         spec.mask = mask
         # Fix NaN pixels
         for i in range(spec.norder):
@@ -197,7 +207,7 @@ def apvisit(filename):
 
 
 # Load APOGEE apStar/asStar spectra    
-def apstar(filename):
+def apstar(filename,badval=20735):
     """
     Read an SDSS APOGEE apStar spectrum.
 
@@ -281,7 +291,7 @@ def apstar(filename):
         #           'LITTROW_GHOST','PERSIST_HIGH','PERSIST_MED','PERSIST_LOW','SIG_SKYLINE','SIG_TELLURIC','NOT_ENOUGH_PSF','']
         #   badflag = [1,1,1,1,1,1,1,1,
         #              0,0,0,0,0,0,1,0]
-        mask = (np.bitwise_and(spec.bitmask,16639)!=0) | (np.isfinite(spec.flux)==False)
+        mask = (np.bitwise_and(spec.bitmask,badval)!=0) | (np.isfinite(spec.flux)==False)
         # Extra masking for bright skylines
         x = np.arange(spec.npix)
         nsky = 4
@@ -309,7 +319,7 @@ def apstar(filename):
 
 
 # Load SDSS BOSS spectra    
-def boss(filename):
+def boss(filename,badval=0):
     """
     Read a SDSS BOSS spectrum.
 
@@ -377,7 +387,7 @@ def boss(filename):
 
 
 # Load SDSS MaStar spectra
-def mastar(filename):
+def mastar(filename,badval=0):
     """
     Read a SDSS MaStar spectrum.
 
@@ -441,7 +451,7 @@ def mastar(filename):
 
 
 # Load IRAF-style spectra
-def iraf(filename):
+def iraf(filename,badval=0):
     """
     Read an IRAF-style spectrum.