Add and use sky <-> det wcs convenience functions and tests

scopesim/optics/fov_manager.py

@@ -49,6 +49,7 @@
 
 import numpy as np
 from astropy import units as u
+from astropy.wcs import WCS
 
 from . import image_plane_utils as ipu
 from ..effects import DetectorList
@@ -121,8 +122,6 @@ def generate_fovs_list(self):
         # ..todo: add catch to split volumes larger than chunk_size
         pixel_scale = from_currsys(self.meta["pixel_scale"])
         plate_scale = from_currsys(self.meta["plate_scale"])
-        pixel_size = pixel_scale / plate_scale
-        plate_scale_deg = plate_scale / 3600.  # ["/mm] / 3600 = [deg/mm]
 
         chunk_size = from_currsys(self.meta["chunk_size"])
         max_seg_size = from_currsys(self.meta["max_segment_size"])
@@ -149,11 +148,9 @@ def generate_fovs_list(self):
                                                 [ys_min, ys_max],
                                                 pixel_scale=pixel_scale / 3600.)
 
-            xy_sky = ipu.calc_footprint(skyhdr)
-            xy_det = xy_sky / plate_scale_deg
-            dethdr = ipu.header_from_list_of_xy(xy_det[:, 0], xy_det[:, 1],
-                                                pixel_size, "D")
-            skyhdr.update(dethdr)
+            dethdr, _ = ipu.det_wcs_from_sky_wcs(
+                WCS(skyhdr), pixel_scale, plate_scale)
+            skyhdr.update(dethdr.to_header())
 
             # useful for spectroscopy mode where slit dimensions is not the same
             # as detector dimensions

scopesim/optics/image_plane_utils.py

@@ -1023,3 +1023,96 @@ def _get_unit_from_headers(*headers, wcs_suffix: str = "") -> str:
         [(i, header[f"CUNIT{i}{wcs_suffix}"])
          for header, i in product(headers, range(1, 3))]
     return unit
+
+
+def det_wcs_from_sky_wcs(sky_wcs: WCS,
+                         pixel_scale: float,
+                         plate_scale: float,
+                         naxis=None) -> Tuple[WCS, np.ndarray]:
+    """
+    Create detector WCS from celestial WCS using pixel and plate scales.
+
+    Parameters
+    ----------
+    sky_wcs : astropy.wcs.WCS
+        Celestial WCS.
+    pixel_scale : float
+        Quantity or float (assumed to be arcsec / pixel).
+    plate_scale : float
+        Quantity or float (assumed to be arcsec / mm).
+    naxis : (int, int), optional
+        Shape of the image, usually ``NAXIS1`` and ``NAXIS2``. If the input WCS
+        holds this information, the default None will use that. Otherwise not
+        providing `naxis` will raise and error.
+
+    Returns
+    -------
+    det_wcs : astropy.wcs.WCS
+        Detector WCS.
+    det_naxis : (int, int)
+        Shape of the image (``NAXIS1``, ``NAXIS2``).
+
+    """
+    # TODO: Using astropy units for now to avoid deg vs. arcsec confusion.
+    #       Once Scopesim is consistent there, remove astropy units.
+    pixel_scale <<= u.arcsec / u.pixel
+    plate_scale <<= u.arcsec / u.mm
+    logging.debug("Pixel scale: %s", pixel_scale)
+    logging.debug("Plate scale: %s", plate_scale)
+
+    pixel_size = pixel_scale / plate_scale
+    # TODO: add check if cunit is consistent along all axes
+    cunit = sky_wcs.wcs.cunit[0]
+    corners = sky_wcs.calc_footprint(center=False, axes=naxis) * cunit
+    logging.debug("WCS sky corners: %s", corners)
+    corners /= plate_scale
+    corners = corners.to(u.mm)
+    logging.debug("WCS det corners: %s", corners)
+
+    return create_wcs_from_points(corners, pixel_size, "D")
+
+
+def sky_wcs_from_det_wcs(det_wcs: WCS,
+                         pixel_scale: float,
+                         plate_scale: float,
+                         naxis=None) -> Tuple[WCS, np.ndarray]:
+    """
+    Create celestial WCS from detector WCS using pixel and plate scales.
+
+    Parameters
+    ----------
+    det_wcs : astropy.wcs.WCS
+        Detector WCS.
+    pixel_scale : float
+        Quantity or float (assumed to be arcsec / pixel).
+    plate_scale : float
+        Quantity or float (assumed to be arcsec / mm).
+    naxis : (int, int), optional
+        Shape of the image, usually ``NAXIS1`` and ``NAXIS2``. If the input WCS
+        holds this information, the default None will use that. Otherwise not
+        providing `naxis` will raise and error.
+
+    Returns
+    -------
+    sky_wcs : astropy.wcs.WCS
+        Celestial WCS.
+    sky_naxis : (int, int)
+        Shape of the image (``NAXIS1``, ``NAXIS2``).
+
+    """
+    # TODO: Using astropy units for now to avoid deg vs. arcsec confusion.
+    #       Once Scopesim is consistent there, remove astropy units.
+    pixel_scale <<= u.arcsec / u.pixel
+    plate_scale <<= u.arcsec / u.mm
+    logging.debug("Pixel scale: %s", pixel_scale)
+    logging.debug("Plate scale: %s", plate_scale)
+
+    # TODO: add check if cunit is consistent along all axes
+    cunit = det_wcs.wcs.cunit[0]
+    corners = det_wcs.calc_footprint(center=False, axes=naxis) * cunit
+    logging.debug("WCS det corners: %s", corners)
+    corners *= plate_scale
+    corners = corners.to(u.arcsec)
+    logging.debug("WCS sky corners: %s", corners)
+
+    return create_wcs_from_points(corners, pixel_scale)

scopesim/tests/tests_optics/test_ImagePlaneUtils.py

@@ -3,6 +3,7 @@
 from copy import deepcopy
 import numpy as np
 from astropy.io import fits
+from astropy.wcs import WCS
 import matplotlib.pyplot as plt
 
 from scopesim.optics import image_plane_utils as imp_utils
@@ -330,3 +331,19 @@ def test_returns_expected_origin_fraction(self, x, y, frac):
     #     x, y = np.array([1.1, 2.9]), np.array([0.0, 0.5])
     #     xs, ys, fracs = imp_utils.sub_pixel_fractions(x, y)
     #     print(xs)
+
+
+class TestSkyDetWCS:
+    def test_wcs_roundtrip(self):
+        # FIXME: This should be a fixture, but everywhere in this module...
+        hdu = imo._image_hdu_rect()
+        sky_wcs = WCS(hdu.header)
+
+        # Scale 1.0 from _image_hdu_rect cdelt, 20 arbitrary plate scale
+        det_wcs, nax = imp_utils.det_wcs_from_sky_wcs(sky_wcs, 1.0, 20)
+        new_wcs, nax = imp_utils.sky_wcs_from_det_wcs(det_wcs, 1.0, 20,
+                                                      naxis=nax)
+
+        # Compare header representations because of missing NAXIS info in new
+        assert new_wcs.to_header() == sky_wcs.to_header()
+        assert all(nax == (hdu.header["NAXIS1"], hdu.header["NAXIS2"]))