Merge pull request #1125 from desihub/improve-sky-sub

Improve sky subtraction

bin/desi_proc

@@ -291,8 +291,6 @@ if ( args.obstype in ['FLAT', 'TESTFLAT', 'SKY', 'TWILIGHT']     )   or \
                     cmd += " --degyx 2 --degyy 0"
                 if args.obstype in ['SCIENCE', 'SKY']:
                     cmd += ' --sky'
-                elif args.obstype in ['ARC', 'TESTARC']:
-                    cmd += ' --arc-lamps'
             else :
                 cmd = "ln -s {} {}".format(inpsf,outpsf)
             runcmd(cmd, inputs=[preprocfile, inpsf], outputs=[outpsf])
@@ -789,8 +787,8 @@ if args.obstype in ['SCIENCE', 'SKY'] and (not args.noskysub ) and (not args.nop
         cmd += " --o {}".format(skyfile)
         if args.no_extra_variance :
             cmd += " --no-extra-variance"
-        if args.adjust_sky_wavelength : cmd += " --adjust-wavelength"
-        if args.adjust_sky_lsf : cmd += " --adjust-lsf"
+        if not args.no_sky_wavelength_adjustment : cmd += " --adjust-wavelength"
+        if not args.no_sky_lsf_adjustment : cmd += " --adjust-lsf"
 
         runcmd(cmd, inputs=[framefile, fiberflatfile], outputs=[skyfile,])
 

py/desispec/scripts/specex.py

@@ -103,11 +103,23 @@ def main(args, comm=None):
     imgfile = args.input_image
     outfile = args.output_psf
 
+
+    nproc = 1
+    rank = 0
+    if comm is not None:
+        nproc = comm.size
+        rank = comm.rank
+
+    hdr=None
+    if rank == 0 :
+        hdr = fits.getheader(imgfile)
+    if comm is not None:
+        hdr = comm.bcast(hdr, root=0)
+
     if args.input_psf is not None:
         inpsffile = args.input_psf
     else:
         from desispec.calibfinder import findcalibfile
-        hdr = fits.getheader(imgfile)
         inpsffile = findcalibfile([hdr,], 'PSF')
 
     optarray = []
@@ -142,11 +154,6 @@ def main(args, comm=None):
 
     # Now we assign bundles to processes
 
-    nproc = 1
-    rank = 0
-    if comm is not None:
-        nproc = comm.size
-        rank = comm.rank
 
     mynbundle = int(nbundle / nproc)
     myfirstbundle = 0
@@ -184,6 +191,9 @@ def main(args, comm=None):
             if not os.path.isdir(outdir):
                 os.makedirs(outdir)
 
+    cam = hdr["camera"].lower().strip()
+    band = cam[0]
+
     failcount = 0
 
     for b in range(myfirstbundle, myfirstbundle+mynbundle):
@@ -197,7 +207,11 @@ def main(args, comm=None):
         com.extend(['--last-bundle', "{}".format(b)])
         com.extend(['--first-fiber', "{}".format(bspecmin[b])])
         com.extend(['--last-fiber', "{}".format(bspecmin[b]+bnspec[b]-1)])
-        com.extend(['--legendre-deg-wave', "{}".format(1)])
+        if band == "z" :
+            com.extend(['--legendre-deg-wave', "{}".format(3)])
+            com.extend(['--fit-continuum'])
+        else :
+            com.extend(['--legendre-deg-wave', "{}".format(1)])
         if args.broken_fibers :
             com.extend(['--broken-fibers', "{}".format(args.broken_fibers)])
         if args.debug :

py/desispec/sky.py

@@ -142,7 +142,7 @@ def _model_variance(frame,cskyflux,cskyivar,skyfibers) :
 
 
 
-def compute_uniform_sky(frame, nsig_clipping=4.,max_iterations=100,model_ivar=False,add_variance=True,adjust_wavelength=True,adjust_lsf=True) :
+def compute_uniform_sky(frame, nsig_clipping=4.,max_iterations=100,model_ivar=False,add_variance=True,adjust_wavelength=True,adjust_lsf=True,only_use_skyfibers_for_adjustments = True) :
     """Compute a sky model.
 
     Sky[fiber,i] = R[fiber,i,j] Flux[j]
@@ -165,6 +165,7 @@ def compute_uniform_sky(frame, nsig_clipping=4.,max_iterations=100,model_ivar=Fa
         add_variance : evaluate calibration error and add this to the sky model variance
         adjust_wavelength : adjust the wavelength of the sky model on sky lines to improve the sky subtraction
         adjust_lsf : adjust the LSF width of the sky model on sky lines to improve the sky subtraction
+        only_use_skyfibers_for_adjustments : interpolate adjustments using sky fibers only (else use median filter across fibers)
 
     returns SkyModel object with attributes wave, flux, ivar, mask
     """
@@ -492,12 +493,23 @@ def compute_uniform_sky(frame, nsig_clipping=4.,max_iterations=100,model_ivar=Fa
         nfibers_for_filter=10 # this number is a bit arbitrary/empirical.
         interpolated_sky_scale = scipy.ndimage.filters.median_filter(interpolated_sky_scale,(nfibers_for_filter,1))
         cskyflux = interpolated_sky_scale*cskyflux
+
+
+        allfibers=np.arange(frame.nspec)
         # the actual median filtering
         if adjust_wavelength :
-            interpolated_sky_dwave = scipy.ndimage.filters.median_filter(interpolated_sky_dwave,(nfibers_for_filter,1))
+            if only_use_skyfibers_for_adjustments : # simple interpolation over fibers
+                for j in range(interpolated_sky_dwave.shape[1]) :
+                    interpolated_sky_dwave[:,j] = np.interp(np.arange(interpolated_sky_dwave.shape[0]),skyfibers,interpolated_sky_dwave[skyfibers,j])
+            else : # median filter
+                interpolated_sky_dwave = scipy.ndimage.filters.median_filter(interpolated_sky_dwave,(nfibers_for_filter,1))
             cskyflux += interpolated_sky_dwave*dskydwave
-        if adjust_lsf :
-            interpolated_sky_dlsf = scipy.ndimage.filters.median_filter(interpolated_sky_dlsf,(nfibers_for_filter,1))
+        if adjust_lsf : # simple interpolation over fibers
+            if only_use_skyfibers_for_adjustments :
+                for j in range(interpolated_sky_dlsf.shape[1]) :
+                    interpolated_sky_dlsf[:,j] = np.interp(np.arange(interpolated_sky_dlsf.shape[0]),skyfibers,interpolated_sky_dlsf[skyfibers,j])
+            else : # median filter
+                interpolated_sky_dlsf = scipy.ndimage.filters.median_filter(interpolated_sky_dlsf,(nfibers_for_filter,1))
             cskyflux += interpolated_sky_dlsf*dskydlsf
 
     # look at chi2 per wavelength and increase sky variance to reach chi2/ndf=1

py/desispec/workflow/desi_proc_funcs.py

@@ -62,8 +62,8 @@ def get_shared_desi_proc_parser():
                         " use the most recent calibrations from *past* nights. If not set, uses default calibs instead.")
     parser.add_argument("--no-model-pixel-variance", action="store_true",
                         help="Do not use a model of the variance in preprocessing")
-    parser.add_argument("--adjust-sky-wavelength", action="store_true", help="Adjust wavelength of sky lines")
-    parser.add_argument("--adjust-sky-lsf", action="store_true", help="Adjust width of sky lines")
+    parser.add_argument("--no-sky-wavelength-adjustment", action="store_true", help="Do not adjust wavelength of sky lines")
+    parser.add_argument("--no-sky-lsf-adjustment", action="store_true", help="Do not adjust width of sky lines")
     parser.add_argument("--starttime", type=str, help='start time; use "--starttime `date +%%s`"')
     parser.add_argument("--timingfile", type=str, help='save runtime info to this json file; augment if pre-existing')