[DO NOT MERGE] Old MOS branch

scopesim/effects/spectral_trace_list_utils.py

@@ -15,7 +15,9 @@
 import numpy as np
 
 from scipy.interpolate import RectBivariateSpline
-from scipy.interpolate import interp1d
+from scipy.interpolate import InterpolatedUnivariateSpline, interp1d
+from scipy.ndimage import zoom
+from matplotlib import pyplot as plt
 
 from astropy.table import Table, vstack
 from astropy.modeling import fitting
@@ -24,8 +26,11 @@
 from astropy.wcs import WCS
 from astropy.modeling.models import Polynomial2D
 
-from ..utils import (power_vector, quantify, from_currsys, close_loop,
-                     figure_factory, get_logger)
+import pandas as pd
+
+from ..optics import image_plane_utils as imp_utils
+from ..utils import (deriv_polynomial2d, power_vector, interp2, check_keys,
+                     from_currsys, quantify, close_loop, figure_factory, get_logger)
 
 
 logger = get_logger(__name__)
@@ -132,7 +137,10 @@ def compute_interpolation_functions(self):
         self.xy2lam = Transform2D.fit(x_arr, y_arr, lam_arr)
         self.xilam2x = Transform2D.fit(xi_arr, lam_arr, x_arr)
         self.xilam2y = Transform2D.fit(xi_arr, lam_arr, y_arr)
+
+        self._xix2y = Transform2D.fit(xi_arr, x_arr, y_arr)
         self._xiy2x = Transform2D.fit(xi_arr, y_arr, x_arr)
+        self._xix2lam = Transform2D.fit(xi_arr, x_arr, lam_arr)
         self._xiy2lam = Transform2D.fit(xi_arr, y_arr, lam_arr)
 
         if self.dispersion_axis == "unknown":
@@ -423,6 +431,12 @@ def footprint(self, wave_min=None, wave_max=None, xi_min=None, xi_max=None):
             Minimum and maximum slit position on the sky.
             If `None`, use the full range that spectral trace is defined on.
             Float values are interpreted as arcsec.
+
+        Returns
+        -------
+        xy_edges : tuple of lists
+            (x_edges, y_edges) in [mm]
+
         """
         # Define the wavelength range of the footprint. This is a compromise
         # between the requested range (by method args) and the definition
@@ -673,8 +687,12 @@ def __init__(self, fov, dlam_per_pix):
         # dimensions are set by the slit aperture
         (n_lam, n_eta, n_xi) = fov.cube.data.shape
 
+        m_xi = int(np.round((fov.meta['xi_max'].to(u.arcsec).value -
+                             fov.meta['xi_min'].to(u.arcsec).value) / d_xi))
+        assert abs(m_xi - n_xi) <= 1
+
         # arrays of cube coordinates
-        cube_xi = d_xi * np.arange(n_xi) + fov.meta["xi_min"].value
+        cube_xi = d_xi * np.arange(n_xi) + fov.meta["xi_min"].to(u.arcsec).value
         cube_eta = d_eta * (np.arange(n_eta) - (n_eta - 1) / 2)
         cube_lam = wcs_lam.all_pix2world(np.arange(n_lam), 1)[0]
         cube_lam *= u.Unit(wcs_lam.wcs.cunit[0]).to(u.um)
@@ -983,8 +1001,6 @@ def make_image_interpolations(hdulist, **kwargs):
     return interps
 
 # ..todo: Check whether the following functions are actually used
-
-
 def rolling_median(x, n):
     """Calculate the rolling median of a sequence for +/- n entries."""
     y = [np.median(x[max(0, i-n):min(len(x), i+n+1)]) for i in range(len(x))]
@@ -1038,3 +1054,152 @@ def get_affine_parameters(coords):
     shears = (np.average(shears, axis=0) * rad2deg) - (90 + rotations)
 
     return rotations, shears
+
+
+class TracesHDUListGenerator:
+    def __init__(self, trace_hdus, aperture_ids, image_plane_ids):
+        """
+
+        Parameters
+        ----------
+        trace_hdus
+        aperture_ids
+        image_plane_ids
+        """
+        self.trace_hdus = trace_hdus
+        self.aperture_ids = aperture_ids
+        self.image_plane_ids = image_plane_ids
+
+    def make_hdulist(self, filename=None, overwrite=False):
+        trace_hdulist = fits.HDUList([self.pri_hdu, self.cat_hdu] +
+                                     self.trace_hdus)
+
+        if filename:
+            trace_hdulist.writeto(filename, overwrite=overwrite)
+
+        return trace_hdulist
+
+    @property
+    def hdulist(self):
+        return self.make_hdulist()
+
+    @property
+    def pri_hdu(self):
+        pri_hdu = fits.PrimaryHDU()
+        pri_hdu.header["ECAT"] = 1
+        pri_hdu.header["EDATA"] = 2
+
+        meta = {"author": "",
+                "source": "",
+                "descript": "Spectral Trace List",
+                "date-cre": "",
+                "date-mod": "",
+                }
+        pri_hdu.header.update(meta)
+
+        return pri_hdu
+
+    @property
+    def cat_hdu(self):
+        trace_names = [f"TRACE_{i}" for i in range(len(self.trace_hdus))]
+        cat_table = Table(names=["description", "extension_id",
+                                 "aperture_id", "image_plane_id"],
+                          data=[trace_names,
+                                2 + np.arange(len(self.trace_hdus)),
+                                self.aperture_ids,
+                                self.image_plane_ids
+                                ])
+        cat_hdu = fits.table_to_hdu(cat_table)
+        cat_hdu.header["EXTNAME"] = "TOC"
+
+        return cat_hdu
+
+
+def make_trace_hdu(const_wave_coords_dicts: list[dict],
+                   n_extra_points: int = 0,
+                   spline_order: int = 1):
+    """
+    Creates the BinTableHDU objects needed to make a spectral trace HDUList
+
+    Parameters
+    ----------
+    const_wave_coords_dicts : list of dicts
+        A list containing dictionaries for every line of constant wavelength
+        [{"wave": w [um], "x": list(xs) [mm],"y": list(ys) [mm]},
+         {..}, ..]
+    n_extra_points : int, 2-tuple
+        Default = 0. How many extra points to put in the spectral table
+        If a 2-tuple is passed, then n extra points will be added such:
+        (wavelength, spatial)
+    spline_order : int
+        Default = 1
+
+    Examples
+    --------
+    ::
+        dict_list = [{"wave": 1.9, "x": [0, 1], "y": [-2, -2]},
+                     {"wave": 2.4, "x": [0, 1], "y": [2, 2]}]
+        tg = TraceHDUGenerator(dict_list, n_extra_points=3)
+        tg.trace_hdu
+
+    """
+    wave, x, y = coords_from_lines_of_const_wavelength(const_wave_coords_dicts,
+                                                       n_extra_points,
+                                                       spline_order)
+    tbl = Table(names=["wavelength", "x", "y"],
+                units=[u.um, u.mm, u.mm],
+                data=[wave, x, y])
+
+    hdu = fits.table_to_hdu(tbl)
+    hdu.header[""]
+    return
+
+
+def coords_from_lines_of_const_wavelength(const_wave_coords_dicts,
+                                          n_extra_points=0, spline_order=1):
+    """
+    Returns expanded coordinates for constant wavelength in a spectral trace
+
+    Parameters
+    ----------
+    const_wave_coords_dicts : list of dicts
+        A list containing dictionaries for every line of constant wavelength
+        [{"wave": w [um], "x": list(xs) [mm],"y": list(ys) [mm]},
+         {..}, ..]
+    n_extra_points : int, 2-tuple
+        Default = 0. How many extra points to put in the spectral table
+        If a 2-tuple is passed, then n extra points will be added such:
+        (wavelength, spatial)
+    spline_order : int
+        Default = 1
+
+    Examples
+    --------
+    The following code creates a vertically bent trace with horizontal lines of
+    constant wavelength, with 20 additional points along the wavelength
+    direction, but 0 additional points along the slit direction
+    ::
+        dict_list = [{"wave": 1.9, "x": [0, 1, 2], "y": [-2.5, -2, -1.5]},
+                     {"wave": 2.15, "x": [1, 2, 3], "y": [0., 0, 0.]},
+                     {"wave": 2.4, "x": [0, 1, 2], "y": [1.5, 2, 2.5]}]
+        w, x, y = coords_from_lines_of_const_wavelength(dict_list,
+                                                        n_extra_points=(20, 0),
+                                                        spline_order=2)
+
+    """
+    # assert spline_order < len(const_wave_coords_dicts)
+    # assert spline_order < len(const_wave_coords_dicts[0]["x"])
+
+    x_im = np.array([dic["x"] for dic in const_wave_coords_dicts])
+    y_im = np.array([dic["y"] for dic in const_wave_coords_dicts])
+    w_im = np.array([[dic["wave"]]*len(dic["x"]) for dic in const_wave_coords_dicts])
+
+    cube = np.dstack([w_im, x_im, y_im])
+
+    xy_scale = 1 + np.array(n_extra_points) / np.array(x_im.shape)
+    if any(s != 1. for s in xy_scale):
+        cube = zoom(cube.astype(float), zoom=(*xy_scale, 1),
+                         order=spline_order)
+
+    ws, xs, ys = cube.reshape(-1, cube.shape[-1]).T
+    return ws, xs, ys

scopesim/optics/image_plane_utils.py

@@ -9,6 +9,7 @@
 from astropy.io import fits
 from astropy.table import Table
 from astropy.utils.exceptions import AstropyWarning
+from astropy.wcs import wcs
 from scipy import ndimage as ndi
 
 from ..utils import (unit_from_table, quantity_from_table, has_needed_keywords,
@@ -345,6 +346,22 @@ def _add_intpixel_sources_to_canvas(canvas_hdu, xpix, ypix, flux, mask):
 
 
 def _add_subpixel_sources_to_canvas(canvas_hdu, xpix, ypix, flux, mask):
+    """
+
+    Parameters
+    ----------
+    canvas_hdu : fits.ImageHDU
+    xpix, ypix : list of floats
+        [pix] Coordinates of sources in units of images pixel
+    flux : list of floats
+    mask : list of bools
+        flags for adding or ignoring a source
+
+    Returns
+    -------
+    canvas_hdu : fits.ImageHDU
+
+    """
     canvas_hdu.header["comment"] = f"Adding {len(flux)} sub-pixel files"
     canvas_shape = canvas_hdu.data.shape
     for xpx, ypx, flx, msk in zip(xpix, ypix, flux, mask):
@@ -1029,6 +1046,51 @@ def split_header(hdr, chunk_size, wcs_suffix=""):
     return hdr_list
 
 
+def extract_region_from_imagehdu(hdu, x_edges, y_edges, wcs_suffix=""):
+    """
+    Returns an ImageHDU for a subsection of the
+
+    Parameters
+    ----------
+    hdu
+    x_edges : list of floats
+        [deg, mm]
+    y_edges : list of floats
+        [deg, mm]
+
+    Returns
+    -------
+
+    """
+
+    s = "D" if wcs_suffix == "D" else " "
+    w = wcs.WCS(hdu.header, key=s, naxis=2)
+    xps, yps = w.all_world2pix([min(x_edges), max(x_edges)],
+                               [min(y_edges), max(y_edges)], 1)
+    #xps, yps = np.round(xps).astype(int), np.round(yps).astype(int)
+    xps = [int(np.floor(min(xps))), int(np.floor(max(xps)))]
+    yps = [int(np.floor(min(yps))), int(np.floor(max(yps)))]
+
+    w, h = hdu.header["NAXIS1"], hdu.header["NAXIS2"]
+    assert xps[0] >= 0 and xps[1] <= w and yps[0] >= 0 and yps[1] <= h, \
+           f"Pixel edges ({xps, yps}) were not within the bounds " \
+           f"of the array shape ({0, hdu.data.shape[-1], 0, hdu.data.shape[-2]})"
+
+    if hdu.header["NAXIS"] == 2:
+        new_data = hdu.data[yps[0]:yps[1]+1, xps[0]:xps[1]+1]
+    elif hdu.header["NAXIS"] == 3:
+        new_data = hdu.data[:, yps[0]:yps[1]+1, xps[0]:xps[1]+1]
+
+    new_hdu = fits.ImageHDU(data=new_data)
+    new_hdu.header.update(hdu.header)
+    new_hdu.header.update({"CRVAL1": np.average(x_edges),
+                           "CRVAL2": np.average(y_edges),
+                           "CRPIX1": np.average(xps),
+                           "CRPIX2": np.average(yps)})
+
+    return new_hdu
+
+
 def _fix_360(arr):
     """Fix the "full circle overflow" that occurs with deg."""
     if isinstance(arr, u.Quantity):

scopesim/rc.py

@@ -34,3 +34,9 @@
     Path(__pkg_dir__).absolute(),
     *[Path(pth).absolute() for pth in __config__["!SIM.file.search_path"]],
 ])
+
+# Needed by kl/mos_branch
+# Can refactor this later to remove it once things work again
+import os
+__basic_inst_path__ = os.path.abspath(os.path.join(__pkg_dir__,
+                                      "tests/mocks/basic_instrument/"))

scopesim/tests/conftest.py

@@ -91,6 +91,12 @@ def patch_mock_path_micado(mock_path_micado):
         yield
 
 
+@pytest.fixture(scope="package")
+def mock_path_basic_instrument():
+    """Path to MICADO mock files."""
+    return MOCK_DIR / "basic_instrument"
+
+
 @pytest.fixture(scope="function")
 def no_file_error():
     """Patch currsys to avoid missing file error."""

scopesim/tests/mocks/basic_instrument/INS_mos_apertures.dat

@@ -0,0 +1,11 @@
+# x_unit : arcsec
+# y_unit : arcsec
+# r_unit : arcsec
+# angle_unit : deg
+
+id  x       y       r       angle   conserve_image  shape
+0   -5      -5      0.4     0       False           round
+1   5       -5      0.6     0       False           rect
+2   0       0       0.8     0       False           hex
+3   -5      5       1       0       False           oct
+4   5       5       1.2     0       False           3

scopesim/tests/mocks/basic_instrument/INS_mos_bundle_ter_corrections.dat

@@ -0,0 +1,20 @@
+# source : MOSAIC MVP inst pkg
+# zenith_dist_unit : deg
+# wavelength_unit : um
+
+wavelength  zenith_dist     TRACE_ALL    TRACE_0    TRACE_1     TRACE_2     TRACE_3     TRACE_4     TRACE_5     TRACE_6
+0.920       0.              0.590        0.165      0.068       0.076       0.069       0.069       0.075       0.068
+1.202       0.              0.638        0.213      0.067       0.078       0.069       0.068       0.076       0.066
+1.638       0.              0.688        0.290      0.062       0.076       0.064       0.063       0.073       0.060
+0.920       15              0.577        0.157      0.068       0.075       0.069       0.068       0.074       0.067
+1.202       15              0.626        0.203      0.067       0.078       0.069       0.068       0.076       0.066
+1.638       15              0.678        0.277      0.063       0.076       0.065       0.064       0.073       0.061
+0.920       30              0.537        0.135      0.065       0.071       0.066       0.066       0.070       0.064
+1.202       30              0.588        0.174      0.066       0.075       0.067       0.067       0.073       0.065
+1.638       30              0.645        0.239      0.064       0.076       0.066       0.065       0.073       0.062
+0.920       45              0.461        0.102      0.058       0.063       0.059       0.059       0.062       0.058
+1.202       45              0.515        0.131      0.062       0.068       0.063       0.063       0.067       0.061
+1.638       45              0.578        0.180      0.063       0.073       0.065       0.064       0.071       0.062
+0.920       60              0.334        0.063      0.045       0.046       0.045       0.045       0.046       0.045
+1.202       60              0.381        0.077      0.050       0.052       0.050       0.050       0.052       0.049
+1.638       60              0.443        0.104      0.055       0.060       0.056       0.055       0.059       0.054

scopesim/tests/mocks/basic_instrument/YAML_basic_instrument.yaml

@@ -51,3 +51,11 @@ effects :
   include: "!OBS.include_slicer"
   kwargs:
     filename: "INS_ifu_apertures.dat"
+    extend_fov_beyond_slit: 5
+
+- name: fibre_list
+  description: collection of apertures corresponding to the MOS fibre heads
+  class: FibreApertureList
+  include: "!OBS.include_fibres"
+  kwargs:
+    filename: "INS_mos_apertures.dat"