Merge branch 'develop'

Guider/GuiderPHD2.py

@@ -467,11 +467,12 @@ def find_star(self, x=None, y=None, width=None, height=None):
                    otherwise returns an error object
            desc. : Select a star
         """
-        params = [int(x), int(y), int(width), int(height)]
+        params = [x, y, width, height]
         if all(map(lambda x: x is None, params)):
             params = []
         else:
             assert not any(map(lambda x: x is None, params)), f"Please provide all parameters, invalid {params}"
+            params = [int(x), int(y), int(width), int(height)]
             params = [params]
         req={"method": "find_star",
              "params": params,

Imaging/Image.py

@@ -1,12 +1,14 @@
 # generic import
 from collections import namedtuple
+import datetime
 import os
 
 # Astropy
 from astropy import units as u
 from astropy import wcs
 from astropy.coordinates import EarthLocation
 from astropy.coordinates import FK5
+from astropy.coordinates import ICRS
 from astropy.coordinates import SkyCoord
 from astropy.io import fits
 from astropy.time import Time
@@ -127,7 +129,7 @@ def pointing_error(self, pointing_reference_coord=None):
         Returns:
             namedtuple: Pointing error information
         """
-        if self._pointing_error is None:
+        if (self._pointing_error is None) or not (pointing_reference_coord is None):
             assert self.pointing is not None, "No world coordinate system (WCS), can't get pointing_error"
             assert self.header_pointing is not None
 
@@ -153,7 +155,9 @@ def pointing_error(self, pointing_reference_coord=None):
 
     def get_center_coordinates(self):
         """
-        Those are 0-based indexing
+        Those are 0-based indexing, also from documentation:
+        https://docs.astropy.org/en/stable/api/astropy.wcs.WCS.html#astropy.wcs.WCS.pixel_to_world_values
+        Note that pixel coordinates are assumed to be 0 at the center of the first pixel in each dimension
         :return:
         """
         return (self.header["NAXIS1"]/2-0.5,
@@ -164,7 +168,7 @@ def get_header_pointing(self):
 
         The header should contain the `RA-FIELD` and `DEC-FIELD` keywords, from
         which the header pointing coordinates are built. Those two entries are
-        written by us as J2000 format, which is also the expected format for the
+        written by us as J2000 format, which might not be the expected format for the
         astrometric resolution
         """
         try:
@@ -183,18 +187,19 @@ def get_header_pointing(self):
             msg = 'Cannot get header pointing information: {}'.format(e)
             self.logger.error(msg)
             raise RuntimeError(msg)
-
 
     def get_wcs_pointing(self):
         """Get the pointing information from the WCS
         Builds the pointing coordinates from the plate-solved WCS. These will be
         compared with the coordinates stored in the header.
-        One should notice that Astrometry.net uses J2000 equinox
+        One should notice that Astrometry.net uses J2000 equinox in general,
+        however, we assume the install related to Gaia DR2 + Tycho2 catalogs, which are J2015.5
+        see https://portal.nersc.gov/project/cosmo/temp/dstn/index-5200/LITE/
         """
         if self.wcs is not None:
             # This was wrong because crval coord relates to crpix reference pixels
-            ra = self.wcs.celestial.wcs.crval[0]
-            dec = self.wcs.celestial.wcs.crval[1]
+            # ra = self.wcs.celestial.wcs.crval[0]
+            # dec = self.wcs.celestial.wcs.crval[1]
 
             # TODO TN trying alternative definition
             cx, cy = self.get_center_coordinates()
@@ -204,16 +209,45 @@ def get_wcs_pointing(self):
                 0,  # 0-based indexing
                 ra_dec_order=True)
 
+            # fk5_J2015_5 = FK5(equinox=Time(datetime.datetime(year=2015, month=6, day=30)))
+            # pointing = SkyCoord(ra=radeg*u.degree,
+            #                     dec=decdeg*u.degree,
+            #                     frame=fk5_J2015_5,
+            #                     obstime=Time(datetime.datetime(year=2015, month=6, day=30)))
+            # icrs_j2k = ICRS()
+            # self.pointing = pointing.transform_to(icrs_j2k)
+
+            icrs_j2k = ICRS()
             self.pointing = SkyCoord(ra=radeg*u.degree,
-                                     dec=decdeg*u.degree,
-                                     frame='icrs', equinox='J2000.0')
+                                dec=decdeg*u.degree,
+                                frame=icrs_j2k)
+
             self.ra = self.pointing.ra.to(u.hourangle)
             self.dec = self.pointing.dec.to(u.degree)
             # Precess to the current equinox otherwise the RA - LST method
             # will be off.
             #self.ha = self.pointing.transform_to(self.FK5_Jnow).ra.to(
             #    u.hourangle) - self.sidereal
 
+    def all_pix2world(self, x, y):
+        coord_ra, coord_dec = self.wcs.all_pix2world(
+            x,
+            y,
+            0,  # 0-based indexing
+            ra_dec_order=True)
+        # fk5_J2015_5 = FK5(equinox=Time(datetime.datetime(year=2015, month=6, day=30)))
+        # coord = SkyCoord(ra=coord_ra*u.deg,
+        #                  dec=coord_dec*u.deg,
+        #                  frame=fk5_J2015_5,
+        #                  obstime=Time(datetime.datetime(year=2015, month=6, day=30)))
+        # icrs_j2k = ICRS()
+        icrs_j2k = ICRS()
+        coord = SkyCoord(ra=coord_ra * u.deg,
+                                 dec=coord_dec * u.deg,
+                                 frame=icrs_j2k)
+
+        return coord.transform_to(icrs_j2k)
+
     def solve_field(self, **kwargs):
         """ Solve field and populate WCS information
             If you use basic catalog for astrometry.net, it is J2K!

Imaging/fits.py

@@ -9,9 +9,10 @@
 import skimage.io as io
 
 # Astropy
+from astropy import units as u
 from astropy.io import fits
+from astropy.time import Time
 from astropy.wcs import WCS
-from astropy import units as u
 
 # Local
 from utils import error

Mount/IndiMount.py

@@ -304,7 +304,8 @@ def get_current_coordinates(self):
                        frame=fk5_now,
                        obstime=Time.now())
         icrs_j2k = ICRS()
-        self.logger.debug(f"Received coordinates in JNOw/CIRS from mount: {ret}")
+        # This was too verbose
+        #self.logger.debug(f"Received coordinates in JNOw/CIRS from mount: {ret}")
         ret = ret.transform_to(icrs_j2k)
 
         # ret = SkyCoord(ra=rahour_decdeg['RA'] * u.hourangle,

Pointer/DifferentialPointer.py

@@ -101,10 +101,10 @@ def points_async(self, mount, camera, observation, fits_headers, pointing_event,
                 self.logger.debug(f"Pointing Coords: {pointing_image.pointing}")
                 self.logger.debug(f"Pointing Error: {pointing_error}")
                 # adjust by slewing again to correct the delta
-                target = mount.get_current_coordinates()
+                current = mount.get_current_coordinates()
                 target = SkyCoord(
-                    ra=target.ra-pointing_error.delta_ra,
-                    dec=target.dec-pointing_error.delta_dec,
+                    ra=current.ra-pointing_error.delta_ra,
+                    dec=current.dec-pointing_error.delta_dec,
                     frame='icrs', equinox='J2000.0')
                 mount.slew_to_coord(target)
                 # update pointing process tracking information

Pointer/InvisibleObjectOffsetPointer.py

@@ -2,18 +2,17 @@
 import threading
 import traceback
 
-# Astropy
-from astropy.coordinates import SkyCoord
-
-# Local
-from utils.error import AstrometrySolverError, ImageAcquisitionError
+# Numerical tools
+import numpy as np
 
 # Astropy
+from astropy.coordinates import SkyCoord
 import astropy.units as u
 
 # Local
 from Imaging.Image import OffsetError
 from Pointer.OffsetPointer import OffsetPointer
+from utils.error import AstrometrySolverError, ImageAcquisitionError
 
 class InvisibleObjectOffsetPointer(OffsetPointer):
     def __init__(self, config=None):
@@ -81,6 +80,7 @@ def offset_points_async(self, mount, camera, guiding_camera, guider, observation
                 msg = f"Going to adjust pointing, need to stop guiding"
                 self.logger.debug(msg)
                 guider.stop_capture()
+                #guider.set_paused(paused=True, full="full")
                 try:
                     exp_time_sec = guiding_camera.is_remaining_exposure_time()
                     guiding_camera.synchronize_with_image_reception(exp_time_sec=exp_time_sec)
@@ -91,7 +91,7 @@ def offset_points_async(self, mount, camera, guiding_camera, guider, observation
             pointing_error = OffsetError(*(np.inf * u.arcsec,) * 3)
             pointing_error_stack = {}
 
-            while (img_num < self.max_iterations and pointing_error.magnitude > self.max_pointing_error.magnitude):
+            while (img_num <= self.max_iterations and pointing_error.magnitude > self.max_pointing_error.magnitude):
                 pointing_image = self.acquire_pointing(camera, guiding_camera, observation, fits_headers,
                                                        img_num)
                 try:
@@ -105,53 +105,42 @@ def offset_points_async(self, mount, camera, guiding_camera, guider, observation
                     # update mount with the actual position
                     if self.sync_mount_upon_solve:
                         mount.sync_to_coord(pointing_image.pointing)
-                    px_identified_target = pointing_image.get_center_coordinates()
+                    # There are some subteleties here: https://astropy-cjhang.readthedocs.io/en/latest/wcs/
+                    current_sky_coord_of_target_sensor_position = pointing_image.all_pix2world(
+                        camera.adjust_center_x,
+                        camera.adjust_center_y)
                 except AstrometrySolverError as e:
                     msg = f"Cannot solve image {pointing_image.fits_file} while in offset_pointing state: {e}"
                     self.logger.error(msg)
                     raise RuntimeError(msg)
-                    # TODO TN We can decide to assume that main camera is already pointing to the expected object
-
-                # There are some subteleties here: https://astropy-cjhang.readthedocs.io/en/latest/wcs/
-                radeg, decdeg = pointing_image.wcs.all_pix2world(
-                    camera.adjust_center_x,
-                    camera.adjust_center_y,
-                    0,  # 0-based indexing
-                    ra_dec_order=True)
-                current_sky_coord_of_target_sensor_position = SkyCoord(
-                    ra=float(radeg) * u.degree,
-                    dec=float(decdeg) * u.degree,
-                    frame='icrs',
-                    equinox='J2000.0')
-                radeg, decdeg = pointing_image.wcs.all_pix2world(
-                    px_identified_target[0],
-                    px_identified_target[1],
-                    0,  # 0-based indexing
-                    ra_dec_order=True)
-                current_sky_coord_of_target_star = SkyCoord(
-                    ra=float(radeg) * u.degree,
-                    dec=float(decdeg) * u.degree,
-                    frame='icrs',
-                    equinox='J2000.0')
 
-                pointing_error = pointing_image.pointing_error(
-                    pointing_reference_coord=current_sky_coord_of_target_star
-                )
-                offset_delta_ra = current_sky_coord_of_target_sensor_position.ra - current_sky_coord_of_target_star.ra
-                offset_delta_dec = current_sky_coord_of_target_sensor_position.dec - current_sky_coord_of_target_star.dec
+                star_pointing_delta = pointing_image.pointing_error()
+                offset_delta_ra = current_sky_coord_of_target_sensor_position.ra - pointing_image.pointing.ra
+                offset_delta_dec = current_sky_coord_of_target_sensor_position.dec - pointing_image.pointing.dec
                 # adjust by slewing to the opposite of the delta
                 current = mount.get_current_coordinates()
                 target = SkyCoord(
-                    ra=current.ra + pointing_error.delta_ra - offset_delta_ra,
-                    dec=current.dec + pointing_error.delta_dec - offset_delta_dec,
+                    ra=current.ra - star_pointing_delta.delta_ra - offset_delta_ra,
+                    dec=current.dec - star_pointing_delta.delta_dec - offset_delta_dec,
                     frame='icrs', equinox='J2000.0')
-                # Now adjust by slewing to the specified counter-offseted coordinates
-                mount.slew_to_coord(target)
-
+                # Virtual target (different from target if we do not do sync on the mount)
+                virtual_target = SkyCoord(
+                    ra=pointing_image.pointing.ra - star_pointing_delta.delta_ra - offset_delta_ra,
+                    dec=pointing_image.pointing.dec - star_pointing_delta.delta_dec - offset_delta_dec,
+                    frame='icrs', equinox='J2000.0')
+                pointing_error = pointing_image.pointing_error(
+                    pointing_reference_coord=virtual_target
+                )
+                self.logger.debug(f"Offset point, current error is: {pointing_error}")
                 # update pointing process tracking information
                 pointing_error_stack[img_num] = pointing_error
                 img_num = img_num + 1
 
+                if (img_num <= self.max_iterations) and (pointing_error.magnitude >
+                                                         self.max_pointing_error.magnitude):
+                    # Now adjust by slewing to the specified counter-offseted coordinates
+                    mount.slew_to_coord(target)
+
             if pointing_error.magnitude > self.max_pointing_error.magnitude:
                 self.logger.error(f"Pointing accuracy was not good enough after "
                                   f"{img_num} iterations, pointing error stack was: "