Merge pull request #28 from gnthibault/develop

Added offset adjustment specifically for slit spectroscopy usecase

Camera/AbstractCamera.py

@@ -27,7 +27,9 @@ def __init__(self, serv_time, *args,  **kwargs):
 
         self.serv_time = serv_time
         self.do_acquisition = kwargs.get("do_acquisition", False)
+        self.do_guiding = kwargs.get("do_guiding", False)
         self.do_pointing = kwargs.get("do_pointing", False)
+        self.do_adjust_pointing = kwargs.get("do_adjust_pointing", False)
         self.do_autofocus = kwargs.get("do_autofocus", False)
 
         self.camera_name = kwargs["camera_name"]

Camera/IndiAbstractCamera.py

@@ -22,26 +22,31 @@ def __init__(self, serv_time, config=None, connect_on_create=True):
                            connect_on_create=connect_on_create)
         self.indi_camera_config = config
 
+
     # TODO TN: setup event based acquisition properly
     def shoot_asyncWithEvent(self, exp_time_sec, filename, exposure_event,
                              **kwargs):
-        # set frame type
-        frame_type = kwargs.get("frame_type", "FRAME_LIGHT")
-        self.set_frame_type(frame_type)
-        # set gain
-        gain = kwargs.get("gain", self.gain)
-        self.set_gain(gain)
-        # set temperature
-        temperature = kwargs.get("temperature", None)
-        if temperature is not None:
-            self.set_cooling_on()
-            self.set_temperature(temperature)
-        # Now shoot
-        self.setExpTimeSec(exp_time_sec)
-        self.logger.debug(f"Camera {self.camera_name}, about to shoot for {self.exp_time_sec}")
-        self.shoot_async()
+        # If there is no external trigger, then we proceed to handle setup on our side
+        external_trigger = kwargs.get("external_trigger", False)
+        if not external_trigger:
+            # set frame type
+            frame_type = kwargs.get("frame_type", "FRAME_LIGHT")
+            self.set_frame_type(frame_type)
+            # set gain
+            gain = kwargs.get("gain", self.gain)
+            self.set_gain(gain)
+            # set temperature
+            temperature = kwargs.get("temperature", None)
+            if temperature is not None:
+                self.set_cooling_on()
+                self.set_temperature(temperature)
+            # Now shoot
+            self.setExpTimeSec(exp_time_sec)
+            self.logger.debug(f"Camera {self.camera_name}, about to shoot for {self.exp_time_sec}")
+            self.shoot_async()
+        # Wether trigger was internal or external, we rely on the last received blob
         self.synchronize_with_image_reception()
-        self.logger.debug(f"Camera {self.camera_name}, done with image reception")
+        self.logger.debug(f"Camera {self.camera_name}, done with image reception, external trigger {external_trigger}")
         image = self.get_received_image()
         try:
             with open(filename, "wb") as f:

Camera/IndiCamera.py

@@ -37,8 +37,12 @@ def __init__(self, logger=None, config=None, connect_on_create=True):
         if config is None:
             config = dict(
                 camera_name='CCD Simulator',
-                autofocus_seconds=5,
                 pointing_seconds=30,
+                adjust_center_x=400,
+                adjust_center_y= 400,
+                adjust_roi_search_size=50,
+                adjust_pointing_seconds=5,
+                autofocus_seconds=5,
                 autofocus_roi_size=500,
                 autofocus_merit_function="half_flux_radius",
                 focuser=dict(
@@ -72,18 +76,19 @@ def __init__(self, logger=None, config=None, connect_on_create=True):
             self.connect()
 
         # Specific initialization
-        self.autofocus_seconds = float(config['autofocus_seconds'])
         self.pointing_seconds = float(config['pointing_seconds'])
+        self.adjust_center_x = float(config['adjust_center_x'])
+        self.adjust_center_y = float(config['adjust_center_y'])
+        self.adjust_roi_search_size = int(config['adjust_roi_search_size'])
+        self.adjust_pointing_seconds = float(config['adjust_pointing_seconds'])
+        self.autofocus_seconds = float(config['autofocus_seconds'])
         self.autofocus_roi_size = int(config['autofocus_roi_size'])
         self.autofocus_merit_function = config['autofocus_merit_function']
         self._setup_focuser(config, connect_on_create)
         self._setup_filter_wheel(config, connect_on_create)
 
         self.logger.debug(f"Indi camera, camera name is: {device_name}")
 
-        # Frame Blob: reference that will be used to receive binary
-        self.last_blob = None
-
         # Default exposureTime, gain
         self.exp_time_sec = 5
         self.gain = 400
@@ -141,8 +146,6 @@ def prepare_shoot(self):
         self.logger.debug('Indi client will register to server in order to '
                           'receive blob CCD1 when it is ready')
         self.indi_client.enable_blob()
-        #self.indi_client.setBLOBMode(PyIndi.B_ALSO, self.device_name, 'CCD1')
-        #self.frame_blob = self.get_prop(propName='CCD1', propType='blob')
 
     def synchronize_with_image_reception(self):
         try:
@@ -170,8 +173,6 @@ def shoot_async(self):
             self.set_number('CCD_EXPOSURE',
                             {'CCD_EXPOSURE_VALUE': self.sanitize_exp_time(self.exp_time_sec)},
                             sync=False)
-            #     self.last_blob = blob_listener.get(timeout=self.exp_time_sec +
-            #                                        self.READOUT_TIME_MARGIN)
         except Exception as e:
             self.logger.error(f"Indi Camera Error in shoot: {e}")
 
@@ -294,7 +295,7 @@ def get_gain(self):
         return gain["GAIN"]
 
     def get_frame_type(self):
-        return self.get_prop('CCD_FRAME_TYPE', 'switch')
+        return self.get_switch('CCD_FRAME_TYPE')
 
     def set_frame_type(self, frame_type):
         """

FilterWheel/IndiFilterWheel.py

@@ -66,12 +66,12 @@ def set_filter_number(self, number):
         self.set_number('FILTER_SLOT', {'FILTER_SLOT_VALUE': number})
 
     def currentFilter(self):
-        ctl = self.get_prop('FILTER_SLOT', 'number')
+        ctl = self.get_number('FILTER_SLOT')
         number = int(ctl[0].value)
         return number, self.filterName(number)
 
     def filters(self):
-        ctl = self.get_prop('FILTER_NAME', 'text')
+        ctl = self.get_text('FILTER_NAME')
         filters = [(x.text, IndiFilterWheel.__name2number(x.name)) for x in ctl]
         return dict(filters)
 

Guider/GuiderPHD2.py

@@ -20,7 +20,7 @@
 
 MAXIMUM_CALIBRATION_TIMEOUT = 6 * 60 * u.second
 FIND_STAR_TIMEOUT           = 60 * u.second
-POSITION_LOCKING_TIMEOUT    = 60 * u.second
+POSITION_LOCKING_TIMEOUT    = 5 * 60 * u.second
 MAXIMUM_DITHER_TIMEOUT      = 45 * u.second
 MAXIMUM_PAUSING_TIMEOUT     = 30 * u.second
 STANDARD_TIMEOUT            = 120 * u.second
@@ -607,6 +607,11 @@ def get_pixel_scale(self):
             self.logger.warning(msg)
             raise RuntimeError(msg)
 
+    def is_lock_position_close_to(self, px_target, max_angle_sep):
+        max_pixel_sep = self.convert_angle_sep_to_pixels(max_angle_sep)
+        lock_pos = self.get_lock_position()
+        return np.linalg.norm(np.array(px_target)-np.array(lock_pos)) < max_pixel_sep
+
     def convert_angle_sep_to_pixels(self, angle_sep):
         pixel_scale = self.get_pixel_scale()
         return angle_sep.to(u.arcsec).value * pixel_scale

Imaging/AutoFocuser.py

@@ -207,7 +207,7 @@ def autofocus(self,
         assert self.camera.focuser.is_connected, self.logger.error(
             f"Focuser {self.camera.focuser} must be connected for autofocus")
 
-        if not autofocus_status:
+        if autofocus_status is None:
             autofocus_status = [False]
 
         if not focus_range:
@@ -335,8 +335,10 @@ def _autofocus(self,
             assert np.isfinite(final_focus)
         except Exception as e:
             self.logger.error(f"Focusing method failed: {e}")
+            autofocus_status[0] = False
             if focus_event is not None:
                 focus_event.set()
+            return
         autofocus_status[0] = True
         if focus_event is not None:
             focus_event.set()

Imaging/SolvedImageAnalysis.py

@@ -0,0 +1,149 @@
+# Generic tools
+import os
+
+# numerical tools
+import numpy as np
+
+# Data viz
+import matplotlib.pyplot as plt
+
+# Astropy
+from astropy.io import fits
+from astropy import units as u
+from astropy.visualization import AsymmetricPercentileInterval, ImageNormalize, MinMaxInterval, SqrtStretch
+from astropy.wcs import WCS
+
+def find_best_candidate_star(pointing_image,
+                             target,
+                             max_identification_error,
+                             make_plots=True):
+    """
+    It is very important to read this page to understand how to preoperly use wcs:
+    https://docs.astropy.org/en/stable/wcs/note_sip.html
+
+    This has not been a strict convention in the past and the default in astropy.wcs is to always include the SIP
+    distortion if the SIP coefficients are present, even if -SIP is not included in CTYPE. The presence of a -SIP
+    suffix in CTYPE is not used as a trigger to initialize the SIP distortion.
+    It is important that headers implement correctly the SIP convention. If the intention is to use the SIP distortion,
+    a header should have the SIP coefficients and the -SIP suffix in CTYPE.
+    astropy.wcs prints INFO messages when inconsistent headers are detected, for example when SIP coefficients are
+    present but CTYPE is missing a -SIP suffix, see examples below. astropy.wcs will print a message about the
+    inconsistent header but will create and use the SIP distortion and it will be used in calls to all_pix2world.
+    If this was not the intended use (e.g. it’s a drizzled image and has no distortions) it is best to remove the SIP
+    coefficients from the header. They can be removed temporarily from a WCS object by:
+    >>> wcsobj.sip = None
+    In addition, if SIP is the only distortion in the header, the two methods, wcs_pix2world and wcs_world2pix,
+    may be used to transform from pixels to world coordinate system while omitting distortions. Another consequence
+    of the inconsistent header is that if to_header() is called with relax=True it will return a header with SIP
+    coefficients and a -SIP suffix in CTYPE and will not reproduce the original header.
+    In conclusion, when astropy.wcs detects inconsistent headers, the recommendation is that the header is inspected
+    and corrected to match the data. Below is an example of a header with SIP coefficients when -SIP is missing from
+    CTYPE. The data is drizzled, i.e. distortion free, so the intention is not to include the SIP distortion.
+
+    Now from: https://docs.astropy.org/en/stable/api/astropy.wcs.WCS.html#astropy.wcs.WCS.all_pix2world
+    Transforms pixel coordinates to world coordinates. Performs all of the following in series:
+        - Detector to image plane correction (if present in the FITS file)
+        - SIP distortion correction (if present in the FITS file)
+        - distortion paper table-lookup correction (if present in the FITS file)
+        - wcslib “core” WCS transformation
+
+    :param pointing_image:
+    :param target:
+    :param max_identification_error:
+    :param make_plots:
+    :return:
+    """
+
+    img_directory = os.path.dirname(pointing_image.fits_file)
+    img_filename = pointing_image.fits_file
+    xy_filename = os.path.splitext(img_filename)[0]+".axy"
+
+    # Open image file
+    hdu = fits.open(img_filename)[0]
+    wcs = WCS(hdu.header)
+    data = hdu.data.astype(np.float32)
+
+    # Open star detection image file
+    detections_data = fits.open(xy_filename)[1].data
+
+    # detection data
+    detections = [((x,y), flux) for x,y,flux,bkg in detections_data]
+    #fluxes = np.array([flux for (x, y), flux in detections])
+    px_centers = np.array([[x, y] for (x, y), flux in detections])            # np.array of shape is (n, 2)
+    rd_centers = [wcs.pixel_to_world(*px.tolist()) for px in px_centers] # list of SkyCoord
+    # rd_centers = list(map(
+    #     lambda x: SkyCoord(x[0]*u.hourangle, x[1]*u.degree),
+    #     wcs.all_pix2world(px_centers, 0).tolist()))
+    # site-packages/astropy/wcs/wcsapi/fitswcs.py:347: UserWarning: 'WCS.all_world2pix' failed to converge to the requested accuracy.
+    # After 20 iterations, the solution is diverging at least for one input point.
+
+    rd_target_star_reference = target
+    px_target_star_reference = wcs.world_to_pixel(rd_target_star_reference)
+
+    i_closest = min(range(len(rd_centers)), key=lambda i: rd_centers[i].separation(rd_target_star_reference))
+    rd_closest_detection = rd_centers[i_closest]
+    closest_distance = rd_closest_detection.separation(rd_target_star_reference)
+    if closest_distance <= max_identification_error:
+        px_candidate_star = px_centers[i_closest] #wcs.world_to_pixel(rd_closest_detection)
+    else:
+        px_candidate_star = None
+        rd_closest_detection = None
+
+    if make_plots:
+        ax = plt.subplot(projection=wcs, label='overlays')
+        fig = ax.get_figure()
+        norm = ImageNormalize(data, interval=MinMaxInterval(), stretch=SqrtStretch())
+        ax.imshow(data, origin='lower', cmap='gray', norm=norm)
+
+        # Show detected star in green, or theoretical star in red
+        radius = max(3, int(np.ceil(min(data.shape)*0.005)))
+        if px_candidate_star is not None:
+            circle = plt.Circle(px_candidate_star, radius, color='g', fill=False)
+        else:
+            circle = plt.Circle(px_target_star_reference, radius, color='r', fill=False)
+        ax.add_patch(circle)
+
+        ra = ax.coords['ra']
+        ra.set_ticks()
+        #ra.set_ticks_position()
+        ra.set_ticklabel(size=12)
+        ra.set_ticklabel_position('l')
+        ra.set_axislabel('RA', minpad=0.3)
+        ra.set_axislabel_position('l')
+        ra.grid(color='yellow', ls='-', alpha=0.3)
+        ra.set_format_unit(u.hourangle)
+        ra.set_major_formatter('hh:mm:ss')
+
+        dec = ax.coords['dec']
+        dec.set_ticks()
+        #dec.set_ticks_position()
+        dec.set_ticklabel(size=12)
+        dec.set_ticklabel_position('b')
+        dec.set_axislabel('DEC', minpad=0.3)
+        dec.set_axislabel_position('b')
+        dec.grid(color='yellow', ls='-', alpha=0.3)
+        dec.set_format_unit(u.degree)
+        dec.set_major_formatter('dd:mm:ss')
+
+        fig.set_tight_layout(True)
+        plot_path = f"{img_directory}/pointing_detection.jpg"
+        fig.savefig(plot_path, transparent=False)
+
+        # explicitly close and delete figure
+        fig.clf()
+        del fig
+    return px_candidate_star
+
+def get_brightest_detection(pointing_image):
+    img_filename = pointing_image.fits_file
+    xy_filename = os.path.splitext(img_filename)[0]+".axy"
+
+    # Open star detection image file
+    detections_data = fits.open(xy_filename)[1].data
+
+    # detection data
+    detections = [((x, y), flux) for x,y,flux,bkg in detections_data]
+    fluxes = np.array([flux for (x, y), flux in detections])
+    px_centers = np.array([[x, y] for (x, y), flux in detections])  # np.array of shape is (n, 2)
+    max_flux_index = np.argmax(fluxes)
+    return px_centers[max_flux_index]