From 707407de21876d5f2d0ff6d381c5babf5db834d6 Mon Sep 17 00:00:00 2001 From: LeMee <95364870+LeonMee@users.noreply.github.com> Date: Wed, 2 Oct 2024 13:35:16 +0200 Subject: [PATCH] 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 | 1 + .../processors/misc/catalog_circular_mask.py | 71 ++++++++++++++ pyobs/images/processors/offsets/__init__.py | 4 +- .../offsets/brighteststar_guiding.py | 94 +++++++++++++++++++ pyobs/utils/offsets/applyaltazoffsets.py | 6 +- .../misc/test_catalog_circular_mask.py | 22 +++++ 6 files changed, 194 insertions(+), 4 deletions(-) create mode 100644 pyobs/images/processors/misc/catalog_circular_mask.py create mode 100644 pyobs/images/processors/offsets/brighteststar_guiding.py create mode 100644 tests/images/processors/misc/test_catalog_circular_mask.py diff --git a/pyobs/images/processors/misc/__init__.py b/pyobs/images/processors/misc/__init__.py index af6b02dc..de4dc6b7 100644 --- a/pyobs/images/processors/misc/__init__.py +++ b/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 diff --git a/pyobs/images/processors/misc/catalog_circular_mask.py b/pyobs/images/processors/misc/catalog_circular_mask.py new file mode 100644 index 00000000..6f4dfa38 --- /dev/null +++ b/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"] diff --git a/pyobs/images/processors/offsets/__init__.py b/pyobs/images/processors/offsets/__init__.py index edff22f2..f4af6f27 100644 --- a/pyobs/images/processors/offsets/__init__.py +++ b/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", ] diff --git a/pyobs/images/processors/offsets/brighteststar_guiding.py b/pyobs/images/processors/offsets/brighteststar_guiding.py new file mode 100644 index 00000000..2143cbb1 --- /dev/null +++ b/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"] diff --git a/pyobs/utils/offsets/applyaltazoffsets.py b/pyobs/utils/offsets/applyaltazoffsets.py index 2306ab20..b389c5f4 100644 --- a/pyobs/utils/offsets/applyaltazoffsets.py +++ b/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() diff --git a/tests/images/processors/misc/test_catalog_circular_mask.py b/tests/images/processors/misc/test_catalog_circular_mask.py new file mode 100644 index 00000000..ec40fc96 --- /dev/null +++ b/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