Add SourceField and subclasses, rework Source (#405)

scopesim/effects/fits_headers.py

@@ -9,6 +9,7 @@
 from astropy.table import Table
 
 from . import Effect
+from ..source.source_fields import HDUSourceField, TableSourceField
 from ..utils import from_currsys, find_file
 
 
@@ -530,19 +531,19 @@ def apply_to(self, hdul, **kwargs):
         if (src := opt_train._last_source) is not None:
             prefix = self.meta["keyword_prefix"]
             for i, field in enumerate(src.fields):
-                src_class = field.__class__.__name__
-                src_dic = deepcopy(src._meta_dicts[i])
-                if isinstance(field, fits.ImageHDU):
+                src_class = field.field.__class__.__name__
+                src_dic = deepcopy(src.meta)
+                if isinstance(field, HDUSourceField):
                     hdr = field.header
                     for key in hdr:
                         src_dic = {key: [hdr[key], hdr.comments[key]]}
 
-                elif isinstance(field, Table):
+                elif isinstance(field, TableSourceField):
                     src_dic.update(field.meta)
                     src_dic["length"] = len(field)
-                    for j, name in enumerate(field.colnames):
+                    for j, name in enumerate(field.field.colnames):
                         src_dic[f"col{j}_name"] = name
-                        src_dic[f"col{j}_unit"] = str(field[name].unit)
+                        src_dic[f"col{j}_unit"] = str(field.field[name].unit)
 
                 self.dict_list = [{"ext_number": self.meta["ext_number"],
                                    "keywords": {

scopesim/effects/spectral_efficiency.py

@@ -106,10 +106,10 @@ def apply_to(self, obj, **kwargs):
             logger.warning("No grating efficiency for trace %s", trace_id)
             return obj
 
-        wcs = WCS(obj.hdu.header).spectral
-        wave_cube = wcs.all_pix2world(np.arange(obj.hdu.data.shape[0]), 0)[0]
-        wave_cube = (wave_cube * u.Unit(wcs.wcs.cunit[0])).to(u.AA)
-        obj.hdu = apply_throughput_to_cube(obj.hdu, effic.throughput)
+        swcs = WCS(obj.hdu.header).spectral
+        with u.set_enabled_equivalencies(u.spectral()):
+            wave = swcs.pixel_to_world(np.arange(swcs.pixel_shape[0])) << u.um
+        obj.hdu = apply_throughput_to_cube(obj.hdu, effic.throughput, wave)
         return obj
 
     def plot(self):

scopesim/effects/ter_curves.py

@@ -16,6 +16,7 @@
 from ..base_classes import SourceBase, FOVSetupBase
 from ..optics.surface import SpectralSurface
 from ..source.source import Source
+from ..source.source_fields import CubeSourceField, SpectrumSourceField
 from ..utils import (from_currsys, quantify, check_keys, find_file,
                      figure_factory, get_logger)
 
@@ -113,13 +114,20 @@ def apply_to(self, obj, **kwargs):
             thru = self.throughput
 
             # apply transmission to source spectra
-            for isp, spec in enumerate(obj.spectra):
-                obj.spectra[isp] = combine_two_spectra(spec, thru, "multiply",
-                                                       wave_min, wave_max)
-
-            # apply transmission to cube fields
-            for icube, cube in enumerate(obj.cube_fields):
-                obj.cube_fields[icube] = apply_throughput_to_cube(cube, thru)
+            for fld in obj.fields:
+                if isinstance(fld, CubeSourceField):
+                    fld.field = apply_throughput_to_cube(
+                        fld.field, thru, fld.wave)
+                elif isinstance(fld, SpectrumSourceField):
+                    fld.spectra = {
+                        isp: combine_two_spectra(spec, thru, "multiply",
+                                                 wave_min, wave_max)
+                        for isp, spec in fld.spectra.items()
+                    }
+                else:
+                    # Rather log than raise here, can still move on
+                    logger.error("Source field is neither Cube nor has "
+                                 "spectra, this shouldn't occur...")
 
             # add the effect background to the source background field
             if self.background_source is not None:

scopesim/effects/ter_curves_utils.py

@@ -8,7 +8,7 @@
 from astropy.table import Table, QTable
 from astropy.io.votable import parse_single_table
 from astropy.io import ascii as ioascii
-from astropy.wcs import WCS
+from astropy.io import fits
 from synphot import SpectralElement, SourceSpectrum, Empirical1D, Observation
 from synphot.units import PHOTLAM
 
@@ -353,25 +353,32 @@ def scale_spectrum(spectrum, filter_name, amplitude):
     return spectrum
 
 
-def apply_throughput_to_cube(cube, thru):
+def apply_throughput_to_cube(
+    cube: fits.ImageHDU,
+    thru: SpectralElement | SourceSpectrum,
+    wave_cube: u.Quantity,
+) -> fits.ImageHDU:
     """
     Apply throughput curve to a spectroscopic cube.
 
     Parameters
     ----------
-    cube : ImageHDU
-         Three-dimensional image, dimension 0 (in python convention) is the
-         spectral dimension. WCS is required.
-    thru : synphot.SpectralElement, synphot.SourceSpectrum
+    cube : fits.ImageHDU
+        Three-dimensional image, dimension 0 (in python convention) is the
+        spectral dimension.
+    thru : SpectralElement | SourceSpectrum
+        Throughput curve, spectrum or filter.
+    wave_cube : u.Quantity["length"]
+        Wavelength axis of the cube.
 
     Returns
     -------
-    cube : ImageHDU, header unchanged, data multiplied with
-         wavelength-dependent throughput.
+    cube : fits.ImageHDU
+        Header unchanged, data multiplied with wavelength-dependent throughput.
+
     """
-    wcs = WCS(cube.header).spectral
-    wave_cube = wcs.all_pix2world(np.arange(cube.data.shape[0]), 0)[0]
-    wave_cube = (wave_cube * u.Unit(wcs.wcs.cunit[0])).to(u.AA)
+    # Note: wave_cube used to be converted to AA, but synphot understands um
+    #       or whatever just as well...
     cube.data *= thru(wave_cube).value[:, None, None]
     return cube
 

scopesim/optics/fov.py

@@ -121,7 +121,7 @@ def extract_from(self, src):
         if not isinstance(src, SourceBase):
             raise ValueError(f"source must be a Source object: {type(src)}")
 
-        fields_in_fov = [field for field in src.fields
+        fields_in_fov = [field.field for field in src.fields
                          if fu.is_field_in_fov(self.header, field)]
         if not fields_in_fov:
             logger.warning("No fields in FOV.")

scopesim/optics/fov_utils.py

@@ -7,6 +7,7 @@
 from synphot import SourceSpectrum, Empirical1D
 
 from . import image_plane_utils as imp_utils
+from ..source.source_fields import SourceField
 from ..utils import from_currsys, quantify, quantity_from_table, get_logger
 
 
@@ -21,8 +22,9 @@ def is_field_in_fov(fov_header, field, wcs_suffix=""):
     ----------
     fov_header : fits.Header
         Header from a FieldOfView object
-    field : [astropy.Table, astropy.ImageHDU]
-        Field object from a Source object
+    field : [SourceField, astropy.Table, astropy.ImageHDU]
+        Field object from a Source object. Should now be SourceField, but bare
+        Table and ImageHDU still supported.
     wcs_suffix : str
         ["S", "D"] Coordinate system: Sky or Detector
 
@@ -31,6 +33,11 @@ def is_field_in_fov(fov_header, field, wcs_suffix=""):
     is_inside_fov : bool
 
     """
+    if isinstance(field, SourceField):
+        field = field.field
+    # TODO: Check if this can always get a SourceField, if yes then just do
+    #       that and remove this.
+
     if isinstance(field, fits.ImageHDU) and \
             field.header.get("BG_SRC") is not None:
         is_inside_fov = True

scopesim/optics/image_plane_utils.py

@@ -52,7 +52,7 @@ def _get_headers(hdus_or_tables):
             else:
                 raise TypeError(
                     "hdu_or_table_list may only contain fits.ImageHDU, Table "
-                    "or fits.Header, found {type(hdu_or_table)}.")
+                    f"or fits.Header, found {type(hdu_or_table)}.")
         if tables:
             yield _make_bounding_header_for_tables(*tables,
                                                    pixel_scale=pixel_scale)

scopesim/optics/optical_train.py

@@ -7,6 +7,7 @@
 import numpy as np
 from scipy.interpolate import interp1d
 from astropy import units as u
+from astropy.wcs import WCS
 
 from tqdm import tqdm
 
@@ -20,7 +21,7 @@
 from ..detector import DetectorManager
 from ..effects import ExtraFitsKeywords
 from ..utils import from_currsys, top_level_catch, get_logger
-from .. import rc, __version__
+from .. import __version__
 
 
 logger = get_logger(__name__)
@@ -275,9 +276,10 @@ def prepare_source(self, source):
         ImageHDU.
         """
         # Convert to PHOTLAM per arcsec2
-        # ..todo: this is not sufficiently general
+        # TODO: this is not sufficiently general
+        # TODO: Maybe move this to source_fields??
 
-        for ispec, spec in enumerate(source.spectra):
+        for ispec, spec in source.spectra.items():
             # Put on fov wavegrid
             wave_min = min(fov.meta["wave_min"] for fov in self.fov_manager.fovs)
             wave_max = max(fov.meta["wave_max"] for fov in self.fov_manager.fovs)
@@ -348,6 +350,8 @@ def prepare_source(self, source):
             cube.header["CDELT3"] = dwave
             cube.header["CUNIT3"] = wave_unit.name
 
+            cube.wcs = WCS(cube.field)
+
         return source
 
     @top_level_catch

scopesim/rc.py

@@ -1,5 +1,5 @@
 # -*- coding: utf-8 -*-
-"""Global configurations for ScopeSim."""
+"""Global configurations for ScopeSim (rc ... runtime configuration)."""
 
 from pathlib import Path
 import yaml