Bench

config/auxtel.ini

@@ -29,7 +29,7 @@ OBS_QUANTUM_EFFICIENCY = calexp_2020031500162-EMPTY_ronchi90lpmm-det000_auxtel_t
 OBS_CAMERA_ROTATION = 0
 # Camera (x,y) flip signs with respect to (north-up, east-left) system
 OBS_CAMERA_DEC_FLIP_SIGN = 1
-OBS_CAMERA_RA_FLIP_SIGN = 1
+OBS_CAMERA_RA_FLIP_SIGN = -1
 
 [CCD]
 # size of the image in pixel
@@ -43,17 +43,17 @@ CCD_MAXADU = 170000
 # electronic gain : elec/ADU
 CCD_GAIN = 1.1
 # rebinning of the image in pixel
-CCD_REBIN = 4
+CCD_REBIN = 2
 
 [dispersers]
 # distance between hologram and CCD in mm
 DISTANCE2CCD = 175
 # uncertainty on distance between hologram and CCD in mm
-DISTANCE2CCD_ERR = 0.5
+DISTANCE2CCD_ERR = 0.75
 # constructor wavelength to make holograms in mm
 LAMBDA_CONSTRUCTOR = 639e-6
 # approximate effective number of lines per millimeter of the hologram
-GROOVES_PER_MM = 100
+GROOVES_PER_MM = 150
 # plate center shift on x in mm in filter frame
 PLATE_CENTER_SHIFT_X = -6.
 # plate center shift on x in mm in filter frame

config/auxtel_quicklook.ini

@@ -29,7 +29,7 @@ OBS_QUANTUM_EFFICIENCY = calexp_2020031500162-EMPTY_ronchi90lpmm-det000_auxtel_t
 OBS_CAMERA_ROTATION = 0
 # Camera (x,y) flip signs with respect to (north-up, east-left) system
 OBS_CAMERA_DEC_FLIP_SIGN = 1
-OBS_CAMERA_RA_FLIP_SIGN = 1
+OBS_CAMERA_RA_FLIP_SIGN = -1
 
 [CCD]
 # size of the image in pixel
@@ -49,11 +49,11 @@ CCD_REBIN = 1
 # distance between hologram and CCD in mm
 DISTANCE2CCD = 175
 # uncertainty on distance between hologram and CCD in mm
-DISTANCE2CCD_ERR = 0.5
+DISTANCE2CCD_ERR = 0.75
 # constructor wavelength to make holograms in mm
 LAMBDA_CONSTRUCTOR = 639e-6
 # approximate effective number of lines per millimeter of the hologram
-GROOVES_PER_MM = 100
+GROOVES_PER_MM = 150
 # plate center shift on x in mm in filter frame
 PLATE_CENTER_SHIFT_X = -6.
 # plate center shift on x in mm in filter frame
@@ -94,7 +94,7 @@ LAMBDA_MAX = 1100
 
 [background subtraction parameters]
 # half transverse width of the signal rectangular window in pixels
-PIXWIDTH_SIGNAL = 20
+PIXWIDTH_SIGNAL = 40
 # distance from dispersion axis to analyse the background in pixels
 PIXDIST_BACKGROUND = 140
 # transverse width of the background rectangular window in pixels

config/lpnhe.ini

@@ -39,9 +39,9 @@ CCD_GAIN = 0.85
 
 [dispersers]
 # distance between hologram and CCD in mm
-DISTANCE2CCD = 198
+DISTANCE2CCD = 204.2
 # uncertainty on distance between hologram and CCD in mm
-DISTANCE2CCD_ERR = 1
+DISTANCE2CCD_ERR = 0.5
 # constructor wavelength to make holograms in mm
 LAMBDA_CONSTRUCTOR = 639e-6
 # approximate effective number of lines per millimeter of the hologram
@@ -84,7 +84,7 @@ LAMBDA_MAX = 1000
 
 [background subtraction parameters]
 # half transverse width of the signal rectangular window in pixels
-PIXWIDTH_SIGNAL = 10
+PIXWIDTH_SIGNAL = 20
 # distance from dispersion axis to analyse the background in pixels
 PIXDIST_BACKGROUND = 80
 # transverse width of the background rectangular window in pixels
@@ -104,15 +104,15 @@ PSF_PIXEL_STEP_TRANSVERSE_FIT = 10
 
 [detection line algorithm parameters]
 # order of the background polynome to fit
-CALIB_BGD_ORDER = 3
+CALIB_BGD_ORDER = 1
 # half range to look for local extrema in pixels around tabulated line values
-CALIB_PEAK_WIDTH = 3
+CALIB_PEAK_WIDTH = 25
 # size of the peak sides to use to fit spectrum base line
 CALIB_BGD_WIDTH = 40
 # window size for the savgol filter in pixels
-CALIB_SAVGOL_WINDOW = 5
+CALIB_SAVGOL_WINDOW = 7
 # polynom order for the savgol filter
-CALIB_SAVGOL_ORDER = 2
+CALIB_SAVGOL_ORDER = 1
 
 [plot settings]
 # paper plot style

runExtractor.py

@@ -1,6 +1,7 @@
 from spectractor import parameters
 from spectractor.extractor.extractor import Spectractor
 from spectractor.logbook import LogBook
+from spectractor.config import load_config
 import sys
 
 if __name__ == "__main__":
@@ -35,21 +36,27 @@
 
     file_names = args.input
 
+    load_config(args.config)
+
     logbook = LogBook(logbook=args.logbook)
     for file_name in file_names:
         disperser_label = args.disperser_label
-        if args.target_xy == "0,0" and args.target_label == "":
+        if parameters.OBS_NAME == "CTIO":
             tag = file_name.split('/')[-1]
             tag = tag.replace('sim_', 'reduc_')
             disperser_label, target_label, xpos, ypos = logbook.search_for_image(tag)
+            guess = [xpos, ypos]
             if target_label is None or xpos is None or ypos is None:
                 continue
         else:
-            xpos, ypos = args.target_xy.split(",")
+            guess = None
+            if args.target_xy != "0,0":
+                xpos, ypos = args.target_xy.split(",")
+                xpos = float(xpos)
+                ypos = float(ypos)
+                guess = [xpos, ypos]
             target_label = args.target_label
-            xpos = float(xpos)
-            ypos = float(ypos)
-            if target_label == "" or (xpos == 0 and ypos == 0):
-                sys.exit("Options --xy and --target must be used together, one of these seems not set.")
-        Spectractor(file_name, args.output_directory, target_label=target_label, guess=[xpos, ypos],
+            # if target_label == "" or (xpos == 0 and ypos == 0):
+            #     sys.exit("Options --xy and --target must be used together, one of these seems not set.")
+        Spectractor(file_name, args.output_directory, target_label=target_label, guess=guess,
                     disperser_label=disperser_label, config=args.config)

spectractor/config.py

@@ -67,11 +67,12 @@ def load_config(config_filename):
     .. doctest:
         :hide:
 
-        >>> load_config("./config/unknown_file.ini")
+        >>> load_config("./config/ctio.ini")
+        >>> load_config("ctio.ini")
+        >>> load_config("./config/unknown_file.ini")  #doctest: +ELLIPSIS
         Traceback (most recent call last):
         ...
-        SystemExit: Config file ./config/unknown_file.ini does not exist.
-        >>> load_config("./config/ctio.ini")
+        FileNotFoundError: Config file ./config/unknown_file.ini does not exist.
 
     """
     if not os.path.isfile(os.path.join(parameters.CONFIG_DIR, "default.ini")):
@@ -82,7 +83,10 @@ def load_config(config_filename):
     from_config_to_parameters(config)
 
     if not os.path.isfile(config_filename):
-        raise FileNotFoundError(f'Config file {config_filename} does not exist.')
+        if not os.path.isfile(os.path.join(parameters.CONFIG_DIR, config_filename)):
+            raise FileNotFoundError(f'Config file {config_filename} does not exist.')
+        else:
+            config_filename = os.path.join(parameters.CONFIG_DIR, config_filename)
     # Load the configuration file
     config = configparser.ConfigParser()
     config.read(config_filename)
@@ -164,4 +168,3 @@ def set_logger(logger):
     import doctest
 
     doctest.testmod()
-

spectractor/extractor/dispersers/holo4_003/N.txt

@@ -1,2 +1,2 @@
-156.49452269170578
--0.275
+150
+1.0

spectractor/extractor/dispersers/holo4_003/hologram_center.txt

@@ -1,2 +1,2 @@
 #x_star y_star grooves_per_mm grooves_per_mm_err 
-400.0 400.0 0.5
+400.0 2750.0 -0.28

spectractor/extractor/extractor.py

@@ -500,8 +500,8 @@ def plot_spectrogram_comparison_simple(self, ax, title='', extent=None, dispersi
             ax[1, 1].get_yaxis().set_label_coords(3.5, 0.5)
             ax[2, 1].get_yaxis().set_label_coords(3.5, 0.5)
             ax[3, 1].remove()
-            ax[3, 0].plot(self.lambdas[sub], data.sum(axis=0)[sub], label='Data')
-            ax[3, 0].plot(self.lambdas[sub], model.sum(axis=0)[sub], label='Model')
+            ax[3, 0].plot(self.lambdas[sub], np.nansum(data, axis=0)[sub], label='Data')
+            ax[3, 0].plot(self.lambdas[sub], np.nansum(model, axis=0)[sub], label='Model')
             ax[3, 0].set_ylabel('Cross spectrum')
             ax[3, 0].set_xlabel(r'$\lambda$ [nm]')
             ax[3, 0].legend(fontsize=7)
@@ -765,6 +765,12 @@ def Spectractor(file_name, output_directory, target_label, guess=None, disperser
     image = Image(file_name, target_label=target_label, disperser_label=disperser_label)
     if guess is not None and image.target_guess is None:
         image.target_guess = np.asarray(guess)
+    if image.target_guess is None:
+        from scipy.signal import medfilt2d
+        data = medfilt2d(image.data.T, kernel_size=9)
+        image.target_guess = np.unravel_index(np.argmax(data), data.shape)
+        my_logger.info(f"\n\tNo guess position of order 0 has been given. Assuming the spectrum to extract comes "
+                       f"from the brightest object, guess position is set as {image.target_guess}.")
     if parameters.DEBUG:
         image.plot_image(scale='symlog', target_pixcoords=image.target_guess)
 

spectractor/extractor/images.py

@@ -14,6 +14,7 @@
 from spectractor.extractor.dispersers import Hologram
 from spectractor.extractor.psf import Moffat
 from spectractor.simulation.adr import hadec2zdpar
+from spectractor.simulation.throughput import TelescopeTransmission
 from spectractor.tools import (plot_image_simple, save_fits, load_fits, fit_poly1d, plot_compass_simple,
                                fit_poly1d_outlier_removal, weighted_avg_and_std,
                                fit_poly2d_outlier_removal, hessian_and_theta,
@@ -68,7 +69,7 @@ def __init__(self, file_name, target_label="", disperser_label="", config=""):
         self.disperser_label = disperser_label
         self.target_label = target_label
         self.target_guess = None
-        self.filter = None
+        self.filter_label = ""
         self.filters = None
         self.header = None
         self.data = None
@@ -111,6 +112,10 @@ def __init__(self, file_name, target_label="", disperser_label="", config=""):
         self.header['D2CCD'] = parameters.DISTANCE2CCD
         self.header.comments["D2CCD"] = "[mm] distance between disperser and CCD"
 
+        if self.filter_label != "" and "empty" not in self.filter_label.lower():
+            t = TelescopeTransmission(filter_label=self.filter_label)
+            t.reset_lambda_range(transmission_threshold=1e-4)
+
         if self.target_label != "":
             self.target = load_target(self.target_label, verbose=parameters.VERBOSE)
             self.header['REDSHIFT'] = str(self.target.redshift)
@@ -440,7 +445,8 @@ def load_CTIO_image(image):
     image.airmass = float(image.header['AIRMASS'])
     image.expo = float(image.header['EXPTIME'])
     image.filters = image.header['FILTERS']
-    image.filter = image.header['FILTER1']
+    if "dia" not in image.header['FILTER1'].lower():
+        image.filter_label = image.header['FILTER1']
     image.disperser_label = image.header['FILTER2']
     image.ra = Angle(image.header['RA'], unit="hourangle")
     image.dec = Angle(image.header['DEC'], unit="deg")
@@ -556,6 +562,18 @@ def load_LPNHE_image(image):  # pragma: no cover
     image.header.comments["ZPOS"] = hdus["XYZ"].header.comments["ZPOS"]
     image.my_logger.info('\n\tImage loaded')
 
+    if image.target_label == "HG-AR":
+        parameters.OBS_OBJECT_TYPE = "HG-AR"
+    elif "QTH" in image.target_label:
+        parameters.OBS_OBJECT_TYPE = "MONOCHROMATOR"
+        parameters.PIXWIDTH_SIGNAL = 30
+        image.target_label += str(hdus["CORNERSTONE"].header["WVLGTH"])
+    elif "LASER" in image.target_label:
+        parameters.OBS_OBJECT_TYPE = "LASER"
+        parameters.PIXWIDTH_SIGNAL = 30
+    else:
+        raise ValueError(f"\nTarget label {image.target_label} unknown for LPNHE optical bench.")
+
 
 def load_AUXTEL_image(image):  # pragma: no cover
     """Specific routine to load AUXTEL fits files and load their data and properties for Spectractor.
@@ -572,6 +590,8 @@ def load_AUXTEL_image(image):  # pragma: no cover
     hdu_list.close()  # need to free allocation for file descripto
     image.date_obs = image.header['DATE']
     image.expo = float(image.header['EXPTIME'])
+    if "empty" not in image.header['FILTER'].lower():
+        image.filter_label = image.header['FILTER']
     # transformations so that stars are like in Stellarium up to a rotation
     # with spectrogram nearly horizontal and on the right of central star
     image.data = image.data.T[::-1, ::-1]
@@ -586,31 +606,41 @@ def load_AUXTEL_image(image):  # pragma: no cover
     image.disperser_label = image.header['GRATING']
     image.ra = Angle(image.header['RA'], unit="deg")
     image.dec = Angle(image.header['DEC'], unit="deg")
-    image.hour_angle = Angle(image.header['HA'], unit="deg")
-    image.temperature = 10  # image.header['OUTTEMP']
-    image.pressure = 730  # image.header['OUTPRESS']
-    image.humidity = 25  # image.header['OUTHUM']
+    image.hour_angle = Angle(image.header['HASTART'], unit="hourangle")
+    if 'AIRTEMP' in image.header:
+        image.temperature = image.header['AIRTEMP']
+    else:
+        image.temperature = 10
+    if 'PRESSURE' in image.header:
+        image.pressure = image.header['PRESSURE']
+    else:
+        image.pressure = 743
+    if 'HUMIDITY' in image.header:
+        image.humidity = image.header['HUMIDITY']
+    else:
+        image.humidity = 40
     if 'adu' in image.header['BUNIT']:
         image.units = 'ADU'
     parameters.OBS_CAMERA_ROTATION = 90 - float(image.header["ROTPA"])
     if parameters.OBS_CAMERA_ROTATION > 360:
         parameters.OBS_CAMERA_ROTATION -= 360
     if parameters.OBS_CAMERA_ROTATION < -360:
         parameters.OBS_CAMERA_ROTATION += 360
-    rotation_wcs = 180 / np.pi * np.arctan2(hdu_list[1].header["CD2_1"], hdu_list[1].header["CD1_1"]) + 90
-    if not np.isclose(rotation_wcs % 360, parameters.OBS_CAMERA_ROTATION % 360, atol=2):
-        image.my_logger.warning(f"\n\tWCS rotation angle is {rotation_wcs} degree while "
-                                f"parameters.OBS_CAMERA_ROTATION={parameters.OBS_CAMERA_ROTATION} degree. "
-                                f"\nBoth differs by more than 2 degree... bug ?")
+    if "CD2_1" in hdu_list[1].header:
+        rotation_wcs = 180 / np.pi * np.arctan2(hdu_list[1].header["CD2_1"], hdu_list[1].header["CD1_1"]) + 90
+        if not np.isclose(rotation_wcs % 360, parameters.OBS_CAMERA_ROTATION % 360, atol=2):
+            image.my_logger.warning(f"\n\tWCS rotation angle is {rotation_wcs} degree while "
+                                    f"parameters.OBS_CAMERA_ROTATION={parameters.OBS_CAMERA_ROTATION} degree. "
+                                    f"\nBoth differs by more than 2 degree... bug ?")
     parameters.OBS_ALTITUDE = float(image.header['OBS-ELEV']) / 1000
     parameters.OBS_LATITUDE = image.header['OBS-LAT']
     image.read_out_noise = 8.5 * np.ones_like(image.data)
     image.target_label = image.header["OBJECT"].replace(" ", "")
-    image.target_guess = [parameters.CCD_IMSIZE - float(image.header["OBJECTY"]),
-                          parameters.CCD_IMSIZE - float(image.header["OBJECTX"])]
-    image.disperser_label = image.header["GRATING"]
+    if "OBJECTX" in image.header:
+        image.target_guess = [parameters.CCD_IMSIZE - float(image.header["OBJECTY"]),
+                              parameters.CCD_IMSIZE - float(image.header["OBJECTX"])]
     image.disperser_label = image.header["GRATING"]
-    parameters.DISTANCE2CCD = 116 + float(image.header["LINSPOS"])  # mm
+    parameters.DISTANCE2CCD = 115 + float(image.header["LINSPOS"])  # mm
     image.compute_parallactic_angle()
 
 
@@ -705,13 +735,14 @@ def find_target(image, guess=None, rotated=False, use_wcs=True, widths=[paramete
             theX = x0 - Dx + sub_image_x0
             theY = y0 - Dy + sub_image_y0
             # crop for next iteration
-            if i < niter-1:
+            if i < niter - 1:
                 Dx = Dx // (i + 2)
                 Dy = Dy // (i + 2)
                 x0 = int(theX)
                 y0 = int(theY)
                 NY, NX = sub_image_subtracted.shape
-                sub_image_subtracted = sub_image_subtracted[max(0, int(sub_image_y0) - Dy):min(NY, int(sub_image_y0) + Dy),
+                sub_image_subtracted = sub_image_subtracted[
+                                       max(0, int(sub_image_y0) - Dy):min(NY, int(sub_image_y0) + Dy),
                                        max(0, int(sub_image_x0) - Dx):min(NX, int(sub_image_x0) + Dx)]
                 sub_errors = sub_errors[max(0, int(sub_image_y0) - Dy):min(NY, int(sub_image_y0) + Dy),
                              max(0, int(sub_image_x0) - Dx):min(NX, int(sub_image_x0) + Dx)]

spectractor/extractor/psf.py

@@ -704,33 +704,24 @@ def evaluate(self, pixels, p=None):
         --------
         >>> from spectractor.extractor.images import Image, find_target
         >>> im = Image('tests/data/reduc_20170605_028.fits', target_label="PNG321.0+3.9")
-        >>> im = Image("/Users/jneveu/20200207/10_CCD1_20200207115122.fz", config="config/lpnhe.ini", target_label="HG-AR")
         >>> im.plot_image()
         >>> guess = [820, 580]
-        >>> guess = [490,360]
         >>> parameters.VERBOSE = True
         >>> parameters.DEBUG = True
         >>> x0, y0 = find_target(im, guess)
 
         >>> p = [1,40,50,1,1e20]
         >>> psf = Order0(target=im.target, p=p)
+
+        2D evaluation:
+
         >>> yy, xx = np.mgrid[:80, :100]
         >>> out = psf.evaluate(pixels=np.array([xx, yy]))
 
-        >>> fig = plt.figure(figsize=(5,5))
-        >>> plt.imshow(out, origin="lower")
-        >>> plt.xlabel("X [pixels]")
-        >>> plt.ylabel("Y [pixels]")
-        >>> plt.show()
+        1D evaluation:
 
         >>> out = psf.evaluate(pixels=np.arange(100))
 
-        >>> fig = plt.figure(figsize=(5,5))
-        >>> plt.plot(out)
-        >>> plt.xlabel("X [pixels]")
-        >>> plt.ylabel("Y [pixels]")
-        >>> plt.show()
-
         .. plot::
 
             import matplotlib.pyplot as plt
@@ -743,7 +734,7 @@ def evaluate(self, pixels, p=None):
             parameters.VERBOSE = True
             parameters.DEBUG = True
             x0, y0 = find_target(im, guess)
-            p = [2,40,30,1,1e20]
+            p = [1,40,50,1,1e20]
             psf = Order0(target=im.target, p=p)
             yy, xx = np.mgrid[:80, :100]
             out = psf.evaluate(pixels=np.array([xx, yy]))