Bright star guiding (#371)

* added processor that filters a catalog (not the image) for central stars

* added an image processor for guiding based on the shift of the brightest star

* automatically correct the radius for the binning

* - caluclate altaz-offsets here
- consider magnification in guiding mode (hard coded)

* added option to set the center of the circular mask manually given by the pixel position

* added logging output and cleaned up _get_radec_ref_target

* corrected docstring of center parameter

* removed accidental division of alt-az offsets by 2

* added a test for the circular catalog mask

pyobs/images/processors/misc/__init__.py

@@ -12,3 +12,4 @@
 from .broadcast import Broadcast
 from .circular_mask import CircularMask
 from .image_source_filter import ImageSourceFilter
+from .catalog_circular_mask import CatalogCircularMask

pyobs/images/processors/misc/catalog_circular_mask.py

@@ -0,0 +1,71 @@
+import logging
+from typing import Any, Tuple, Union, Type
+
+import numpy as np
+
+from pyobs.images.processor import ImageProcessor
+from pyobs.images import Image
+
+log = logging.getLogger(__name__)
+
+
+class CatalogCircularMask(ImageProcessor):
+    """Mask catalog for central circle with given radius."""
+
+    __module__ = "pyobs.images.processors.misc"
+
+    def __init__(
+        self,
+        radius: float,
+        center: Union[Tuple[int, int], Tuple[float, float], Tuple[str, str]] = ("CRPIX1", "CRPIX2"),
+        **kwargs: Any
+    ):
+        """Init an image processor that masks out everything except for a central circle.
+
+        Args:
+            radius: radius of the central circle in pixels
+            center: fits-header keywords or pixel coordinates defining the center of the circle
+        """
+        ImageProcessor.__init__(self, **kwargs)
+
+        # init
+        self._radius = radius
+        self._radius_is_corrected = False
+        self._center = center
+
+    async def __call__(self, image: Image) -> Image:
+        """Remove everything outside the given radius from the image.
+
+        Args:
+            image: Image to mask.
+
+        Returns:
+            Image with masked Catalog.
+        """
+        if not self._radius_is_corrected:
+            self._correct_radius_for_binning(image)
+
+        center_x, center_y = self._get_center(image)
+
+        catalog = image.safe_catalog
+        mask = (catalog["x"] - center_x) ** 2 + (catalog["y"] - center_y) ** 2 <= self._radius**2
+        image.catalog = catalog[mask]
+
+        return image
+
+    def _correct_radius_for_binning(self, image):
+        binning_x, binning_y = image.header["XBINNING"], image.header["YBINNING"]
+        if binning_x == binning_y:
+            self._radius /= binning_x
+        else:
+            log.warning("Binning factor is not the same for x and y axis. Filter radius remains uncorrected ...")
+        self._radius_is_corrected = True
+
+    def _get_center(self, image):
+        if isinstance(self._center[0], str) and isinstance(self._center[1], str):
+            return image.header[self._center[0]], image.header[self._center[1]]
+        else:
+            return self._center
+
+
+__all__ = ["CatalogCircularMask"]

pyobs/images/processors/offsets/__init__.py

@@ -2,6 +2,7 @@
 Offsets
 -------
 """
+from .brighteststar_guiding import BrightestStarGuiding
 from .offsets import Offsets
 from .astrometry import AstrometryOffsets
 from .brighteststar import BrightestStarOffsets
@@ -18,6 +19,7 @@
     "ProjectedOffsets",
     "FitsHeaderOffsets",
     "BrightestStarOffsets",
+    "BrightestStarGuiding",
     "DummyOffsets",
-    "DummySkyOffsets"
+    "DummySkyOffsets",
 ]

pyobs/images/processors/offsets/brighteststar_guiding.py

@@ -0,0 +1,94 @@
+import logging
+from typing import Tuple, Any, Optional
+
+from astropy.coordinates import Angle, AltAz, EarthLocation
+from astropy.table import Table, Row
+from astropy.wcs import WCS
+from pandas._typing import npt
+from pyobs.utils.time import Time
+import astropy.units as u
+
+from pyobs.images import Image
+from pyobs.images.meta import PixelOffsets, OnSkyDistance, AltAzOffsets
+from .offsets import Offsets
+
+log = logging.getLogger(__name__)
+
+
+class BrightestStarGuiding(Offsets):
+    """Calculates offsets from the center of the image to the brightest star."""
+
+    __module__ = "pyobs.images.processors.offsets"
+
+    def __init__(self, center_header_cards: Tuple[str, str] = ("CRPIX1", "CRPIX2"), **kwargs: Any):
+        """Initializes a new auto guiding system."""
+        Offsets.__init__(self, **kwargs)
+
+        self._center_header_cards = center_header_cards
+        self._ref_image: Optional[Tuple[npt.NDArray[float], npt.NDArray[float]]] = None
+        self._ref_pos: Optional[Tuple[float, float]] = None
+
+    async def __call__(self, image: Image) -> Image:
+        """Processes an image and sets x/y pixel offset to reference in offset attribute.
+
+        Args:
+            image: Image to process.
+
+        Returns:
+            Original image.
+
+        Raises:
+            ValueError: If offset could not be found.
+        """
+
+        catalog = image.safe_catalog
+        if catalog is None or len(catalog) < 1:
+            log.warning("No catalog found in image.")
+            return image
+
+        if not self._reference_initialized():
+            log.info("Initialising auto-guiding with new image...")
+            self._ref_image = image
+            self._ref_pos = self._get_brightest_star_position(catalog)
+            return image
+
+        star_pos = self._get_brightest_star_position(catalog)
+
+        offset = (star_pos[0] - self._ref_pos[0], star_pos[1] - self._ref_pos[1])
+        image.set_meta(PixelOffsets(*offset))
+        log.info("Found pixel offset of dx=%.2f, dy=%.2f", offset[0], offset[1])
+
+        altaz_offset = self._calc_altaz_offset(image, star_pos)
+        image.set_meta(AltAzOffsets(*altaz_offset))
+        return image
+
+    def _reference_initialized(self):
+        return self._ref_image is not None
+
+    @staticmethod
+    def _get_brightest_star_position(catalog: Table) -> Tuple[float, float]:
+        brightest_star: Row = max(catalog, key=lambda row: row["flux"])
+        log.info("Found brightest star at x=%.2f, y=%.2f", brightest_star["x"], brightest_star["y"])
+        return brightest_star["x"], brightest_star["y"]
+
+    def _calc_altaz_offset(self, image: Image, star_pos: Tuple[float, float]):
+        radec_ref, radec_target = self._get_radec_ref_target(image, star_pos)
+        hdr = image.header
+        location = EarthLocation(lat=hdr["LATITUDE"] * u.deg, lon=hdr["LONGITUD"] * u.deg, height=hdr["HEIGHT"] * u.m)
+        frame = AltAz(obstime=Time(image.header["DATE-OBS"]), location=location)
+        altaz_ref = radec_ref.transform_to(frame)
+        altaz_target = radec_target.transform_to(frame)
+
+        # get offset
+        daz, dalt = altaz_ref.spherical_offsets_to(altaz_target)
+
+        return dalt.arcsec, daz.arcsec
+
+    def _get_radec_ref_target(self, image: Image, star_pos):
+        wcs = WCS(image.header)
+        ref = wcs.pixel_to_world(*self._ref_pos)
+        target = wcs.pixel_to_world(*star_pos)
+        return ref, target
+
+
+__all__ = ["BrightestStarGuiding"]

pyobs/utils/offsets/applyaltazoffsets.py

@@ -1,7 +1,7 @@
 import logging
 from typing import Any
 import numpy as np
-from astropy.coordinates import EarthLocation, AltAz
+from astropy.coordinates import EarthLocation, AltAz, Angle
 import astropy.units as u
 
 from pyobs.images import Image
@@ -53,8 +53,8 @@ async def __call__(self, image: Image, telescope: ITelescope, location: EarthLoc
         if image.has_meta(AltAzOffsets):
             # offsets are Alt/Az, so get them directly
             offsets = image.get_meta(AltAzOffsets)
-            log.info("Found Alt/Az shift of dra=%.2f, ddec=%.2f.", offsets.dalt, offsets.daz)
-            dalt, daz = offsets.dalt * u.arcsec, offsets.daz * u.arcsec
+            log.info("Found Alt/Az shift of dalt=%.2f, daz=%.2f.", offsets.dalt, offsets.daz)
+            dalt, daz = Angle(offsets.dalt * u.arcsec), Angle(offsets.daz * u.arcsec)
 
         elif image.has_meta(RaDecOffsets):
             raise NotImplementedError()

tests/images/processors/misc/test_catalog_circular_mask.py

@@ -0,0 +1,22 @@
+import numpy as np
+import pytest
+
+from astropy.table import Table
+from pyobs.images import Image
+from pyobs.images.processors.misc import CatalogCircularMask
+
+
+@pytest.mark.asyncio
+async def test_call():
+    radius = 1
+    mask = CatalogCircularMask(radius=radius)
+    data = np.ones((4, 4))
+    catalog = Table([[1, 3], [1, 3]], names=("x", "y"))
+    image = Image(data, catalog=catalog)
+    image.header["CRPIX1"] = 1.5
+    image.header["CRPIX2"] = 1.5
+    image.header["XBINNING"] = 1
+    image.header["YBINNING"] = 1
+    masked_image = await mask(image)
+    expected_output_catalog = Table([[1], [1]], names=("x", "y"))
+    assert masked_image.catalog == expected_output_catalog