Merge pull request #307 from AstarVienna/fh/header

Add basic sky coordinates (WCS) to ScopeSim output

scopesim/detector/detector.py

@@ -1,11 +1,12 @@
-from copy import deepcopy
+import logging
 import numpy as np
 
 from ..base_classes import ImagePlaneBase, DetectorBase
 from ..optics import image_plane_utils as imp_utils
 from .. import utils
 
 from astropy.io import fits
+from astropy.wcs import WCS
 
 
 class Detector(DetectorBase):
@@ -28,12 +29,22 @@ def extract_from(self, image_plane, spline_order=1, reset=True):
                                                        wcs_suffix="D")
 
     def reset(self):
-        self._hdu.data = np.zeros(self._hdu.data.shape)
+        self._hdu.data = np.zeros_like(self._hdu.data)
 
     @property
     def hdu(self):
         new_meta = utils.stringify_dict(self.meta)
         self._hdu.header.update(new_meta)
+
+        pixel_scale = utils.from_currsys("!INST.pixel_scale")
+        plate_scale = utils.from_currsys("!INST.plate_scale")
+        if pixel_scale == 0 or plate_scale == 0:
+            logging.warning("Could not create sky WCS.")
+        else:
+            sky_wcs, _ = imp_utils.sky_wcs_from_det_wcs(
+                WCS(self._hdu.header, key="D"), pixel_scale, plate_scale)
+            self._hdu.header.update(sky_wcs.to_header())
+
         return self._hdu
 
     @property

scopesim/detector/detector_array.py

@@ -1,27 +1,52 @@
+# -*- coding: utf-8 -*-
+"""Contains DetectorArray and aux functions."""
+
 import logging
 
 from astropy.io import fits
 
 from .detector import Detector
 
-from ..effects.effects_utils import get_all_effects
-from .. import effects as efs
 from .. import utils
 
 
 class DetectorArray:
+    """Manages the individual Detectors, mostly used for readout."""
+
+    # TODO: Should this class be called DetectorManager analogous to the
+    #       FOVManager and OpticsManager classes?
+    # TODO: Either way this could inherit from some collections.abc
+
     def __init__(self, detector_list=None, **kwargs):
         self.meta = {}
         self.meta.update(kwargs)
-        self.detector_list = detector_list
+
+        # The effect from which the instance is constructed
+        self._detector_list = detector_list
+
         self.array_effects = []
         self.dtcr_effects = []
         self.detectors = []
+
+        # The (initially empty) HDUList produced by the readout
         self.latest_exposure = None
 
-    def readout(self, image_planes, array_effects=None, dtcr_effects=None, **kwargs):
+    def readout(self, image_planes, array_effects=None, dtcr_effects=None,
+                **kwargs) -> fits.HDUList:
         """
-        Read out the detector array into a FITS file.
+        Read out the detector array into a FITS HDU List.
+
+        1. Select the relevant image plane to extract images from.
+        2. Apply detector array effects (apply to the entire image plane)
+        3. Make a series of Detectors for each row in a DetectorList object.
+        4. Iterate through all Detectors, extract image from image_plane.
+        5. Apply all effects (to all Detectors).
+        6. Add necessary header keywords (not implemented).
+        7. Generate a HDUList with the ImageHDUs and any extras:
+          - add ``PrimaryHDU`` with meta data regarding observation in header
+          - add ``ImageHDU`` objects
+          - add ``ASCIITableHDU`` with Effects meta data in final table
+            extension (not implemented)
 
         Parameters
         ----------
@@ -36,52 +61,43 @@ def readout(self, image_planes, array_effects=None, dtcr_effects=None, **kwargs)
 
         Returns
         -------
-        self.latest_exposure : fits.HDUList
+        latest_exposure : fits.HDUList
+            Output FITS HDU List.
 
         """
         # .. note:: Detector is what used to be called Chip
         #           DetectorArray is the old Detector
 
-        # 0. Select the relevant image plane to extract images from
-        # 1. make a series of Detectors for each row in a DetectorList object
-        # 2. iterate through all Detectors, extract image from image_plane
-        # 3. apply all effects (to all Detectors)
-        # 4. add necessary header keywords
-
-        # 5. Generate a HDUList with the ImageHDUs and any extras:
-        # - add PrimaryHDU with meta data regarding observation in header
-        # - add ImageHDUs
-        # - add ASCIITableHDU with Effects meta data in final table extension
-
         self.array_effects = array_effects or []
         self.dtcr_effects = dtcr_effects or []
         self.meta.update(kwargs)
 
-        # 0. Get the image plane that corresponds to this detector array
-        image_plane_id = self.detector_list.meta["image_plane_id"]
-        image_plane = [implane for implane in image_planes if
-                       implane.id == image_plane_id][0]
+        # 1. Get the image plane that corresponds to this detector array
+        # TODO: This silently only returns the first match, is that intended??
+        image_plane = next(implane for implane in image_planes if
+                           implane.id == self._detector_list.image_plane_id)
 
-        # 0a. Apply detector array effects (apply to the entire image plane)
+        # 2. Apply detector array effects (apply to the entire image plane)
         for effect in self.array_effects:
             image_plane = effect.apply_to(image_plane, **self.meta)
 
-        # 1. make a series of Detectors for each row in a DetectorList object
+        # 3. make a series of Detectors for each row in a DetectorList object
         self.detectors = [Detector(hdr, **self.meta)
-                          for hdr in self.detector_list.detector_headers()]
+                          for hdr in self._detector_list.detector_headers()]
 
-        # 2. iterate through all Detectors, extract image from image_plane
+        # 4. iterate through all Detectors, extract image from image_plane
+        logging.info("Extracting from %d detectors...", len(self.detectors))
         for detector in self.detectors:
             detector.extract_from(image_plane)
 
-            # 3. apply all effects (to all Detectors)
+            # 5. apply all effects (to all Detectors)
             for effect in self.dtcr_effects:
                 detector = effect.apply_to(detector)
 
-            # 4. add necessary header keywords
+            # 6. add necessary header keywords
             # .. todo: add keywords
 
-        # 5. Generate a HDUList with the ImageHDUs and any extras:
+        # 7. Generate a HDUList with the ImageHDUs and any extras:
         pri_hdu = make_primary_hdu(self.meta)
 
         # ..todo: effects_hdu unnecessary as long as make_effects_hdu does not do anything
@@ -100,11 +116,11 @@ def readout(self, image_planes, array_effects=None, dtcr_effects=None, **kwargs)
 
     def __repr__(self):
         msg = (f"{self.__class__.__name__}"
-               f"({self.detector_list!r}, **{self.meta!r})")
+               f"({self._detector_list!r}, **{self.meta!r})")
         return msg
 
     def __str__(self):
-        return f"{self.__class__.__name__} with {self.detector_list!s}"
+        return f"{self.__class__.__name__} with {self._detector_list!s}"
 
     def _repr_pretty_(self, p, cycle):
         """For ipython."""
@@ -116,23 +132,11 @@ def _repr_pretty_(self, p, cycle):
 
 def make_primary_hdu(meta):
     """Create the primary header from meta data."""
-    new_meta = utils.stringify_dict(meta)
     prihdu = fits.PrimaryHDU()
-    prihdu.header.update(new_meta)
-
+    prihdu.header.update(utils.stringify_dict(meta))
     return prihdu
 
 
 def make_effects_hdu(effects):
     # .. todo:: decide what goes into the effects table of meta data
     return fits.TableHDU()
-
-
-# def get_detector_list(effects):
-#     detector_lists = get_all_effects(effects, efs.DetectorList)
-#
-#     if len(detector_lists) != 1:
-#         logging.warning("None or more than one DetectorList found. Using the"
-#                       " first instance.{}".format(detector_lists))
-#
-#     return detector_lists[0]

scopesim/effects/detector_list.py

@@ -206,7 +206,7 @@ def image_plane_header(self):
 
         pixel_size = pixel_size.to(u.mm).value
         hdr = header_from_list_of_xy(x_det, y_det, pixel_size, "D")
-        hdr["IMGPLANE"] = self.meta["image_plane_id"]
+        hdr["IMGPLANE"] = self.image_plane_id
 
         return hdr
 
@@ -281,6 +281,11 @@ def plot(self, axes=None):
 
         return axes
 
+    @property
+    def image_plane_id(self) -> int:
+        """Get ID of the corresponding image plane."""
+        return self.meta["image_plane_id"]
+
 
 class DetectorWindow(DetectorList):
     """

scopesim/effects/psf_utils.py

@@ -97,9 +97,10 @@ def make_strehl_map_from_table(tbl, pixel_scale=1*u.arcsec):
                                          np.arange(-25, 26))).T,
                     method="nearest")
 
-    new_wcs, _ = imp_utils.create_wcs_from_points(np.array([[-25, -25],
-                                                            [25, 25]]),
-                                                  pixel_scale=1, arcsec=True)
+    new_wcs, _ = imp_utils.create_wcs_from_points(
+        np.array([[-25, -25],
+                  [25, 25]]) * u.arcsec,
+        pixel_scale=1*u.arcsec/u.pixel)
 
     map_hdu = fits.ImageHDU(header=new_wcs.to_header(), data=smap)