Merge pull request #2444 from desihub/traceshift2025b

Also scan for large trace shifts in continuum/LED exposures

Author: julienguy

bin/desi_ccd_xy

@@ -2,6 +2,7 @@
 
 import numpy as np
 import astropy.io.fits as pyfits
+from astropy.table import Table
 import specter.psf
 import sys
 import argparse
@@ -18,16 +19,29 @@ def readpsf(filename) :
         return specter.psf.SpotGridPSF(filename)
 
 parser = argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter)
-parser.add_argument('--psf', type = str, default = None, required = True,
-                    help = 'path of psf file')
+parser.add_argument('--psf', type = str, default = None, required = True, nargs = "*",
+                    help = 'path of psf files')
 parser.add_argument('--fiber', type = int, default = None, required = True,
                     help = 'fiber for psf1')
 parser.add_argument('--wavelength', type = float, default = 6000., required = False,
                     help = 'wavelength')
+parser.add_argument('-o','--outtable', type = str, default = None, required = False,
+                    help = 'output table')
 
 args        = parser.parse_args()
 
+xx=[]
+yy=[]
+for filename in args.psf :
+    psf=readpsf(filename)
+    x,y=psf.xy(args.fiber%500,args.wavelength)
+    print("x,y",x,y)
+    xx.append(x)
+    yy.append(y)
 
-psf=readpsf(args.psf)
-xy=psf.xy(args.fiber%500,args.wavelength)
-print("xy=",xy)
+if args.outtable is not None :
+    t=Table()
+    t["X"]=np.array(xx)
+    t["Y"]=np.array(yy)
+    print(t)
+    t.write(args.outtable,overwrite=True)

py/desispec/scripts/proc.py

@@ -558,7 +558,7 @@ def main(args=None, comm=None):
                     cmd += " --psf {}".format(inpsf)
                     cmd += " --degxx 2 --degxy 0"
                     if args.obstype in ['FLAT', 'TESTFLAT', 'TWILIGHT'] :
-                        cmd += " --continuum"
+                        cmd += " --continuum --no-large-shift-scan"
                     else :
                         cmd += " --degyx 2 --degyy 0"
                     if args.obstype in ['SCIENCE', 'SKY']:

py/desispec/scripts/trace_shifts.py

@@ -22,7 +22,7 @@
 from desispec.io import shorten_filename
 from desiutil.log import get_logger
 from desispec.util import parse_int_args
-from desispec.trace_shifts import write_traces_in_psf,compute_dx_from_cross_dispersion_profiles,compute_dy_from_spectral_cross_correlation,monomials,polynomial_fit,compute_dy_using_boxcar_extraction,compute_dx_dy_using_psf,shift_ycoef_using_external_spectrum,recompute_legendre_coefficients,recompute_legendre_coefficients_for_x,recompute_legendre_coefficients_for_y,list_of_expected_spot_positions
+from desispec.trace_shifts import write_traces_in_psf,compute_dx_from_cross_dispersion_profiles,compute_dy_from_spectral_cross_correlation,monomials,polynomial_fit,compute_dy_using_boxcar_extraction,compute_dx_dy_using_psf,shift_ycoef_using_external_spectrum,recompute_legendre_coefficients,recompute_legendre_coefficients_for_x,recompute_legendre_coefficients_for_y,list_of_expected_spot_positions,compute_x_offset_from_central_band_cross_dispersion_profile
 from desispec.large_trace_shifts import detect_spots_in_image,match_same_system
 
 def parse(options=None):
@@ -60,7 +60,7 @@ def parse(options=None):
     parser.add_argument('--degyy', type = int, default = 0, required=False,
                         help = 'polynomial degree for y shifts along y')
     parser.add_argument('--continuum', action='store_true',
-                        help = 'only fit shifts along x for continuum input image')
+                        help = 'only fit shifts along x for continuum or LED input image')
     parser.add_argument('--auto', action='store_true',
                         help = 'choose best method (sky,continuum or just internal calib) from the FLAVOR keyword in the input image header')
 
@@ -72,7 +72,8 @@ def parse(options=None):
                         help = "width of cross-dispersion profile")
     parser.add_argument('--ccd-rows-rebin', type = int, default = 4 , required=False,
                         help = "rebinning of CCD rows to run faster")
-
+    parser.add_argument('--no-large-shift-scan', action="store_true",
+                        help = "do not perform a large shift scan for arc lamp or continuum exposures")
     args = parser.parse_args(options)
 
     return args
@@ -459,15 +460,21 @@ def fit_trace_shifts(image, args):
             args.sky = True
         log.info("wavelength calib, internal={}, sky={} , arc_lamps={}".format(internal_wavelength_calib,args.sky,args.arc_lamps))
 
-    if args.arc_lamps :
+    cfinder = CalibFinder([image.meta])
+    fibers=np.arange(tset.nspec)
+    if cfinder.haskey("BROKENFIBERS") :
+        brokenfibers=parse_int_args(cfinder.value("BROKENFIBERS"))%500
+        log.debug(f"brokenfibers={brokenfibers}")
+        fibers=fibers[np.isin(fibers, brokenfibers, invert=True)]
+
+    # for arc lamp images, we can predict where we expect to see spots on the CCD image
+    # given an input traceset (variable tset), and use the comparison between the prediction
+    # and actual spots locations to derive a first correction to the coordinates saved in the traceset
+    # this is disabled by the option --no-large-shift-scan
+    if args.arc_lamps and (not args.no_large_shift_scan) :
 
         log.info("for arc lamps, find a first solution by comparing expected spots positions with detections over the whole CCD")
-        cfinder = CalibFinder([image.meta])
-        fibers=np.arange(tset.nspec)
-        if cfinder.haskey("BROKENFIBERS") :
-            brokenfibers=parse_int_args(cfinder.value("BROKENFIBERS"))%500
-            log.debug(f"brokenfibers={brokenfibers}")
-            fibers=fibers[np.isin(fibers, brokenfibers, invert=True)]
+
         xref,yref = list_of_expected_spot_positions(tset,fibers)
         xin,yin   = detect_spots_in_image(image)
 
@@ -510,6 +517,16 @@ def fit_trace_shifts(image, args):
         xcoef[:,0] += delta_xref
         ycoef[:,0] += delta_yref
 
+    # for continuum images, we can predict where we expect to see spectral traces on the CCD image
+    # given an input traceset (variable tset), and use the comparison between the prediction
+    # and actual spectral traces to derive a first correction to the X coordinates saved in the traceset
+    # this is disabled by the option --no-large-shift-scan
+    if args.continuum  and (not args.no_large_shift_scan) :
+        log.info("for continuum or LED exposures, find a first solution on delta_X by comparing a wide cross-dispersion profile with expectations")
+        delta_xref = compute_x_offset_from_central_band_cross_dispersion_profile(tset, image, fibers=fibers)
+        log.info(f"apply best shift delta x = {delta_xref} to traceset")
+        xcoef[:,0] += delta_xref
+
     spectrum_filename = args.spectrum
     continuum_subtract = False
     if args.sky :
@@ -557,6 +574,7 @@ def fit_trace_shifts(image, args):
 
     else :
 
+
         # internal calibration method that does not use the psf
         # nor a prior set of lines. this method is much faster
 

py/desispec/trace_shifts.py

@@ -1507,3 +1507,42 @@ def list_of_expected_spot_positions(traceset,fibers=None,max_number_of_lines=50)
     xspot=np.hstack(xspot)
     yspot=np.hstack(yspot)
     return xspot,yspot
+
+
+def compute_x_offset_from_central_band_cross_dispersion_profile(tset, image, fibers=None,halfwidth=50) :
+    log=get_logger()
+    bands=7 # measure delta_x in 7 horizontal bands of 100 pixel wide each and return the median offset
+    bandwidth=100
+    hb=bands//2
+    n0=image.pix.shape[0]
+    n1=image.pix.shape[1]
+    ivar=image.ivar*(image.mask==0)
+    if fibers is None :
+        fibers = np.arange(tset.nspec)
+    delta_x_array=[]
+    for b in range(-hb,hb+1) :
+        begin=int(n0//2+(b-0.5)*bandwidth)
+        end=begin+bandwidth
+
+        sw=np.sum(ivar[begin:end,:],axis=0)
+        swf=np.sum(image.pix[begin:end,:]*ivar[begin:end,:])
+        prof=(swf/(sw+(sw==0)))
+        y=int(n0//2+b*bandwidth)
+        xc=np.array([tset.x_vs_y(fiber,y) for fiber in fibers])
+        oversample=4
+        ox=np.arange(n1*oversample)/oversample
+        xx=np.arange(n1)
+        sigma=1.5 # pixel, approximate
+        mprof=np.exp(-(ox[:,None]-xc[None,:])**2/2./sigma**2).sum(axis=1)
+        norm = 1./np.sqrt(np.sum(prof**2)*(np.sum(mprof**2)/oversample))
+        ddx=np.linspace(-halfwidth,halfwidth,4*halfwidth+1)
+        xcorr=np.zeros(ddx.size)
+        for i,dx in enumerate(ddx) :
+            xcorr[i]=np.sum(prof*np.interp(xx,ox+dx,mprof))
+        xcorr *= norm
+        imax=np.argmax(xcorr)
+        log.info("horizontal band {}:{} delta x = {:.1f} pixels".format(begin,end,ddx[imax]))
+        delta_x_array.append(ddx[imax])
+    delta_x = np.median(delta_x_array)
+    log.info("median delta x = {:.1f} pixels".format(delta_x))
+    return delta_x