From 727716a02e48b19c6728071a8f2123dfd11312d2 Mon Sep 17 00:00:00 2001 From: oczoske Date: Wed, 16 Jul 2025 18:39:09 +0200 Subject: [PATCH 01/11] start on MOSAIC spectral traces --- scopesim/effects/__init__.py | 1 + scopesim/effects/mosaic_trace_list.py | 56 +++++++++++++++++++++++++++ 2 files changed, 57 insertions(+) create mode 100644 scopesim/effects/mosaic_trace_list.py diff --git a/scopesim/effects/__init__.py b/scopesim/effects/__init__.py index b490c48d..14ab0fc7 100644 --- a/scopesim/effects/__init__.py +++ b/scopesim/effects/__init__.py @@ -8,6 +8,7 @@ from .spectral_trace_list import * from .spectral_efficiency import * from .metis_lms_trace_list import * +from .mosaic_trace_list import * from .surface_list import * from .ter_curves import * from . import ter_curves_utils diff --git a/scopesim/effects/mosaic_trace_list.py b/scopesim/effects/mosaic_trace_list.py new file mode 100644 index 00000000..e3ce2f28 --- /dev/null +++ b/scopesim/effects/mosaic_trace_list.py @@ -0,0 +1,56 @@ +# -*- coding: utf-8 -*- +"""SpectralTraceList and SpectralTrace for MOSAIC""" +import numpy as np +from astropy.table import Table +from astropy import units as u + +from .spectral_trace_list import SpectralTraceList +from .spectral_trace_list_utils import SpectralTrace + +from ..utils import get_logger, quantify + +logger = get_logger(__name__) + +class MosaicSpectralTraceList(SpectralTraceList): + """SpectralTraceList for MOSAIC""" + + def __init__(self, **kwargs): + super().__init__(**kwargs) + + self.aplist = self._file["Aperture List"].data + + self.view = np.array([(self.aplist["right"].max() - + self.aplist["left"].min()), + (self.aplist["top"].max() - + self.aplist["bottom"].min())]) + + def make_spectral_traces(self): + """Return a dictionary of spectral traces read in from a file.""" + self.ext_data = self._file[0].header["EDATA"] + self.ext_cat = self._file[0].header["ECAT"] + self.catalog = Table(self._file[self.ext_cat].data) + spec_traces = {} + for row in self.catalog: + if row["image_plane_id"] == -1: + continue + params = {col: row[col] for col in row.colnames} + params.update(self.meta) + hdu = self._file[row["extension_id"]] + spec_traces[row["description"]] = MosaicSpectralTrace(hdu, **params) + + self.spectral_traces = spec_traces + +class MosaicSpectralTrace(SpectralTrace): + + def __init__(self, trace_tbl, **kwargs): + super().__init__(trace_tbl, **kwargs) + + def compute_interpolation_functions(self): + print(self.table) + x_arr = self.table[self.meta["x_colname"]] + y_arr = self.table[self.meta["y_colname"]] + xi_arr = self.table[self.meta["s_colname"]] + lam_arr = self.table[self.meta["wave_colname"]] + + self.wave_min = quantify(np.min(lam_arr), u.um).value + self.wave_max = quantify(np.max(lam_arr), u.um).value From f520082cf37e5ab57dd5530b284c65d8b8d025ce Mon Sep 17 00:00:00 2001 From: oczoske Date: Mon, 21 Jul 2025 13:13:33 +0200 Subject: [PATCH 02/11] remove stray print --- scopesim/effects/mosaic_trace_list.py | 1 - 1 file changed, 1 deletion(-) diff --git a/scopesim/effects/mosaic_trace_list.py b/scopesim/effects/mosaic_trace_list.py index e3ce2f28..3e417cc0 100644 --- a/scopesim/effects/mosaic_trace_list.py +++ b/scopesim/effects/mosaic_trace_list.py @@ -46,7 +46,6 @@ def __init__(self, trace_tbl, **kwargs): super().__init__(trace_tbl, **kwargs) def compute_interpolation_functions(self): - print(self.table) x_arr = self.table[self.meta["x_colname"]] y_arr = self.table[self.meta["y_colname"]] xi_arr = self.table[self.meta["s_colname"]] From 309eb38fbc883bd21d026a24e646bd6c9dacd5da Mon Sep 17 00:00:00 2001 From: oczoske Date: Fri, 25 Jul 2025 17:58:33 +0200 Subject: [PATCH 03/11] Implement Transform1D, with tests --- scopesim/effects/mosaic_trace_list.py | 152 +++++++++++++++++- .../tests_effects/test_mosaic_trace_list.py | 64 ++++++++ 2 files changed, 214 insertions(+), 2 deletions(-) create mode 100644 scopesim/tests/tests_effects/test_mosaic_trace_list.py diff --git a/scopesim/effects/mosaic_trace_list.py b/scopesim/effects/mosaic_trace_list.py index 3e417cc0..762ad00f 100644 --- a/scopesim/effects/mosaic_trace_list.py +++ b/scopesim/effects/mosaic_trace_list.py @@ -1,13 +1,20 @@ # -*- coding: utf-8 -*- """SpectralTraceList and SpectralTrace for MOSAIC""" +from tqdm.auto import tqdm + import numpy as np from astropy.table import Table from astropy import units as u +from astropy.wcs import WCS +from astropy.modeling import fitting +from astropy.modeling.models import Polynomial1D from .spectral_trace_list import SpectralTraceList from .spectral_trace_list_utils import SpectralTrace -from ..utils import get_logger, quantify +from ..utils import get_logger, quantify, power_vector +from ..optics.fov import FieldOfView +from ..optics.fov_volume_list import FovVolumeList logger = get_logger(__name__) @@ -18,12 +25,87 @@ def __init__(self, **kwargs): super().__init__(**kwargs) self.aplist = self._file["Aperture List"].data - self.view = np.array([(self.aplist["right"].max() - self.aplist["left"].min()), (self.aplist["top"].max() - self.aplist["bottom"].min())]) + def apply_to(self, obj, **kwargs): + """See parent docstring.""" + ### This is copied from MetisSpectralTraceList, make less redundant? + if isinstance(obj, FovVolumeList): + logger.debug("Executing %s, FoV setup", self.meta['name']) + # Create a single volume that covers the aperture and + # the maximum wavelength range of the grating + volumes = [spectral_trace.fov_grid() + for spectral_trace in self.spectral_traces.values()] + wave_min = min(vol["wave_min"] for vol in volumes) + wave_max = max(vol["wave_max"] for vol in volumes) + extracted_vols = obj.extract(axes=["wave"], + edges=([[wave_min, wave_max]])) + obj.volumes = extracted_vols + print("Mosaic Spectral Trace List:", obj) + + if isinstance(obj, FieldOfView): + # Application to field of view + logger.debug("Executing %s, FoV", self.meta['name']) + if obj.hdu is not None and obj.hdu.header["NAXIS"] == 3: + obj.cube = obj.hdu + elif obj.hdu is None and obj.cube is None: + print("MosSpTrL: Making cube") + obj.cube = obj.make_cube_hdu() + + fovcube = obj.cube.data + n_z, n_y, n_x = fovcube.shape + fovwcs = WCS(obj.cube.header) + # Make this linear to avoid jump at RA 0 deg + fovwcs.wcs.ctype = ["LINEAR", "LINEAR", fovwcs.wcs.ctype[2]] + fovwcs_spat = fovwcs.sub(2) + + ## This is the place where we need to look at the apertures + ## - collapse each aperture to 1D spectrum by integrating spatially + ## - map each 1D spectrum to detector/fov + fovimage = np.zeros((obj.detector_header["NAXIS2"], + obj.detector_header["NAXIS1"]), + dtype=np.float32) + + for sptid, spt in tqdm(self.spectral_traces.items(), + desc="Fiber traces", position=2): + nx_slice = (self.aplist[sptid]["right"] - self.aplist[sptid]["left"])/pixscale + nx_slice = (self.aplist[sptid]["top"] - self.aplist[sptid]["bottom"])/pixscale + + ymin = spt.meta["fov"]["y_min"] + ymax = spt.meta["fov"]["y_max"] + + slicewcs = fovwcs.deepcopy() + slicewcs.wcs.ctype = ["LINEAR", "LINEAR", + slicewcs.wcs.ctype[2]] + slicewcs.wcs.crpix[0] = (nx_slice + 1) / 2 + slicewcs.wcs.crpix[1] = (ny_slice + 1) / 2 + slicewcs.wcs.crval[0] = (xmin + xmax) / 2 / 3600 + slicewcs.wcs.crval[1] = (ymin + ymax) / 2 / 3600 + slicewcs.wcs.cdelt[0] = (ymax - ymin) / ny_slice / 3600 + slicewcs.wcs.cdelt[1] = (ymax - ymin) / ny_slice / 3600 + slicewcs_spat = slicewcs.sub(2) + + # World coordinates for the slice + xworld, yworld = slicewcs_spat.app_pix2world(xslice, yslice, 0) + # FOV pixel coordinates for the slice + xfov, yfov = fovwcs_spat.all_world2pix(xworld, yworld, 0) + + for islice in range(n_z): + # this should not be necessary + ifov = RectBivariateSpline(np.arange(n_y), + np.arange(n_x), + fovcube[islice], kx=1, ky=1) + fovcube[islice, yfov, xfov].sum() + # build slicefov = FieldOfView1D(...) + # fovimage[] += slicefov.hdu.data (better integer row = slicefov...) + + return obj + + + def make_spectral_traces(self): """Return a dictionary of spectral traces read in from a file.""" self.ext_data = self._file[0].header["EDATA"] @@ -40,11 +122,14 @@ def make_spectral_traces(self): self.spectral_traces = spec_traces + class MosaicSpectralTrace(SpectralTrace): def __init__(self, trace_tbl, **kwargs): super().__init__(trace_tbl, **kwargs) + + def compute_interpolation_functions(self): x_arr = self.table[self.meta["x_colname"]] y_arr = self.table[self.meta["y_colname"]] @@ -53,3 +138,66 @@ def compute_interpolation_functions(self): self.wave_min = quantify(np.min(lam_arr), u.um).value self.wave_max = quantify(np.max(lam_arr), u.um).value + + + +class Transform1D(): + """ + 1-dimensional polynomial transform. + """ + + def __init__(self, coeffs, pretransform=None, + posttransform=None): + self.coeffs = np.asarray(coeffs) + self.nx = self.coeffs.shape[0] + self.pretransform = self._repackage(pretransform) + self.posttransform = self._repackage(posttransform) + + def _repackage(self, trafo): + """Make sure `trafo` is a tuple.""" + if trafo is not None and not isinstance(trafo, tuple): + trafo = (trafo, {}) + return trafo + + def __call__(self, x, **kwargs): + """ + Apply the polynomial transform. + + The transformation is a polynomial based on the simple monomials x^i. + """ + + if "pretransform" in kwargs: + self.pretransform = self._repackage(kwargs["pretransform"]) + if "postransform" in kwargs: + self.posttransform = self._repackage(kwargs["posttransform"]) + + x = np.array(x) + + # Apply pre transform + if self.pretransform is not None: + x = self.pretransform[0](x, **self.pretransform[1]) + + xvec = power_vector(x, self.nx - 1) + + result = self.coeffs @ xvec + + # Apply posttransform + if self.posttransform is not None: + result = self.posttransform[0](result, **self.posttransform[1]) + + return result + + @classmethod + def fit(cls, xin, xout, degree=4): + """Determine polynomial fit""" + pinit = Polynomial1D(degree=degree) + fitter = fitting.LinearLSQFitter() + fit = fitter(pinit, xin, xout) + return Transform1D(fit.parameters) + + def gradient(self): + """Compute the gradient of a 1d polynomial transformation""" + coeffs = self.coeffs + + dcoeffs = (coeffs * np.arange(self.nx))[1:] + return Transform1D(dcoeffs) diff --git a/scopesim/tests/tests_effects/test_mosaic_trace_list.py b/scopesim/tests/tests_effects/test_mosaic_trace_list.py new file mode 100644 index 00000000..96187e2d --- /dev/null +++ b/scopesim/tests/tests_effects/test_mosaic_trace_list.py @@ -0,0 +1,64 @@ +"""Unit tests for mosaic_trace_list.py""" + +# pylint: disable=missing-function-docstring +# pylint: disable=invalid-name +# pylint: disable=too-few-public-methods +import pytest + +import numpy as np + +from astropy.io import fits + +from scopesim.utils import power_vector +from scopesim.effects.mosaic_trace_list import Transform1D + +@pytest.fixture(name="tf1d", scope="class") +def fixture_tf1d(): + """Instantiate a Transform1D""" + coeffs = np.array([2, -1, 1]) + return Transform1D(coeffs) + +@pytest.fixture(name="quadratic", scope="class") +def fixture_quadratic(): + """Quadratic model, analytic and coeffients""" + coeffs = np.array([1, -1, 2]) + + def quadfunc(x): + z_a = 1 - 1 * x + 2 * x**2 + return z_a + + def dquad_dx(x): + return -1 + 4 * x + + return {'coeffs': coeffs, + 'function': quadfunc, + 'gradient': dquad_dx} + +class TestTransform1D: + """Tests for Transform1D()""" + def test_initialises_with_coeffs(self, tf1d): + assert isinstance(tf1d, Transform1D) + + def test_call_gives_correct_result(self, quadratic): + x = np.random.randn() + + # coefficients and explicit function + tf1d = Transform1D(quadratic['coeffs']) + assert tf1d(x) == quadratic['function'](x) + + def test_gradient_gives_correct_result(self, quadratic): + x = np.random.randn() + + tf2d = Transform1D(quadratic['coeffs']) + tf2d_grad = tf2d.gradient() + + assert tf2d_grad(x) == quadratic['gradient'](x) + + def test_fit_gives_correct_coeffs(self): + x = np.linspace(0, 1, 10) + y = 1 - 0.5 * x + 2.3 * x**2 - 3 * x**3 + + coeffs = np.array([1, -0.5, 2.3, -3]) + tf1d = Transform1D.fit(x, y, degree=3) + + assert tf1d.coeffs == pytest.approx(coeffs) From 32f623f6ff2779682046672703e27e7bb748e8e0 Mon Sep 17 00:00:00 2001 From: oczoske Date: Mon, 28 Jul 2025 17:47:55 +0200 Subject: [PATCH 04/11] first working version... --- scopesim/effects/mosaic_trace_list.py | 80 ++++++++++++++++----------- 1 file changed, 48 insertions(+), 32 deletions(-) diff --git a/scopesim/effects/mosaic_trace_list.py b/scopesim/effects/mosaic_trace_list.py index 762ad00f..a1b281b7 100644 --- a/scopesim/effects/mosaic_trace_list.py +++ b/scopesim/effects/mosaic_trace_list.py @@ -5,6 +5,7 @@ import numpy as np from astropy.table import Table from astropy import units as u +from astropy.io import fits from astropy.wcs import WCS from astropy.modeling import fitting from astropy.modeling.models import Polynomial1D @@ -61,6 +62,10 @@ def apply_to(self, obj, **kwargs): # Make this linear to avoid jump at RA 0 deg fovwcs.wcs.ctype = ["LINEAR", "LINEAR", fovwcs.wcs.ctype[2]] fovwcs_spat = fovwcs.sub(2) + fovwcs_spec = fovwcs.spectral + + det_header = obj.detector_header + naxis1d, naxis2d = det_header["NAXIS1"], det_header["NAXIS2"] ## This is the place where we need to look at the apertures ## - collapse each aperture to 1D spectrum by integrating spatially @@ -68,40 +73,46 @@ def apply_to(self, obj, **kwargs): fovimage = np.zeros((obj.detector_header["NAXIS2"], obj.detector_header["NAXIS1"]), dtype=np.float32) + pixscale = self.meta['pixel_scale'] + + image = np.zeros((naxis2d, naxis1d), dtype=np.float32) + imgwcs = WCS(naxis=2) + imgwcs.wcs.crpix = [fovwcs_spec.wcs.crpix[0], 5] + imgwcs.wcs.crval = [fovwcs_spec.wcs.crval[0], 1] + imgwcs.wcs.cdelt = [fovwcs_spec.wcs.cdelt[0], 1/5] + imgwcs.wcs.ctype = [fovwcs_spec.wcs.ctype[0], "LINEAR"] + imgwcs.wcs.cunit = [fovwcs_spec.wcs.cunit[0], ""] + ifib = 1 for sptid, spt in tqdm(self.spectral_traces.items(), desc="Fiber traces", position=2): - nx_slice = (self.aplist[sptid]["right"] - self.aplist[sptid]["left"])/pixscale - nx_slice = (self.aplist[sptid]["top"] - self.aplist[sptid]["bottom"])/pixscale - - ymin = spt.meta["fov"]["y_min"] - ymax = spt.meta["fov"]["y_max"] - - slicewcs = fovwcs.deepcopy() - slicewcs.wcs.ctype = ["LINEAR", "LINEAR", - slicewcs.wcs.ctype[2]] - slicewcs.wcs.crpix[0] = (nx_slice + 1) / 2 - slicewcs.wcs.crpix[1] = (ny_slice + 1) / 2 - slicewcs.wcs.crval[0] = (xmin + xmax) / 2 / 3600 - slicewcs.wcs.crval[1] = (ymin + ymax) / 2 / 3600 - slicewcs.wcs.cdelt[0] = (ymax - ymin) / ny_slice / 3600 - slicewcs.wcs.cdelt[1] = (ymax - ymin) / ny_slice / 3600 - slicewcs_spat = slicewcs.sub(2) - - # World coordinates for the slice - xworld, yworld = slicewcs_spat.app_pix2world(xslice, yslice, 0) - # FOV pixel coordinates for the slice - xfov, yfov = fovwcs_spat.all_world2pix(xworld, yworld, 0) - - for islice in range(n_z): - # this should not be necessary - ifov = RectBivariateSpline(np.arange(n_y), - np.arange(n_x), - fovcube[islice], kx=1, ky=1) - fovcube[islice, yfov, xfov].sum() - # build slicefov = FieldOfView1D(...) - # fovimage[] += slicefov.hdu.data (better integer row = slicefov...) - + theap = self.aplist[self.aplist['id'] == sptid] + nx_slice = (theap["right"] - theap["left"] ) / pixscale + ny_slice = (theap["top"] - theap["bottom"]) / pixscale + + print("SLICE:", theap["left"], theap["right"], theap["bottom"], theap["top"]) + # apertures are defined in arcsec. fovwcs is in degrees + xmin, xmax, ymin, ymax = (theap["left"]/3600, theap["right"]/3600, + theap["bottom"]/3600, theap["top"]/3600) + + imin = max(0, int(np.round(fovwcs_spat.all_world2pix(xmin, 0, 0)[0][0]))) + imax = int(np.round(fovwcs_spat.all_world2pix(xmax, 0, 0)[0][0])) + jmin = max(0, int(np.round(fovwcs_spat.all_world2pix(0, ymin, 0)[1][0]))) + jmax = int(np.round(fovwcs_spat.all_world2pix(0, ymax, 0)[1][0])) + print("SLICE:", imin, imax, jmin, jmax) + + # Sum over the spatial dimensions of the aperture + spec = fovcube[:, jmin:jmax, imin:imax].sum(axis=(1,2)) + nlam = len(spec) + # Need to interpolate this to the output wavelength grid + jfib = int(imgwcs.all_world2pix(fovwcs_spec.wcs.crval[0], ifib, 1)[1]) + print(sptid, ": Row ", jfib, ", Flux ", spec.sum()) + image[jfib,:nlam] += spec + ifib += 1 + + image_hdr = imgwcs.to_header() + image_hdr.extend(det_header) + obj.hdu = fits.ImageHDU(data=image, header=image_hdr) return obj @@ -133,13 +144,18 @@ def __init__(self, trace_tbl, **kwargs): def compute_interpolation_functions(self): x_arr = self.table[self.meta["x_colname"]] y_arr = self.table[self.meta["y_colname"]] - xi_arr = self.table[self.meta["s_colname"]] + #xi_arr = self.table[self.meta["s_colname"]] lam_arr = self.table[self.meta["wave_colname"]] self.wave_min = quantify(np.min(lam_arr), u.um).value self.wave_max = quantify(np.max(lam_arr), u.um).value + print("WAVEMIN:", self.wave_min) + print("WAVEMAX:", self.wave_max) + self.lam2x = Transform1D.fit(lam_arr, x_arr, degree=2) + self.x2lam = Transform1D.fit(x_arr, lam_arr, degree=2) + self.lam2y = Transform1D.fit(lam_arr, y_arr, degree=2) class Transform1D(): """ From 2e85f94c470bf851612cf9accb9bbe1fea332690 Mon Sep 17 00:00:00 2001 From: oczoske Date: Wed, 30 Jul 2025 10:52:21 +0200 Subject: [PATCH 05/11] just a debug message --- scopesim/effects/metis_lms_trace_list.py | 1 + 1 file changed, 1 insertion(+) diff --git a/scopesim/effects/metis_lms_trace_list.py b/scopesim/effects/metis_lms_trace_list.py index a07691b2..a9600fb5 100644 --- a/scopesim/effects/metis_lms_trace_list.py +++ b/scopesim/effects/metis_lms_trace_list.py @@ -84,6 +84,7 @@ def apply_to(self, obj, **kwargs): if isinstance(obj, FieldOfView): # Application to field of view + logger.debug("Executing %s, FoV", self.meta['name']) if obj.hdu is not None and obj.hdu.header["NAXIS"] == 3: obj.cube = obj.hdu elif obj.hdu is None and obj.cube is None: From 2136969bb8192fffa7e669cf58769229081ca95e Mon Sep 17 00:00:00 2001 From: oczoske Date: Wed, 30 Jul 2025 17:41:44 +0200 Subject: [PATCH 06/11] use interpolation to detector; average over apertures --- scopesim/effects/mosaic_trace_list.py | 27 ++++++++++++++------------- 1 file changed, 14 insertions(+), 13 deletions(-) diff --git a/scopesim/effects/mosaic_trace_list.py b/scopesim/effects/mosaic_trace_list.py index a1b281b7..7689fb79 100644 --- a/scopesim/effects/mosaic_trace_list.py +++ b/scopesim/effects/mosaic_trace_list.py @@ -9,7 +9,7 @@ from astropy.wcs import WCS from astropy.modeling import fitting from astropy.modeling.models import Polynomial1D - +from synphot import SourceSpectrum, Empirical1D from .spectral_trace_list import SpectralTraceList from .spectral_trace_list_utils import SpectralTrace @@ -45,7 +45,6 @@ def apply_to(self, obj, **kwargs): extracted_vols = obj.extract(axes=["wave"], edges=([[wave_min, wave_max]])) obj.volumes = extracted_vols - print("Mosaic Spectral Trace List:", obj) if isinstance(obj, FieldOfView): # Application to field of view @@ -53,7 +52,6 @@ def apply_to(self, obj, **kwargs): if obj.hdu is not None and obj.hdu.header["NAXIS"] == 3: obj.cube = obj.hdu elif obj.hdu is None and obj.cube is None: - print("MosSpTrL: Making cube") obj.cube = obj.make_cube_hdu() fovcube = obj.cube.data @@ -63,8 +61,11 @@ def apply_to(self, obj, **kwargs): fovwcs.wcs.ctype = ["LINEAR", "LINEAR", fovwcs.wcs.ctype[2]] fovwcs_spat = fovwcs.sub(2) fovwcs_spec = fovwcs.spectral + fovlam = fovwcs_spec.all_pix2world(np.arange(n_z), 0)[0] + fovlam <<= u.Unit(fovwcs_spec.wcs.cunit[0]) det_header = obj.detector_header + detwcs = WCS(det_header, key='D') naxis1d, naxis2d = det_header["NAXIS1"], det_header["NAXIS2"] ## This is the place where we need to look at the apertures @@ -90,7 +91,6 @@ def apply_to(self, obj, **kwargs): nx_slice = (theap["right"] - theap["left"] ) / pixscale ny_slice = (theap["top"] - theap["bottom"]) / pixscale - print("SLICE:", theap["left"], theap["right"], theap["bottom"], theap["top"]) # apertures are defined in arcsec. fovwcs is in degrees xmin, xmax, ymin, ymax = (theap["left"]/3600, theap["right"]/3600, theap["bottom"]/3600, theap["top"]/3600) @@ -99,15 +99,19 @@ def apply_to(self, obj, **kwargs): imax = int(np.round(fovwcs_spat.all_world2pix(xmax, 0, 0)[0][0])) jmin = max(0, int(np.round(fovwcs_spat.all_world2pix(0, ymin, 0)[1][0]))) jmax = int(np.round(fovwcs_spat.all_world2pix(0, ymax, 0)[1][0])) - print("SLICE:", imin, imax, jmin, jmax) - # Sum over the spatial dimensions of the aperture - spec = fovcube[:, jmin:jmax, imin:imax].sum(axis=(1,2)) - nlam = len(spec) + # Average over the spatial dimensions of the aperture (still per arcsec2) + fovflux = fovcube[:, jmin:jmax, imin:imax].mean(axis=(1,2)) + spec = SourceSpectrum(Empirical1D, points=fovlam.to(u.um), + lookup_table=fovflux) + # Need to interpolate this to the output wavelength grid + detlam = spt.x2lam(detwcs.all_pix2world(np.arange(naxis1d), 0, 0)[0]) + detlam <<= u.um jfib = int(imgwcs.all_world2pix(fovwcs_spec.wcs.crval[0], ifib, 1)[1]) - print(sptid, ": Row ", jfib, ", Flux ", spec.sum()) - image[jfib,:nlam] += spec + logger.debug("Flux from %s: %.4g", spt.trace_id, spec(detlam).value.sum()) + ### TODO: WHAT ABOUT THEM UNITS? + image[jfib,] += spec(detlam).value ifib += 1 image_hdr = imgwcs.to_header() @@ -150,9 +154,6 @@ def compute_interpolation_functions(self): self.wave_min = quantify(np.min(lam_arr), u.um).value self.wave_max = quantify(np.max(lam_arr), u.um).value - print("WAVEMIN:", self.wave_min) - print("WAVEMAX:", self.wave_max) - self.lam2x = Transform1D.fit(lam_arr, x_arr, degree=2) self.x2lam = Transform1D.fit(x_arr, lam_arr, degree=2) self.lam2y = Transform1D.fit(lam_arr, y_arr, degree=2) From 041a27848baa275221ca7f0d55bcf8919d485a31 Mon Sep 17 00:00:00 2001 From: oczoske Date: Thu, 31 Jul 2025 14:44:08 +0200 Subject: [PATCH 07/11] Units --- scopesim/effects/mosaic_trace_list.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/scopesim/effects/mosaic_trace_list.py b/scopesim/effects/mosaic_trace_list.py index 7689fb79..ce71095e 100644 --- a/scopesim/effects/mosaic_trace_list.py +++ b/scopesim/effects/mosaic_trace_list.py @@ -26,6 +26,7 @@ def __init__(self, **kwargs): super().__init__(**kwargs) self.aplist = self._file["Aperture List"].data + # TODO: check units or normalise to arcsec self.view = np.array([(self.aplist["right"].max() - self.aplist["left"].min()), (self.aplist["top"].max() - @@ -88,6 +89,9 @@ def apply_to(self, obj, **kwargs): for sptid, spt in tqdm(self.spectral_traces.items(), desc="Fiber traces", position=2): theap = self.aplist[self.aplist['id'] == sptid] + # solid angle in arcsec**2 + solid_angle = ((theap["right"] - theap["left"]) * + (theap["top"] - theap["bottom"])) nx_slice = (theap["right"] - theap["left"] ) / pixscale ny_slice = (theap["top"] - theap["bottom"]) / pixscale @@ -101,7 +105,7 @@ def apply_to(self, obj, **kwargs): jmax = int(np.round(fovwcs_spat.all_world2pix(0, ymax, 0)[1][0])) # Average over the spatial dimensions of the aperture (still per arcsec2) - fovflux = fovcube[:, jmin:jmax, imin:imax].mean(axis=(1,2)) + fovflux = fovcube[:, jmin:jmax, imin:imax].mean(axis=(1,2)) * solid_angle spec = SourceSpectrum(Empirical1D, points=fovlam.to(u.um), lookup_table=fovflux) @@ -111,10 +115,12 @@ def apply_to(self, obj, **kwargs): jfib = int(imgwcs.all_world2pix(fovwcs_spec.wcs.crval[0], ifib, 1)[1]) logger.debug("Flux from %s: %.4g", spt.trace_id, spec(detlam).value.sum()) ### TODO: WHAT ABOUT THEM UNITS? - image[jfib,] += spec(detlam).value + detdisp = np.diff(detlam, prepend=detlam[0]) + image[jfib,] += (spec(detlam) * detdisp).value ifib += 1 image_hdr = imgwcs.to_header() + image_hdr["BUNIT"] = "ph s-1" image_hdr.extend(det_header) obj.hdu = fits.ImageHDU(data=image, header=image_hdr) return obj From 489b22a3b9489ff344612dc36fcbe1aa2519a1aa Mon Sep 17 00:00:00 2001 From: oczoske Date: Mon, 11 Aug 2025 17:55:58 +0200 Subject: [PATCH 08/11] remove unused variables --- scopesim/effects/mosaic_trace_list.py | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/scopesim/effects/mosaic_trace_list.py b/scopesim/effects/mosaic_trace_list.py index ce71095e..b11c01c8 100644 --- a/scopesim/effects/mosaic_trace_list.py +++ b/scopesim/effects/mosaic_trace_list.py @@ -72,10 +72,6 @@ def apply_to(self, obj, **kwargs): ## This is the place where we need to look at the apertures ## - collapse each aperture to 1D spectrum by integrating spatially ## - map each 1D spectrum to detector/fov - fovimage = np.zeros((obj.detector_header["NAXIS2"], - obj.detector_header["NAXIS1"]), - dtype=np.float32) - pixscale = self.meta['pixel_scale'] image = np.zeros((naxis2d, naxis1d), dtype=np.float32) imgwcs = WCS(naxis=2) @@ -89,11 +85,10 @@ def apply_to(self, obj, **kwargs): for sptid, spt in tqdm(self.spectral_traces.items(), desc="Fiber traces", position=2): theap = self.aplist[self.aplist['id'] == sptid] + # solid angle in arcsec**2 solid_angle = ((theap["right"] - theap["left"]) * (theap["top"] - theap["bottom"])) - nx_slice = (theap["right"] - theap["left"] ) / pixscale - ny_slice = (theap["top"] - theap["bottom"]) / pixscale # apertures are defined in arcsec. fovwcs is in degrees xmin, xmax, ymin, ymax = (theap["left"]/3600, theap["right"]/3600, From 6a99a7b69219ae796bd71d2e1aaa9dc6d6fdf2d9 Mon Sep 17 00:00:00 2001 From: oczoske Date: Wed, 20 Aug 2025 14:57:56 +0200 Subject: [PATCH 09/11] Correctly use detwcs --- scopesim/effects/mosaic_trace_list.py | 13 ++++--------- 1 file changed, 4 insertions(+), 9 deletions(-) diff --git a/scopesim/effects/mosaic_trace_list.py b/scopesim/effects/mosaic_trace_list.py index b11c01c8..82c8e592 100644 --- a/scopesim/effects/mosaic_trace_list.py +++ b/scopesim/effects/mosaic_trace_list.py @@ -74,12 +74,6 @@ def apply_to(self, obj, **kwargs): ## - map each 1D spectrum to detector/fov image = np.zeros((naxis2d, naxis1d), dtype=np.float32) - imgwcs = WCS(naxis=2) - imgwcs.wcs.crpix = [fovwcs_spec.wcs.crpix[0], 5] - imgwcs.wcs.crval = [fovwcs_spec.wcs.crval[0], 1] - imgwcs.wcs.cdelt = [fovwcs_spec.wcs.cdelt[0], 1/5] - imgwcs.wcs.ctype = [fovwcs_spec.wcs.ctype[0], "LINEAR"] - imgwcs.wcs.cunit = [fovwcs_spec.wcs.cunit[0], ""] ifib = 1 for sptid, spt in tqdm(self.spectral_traces.items(), @@ -107,14 +101,15 @@ def apply_to(self, obj, **kwargs): # Need to interpolate this to the output wavelength grid detlam = spt.x2lam(detwcs.all_pix2world(np.arange(naxis1d), 0, 0)[0]) detlam <<= u.um - jfib = int(imgwcs.all_world2pix(fovwcs_spec.wcs.crval[0], ifib, 1)[1]) + yvals = spt.lam2y(detlam.value) + jfib = detwcs.all_world2pix(0, yvals.mean(), 0)[1].astype(int) logger.debug("Flux from %s: %.4g", spt.trace_id, spec(detlam).value.sum()) - ### TODO: WHAT ABOUT THEM UNITS? + detdisp = np.diff(detlam, prepend=detlam[0]) image[jfib,] += (spec(detlam) * detdisp).value ifib += 1 - image_hdr = imgwcs.to_header() + image_hdr = detwcs.to_header() image_hdr["BUNIT"] = "ph s-1" image_hdr.extend(det_header) obj.hdu = fits.ImageHDU(data=image, header=image_hdr) From 0de1720947093e406e2dcab79d3592f1ff53bb16 Mon Sep 17 00:00:00 2001 From: oczoske Date: Wed, 20 Aug 2025 17:06:11 +0200 Subject: [PATCH 10/11] Add effect to collapse fibre bundles to 1D spectrum --- scopesim/effects/mosaic_trace_list.py | 43 +++++++++++++++++++++++---- 1 file changed, 38 insertions(+), 5 deletions(-) diff --git a/scopesim/effects/mosaic_trace_list.py b/scopesim/effects/mosaic_trace_list.py index 82c8e592..a06a0b1d 100644 --- a/scopesim/effects/mosaic_trace_list.py +++ b/scopesim/effects/mosaic_trace_list.py @@ -1,6 +1,7 @@ # -*- coding: utf-8 -*- """SpectralTraceList and SpectralTrace for MOSAIC""" from tqdm.auto import tqdm +from typing import ClassVar import numpy as np from astropy.table import Table @@ -16,9 +17,13 @@ from ..utils import get_logger, quantify, power_vector from ..optics.fov import FieldOfView from ..optics.fov_volume_list import FovVolumeList +from ..detector import Detector logger = get_logger(__name__) + + + class MosaicSpectralTraceList(SpectralTraceList): """SpectralTraceList for MOSAIC""" @@ -56,7 +61,7 @@ def apply_to(self, obj, **kwargs): obj.cube = obj.make_cube_hdu() fovcube = obj.cube.data - n_z, n_y, n_x = fovcube.shape + n_z = fovcube.shape[0] fovwcs = WCS(obj.cube.header) # Make this linear to avoid jump at RA 0 deg fovwcs.wcs.ctype = ["LINEAR", "LINEAR", fovwcs.wcs.ctype[2]] @@ -75,7 +80,6 @@ def apply_to(self, obj, **kwargs): image = np.zeros((naxis2d, naxis1d), dtype=np.float32) - ifib = 1 for sptid, spt in tqdm(self.spectral_traces.items(), desc="Fiber traces", position=2): theap = self.aplist[self.aplist['id'] == sptid] @@ -107,7 +111,6 @@ def apply_to(self, obj, **kwargs): detdisp = np.diff(detlam, prepend=detlam[0]) image[jfib,] += (spec(detlam) * detdisp).value - ifib += 1 image_hdr = detwcs.to_header() image_hdr["BUNIT"] = "ph s-1" @@ -135,12 +138,11 @@ def make_spectral_traces(self): class MosaicSpectralTrace(SpectralTrace): + """A single spectral trace for MOSAIC""" def __init__(self, trace_tbl, **kwargs): super().__init__(trace_tbl, **kwargs) - - def compute_interpolation_functions(self): x_arr = self.table[self.meta["x_colname"]] y_arr = self.table[self.meta["y_colname"]] @@ -214,3 +216,34 @@ def gradient(self): dcoeffs = (coeffs * np.arange(self.nx))[1:] return Transform1D(dcoeffs) + + +class MosaicCollapseSpectralTraces(MosaicSpectralTraceList): + """Collapse SpectralTraces to 1D spectrum""" + required_keys = {"filename"} + z_order: ClassVar[tuple[int, ...]] = (899,) + + def __init__(self, **kwargs): + super().__init__(**kwargs) + + def apply_to(self, det, **kwargs): + """Apply to detector readout""" + if not isinstance(det, Detector): + print("Type of det:", type(det)) + return det + + image = det._hdu.data + detwcs = WCS(det._hdu.header, key='D') + spec = np.zeros(image.shape[1], dtype=np.float32) + for sptid, spt in tqdm(self.spectral_traces.items(), + desc="Fiber traces", position=2): + y_mm = spt.table['y'][0] + jfib = int(detwcs.all_world2pix(0, y_mm, 0)[1]) + spec += image[jfib,] + + x_mm = detwcs.all_pix2world(np.arange(image.shape[1]), 1, 0)[0] + lam = spt.x2lam(x_mm) + det._hdu = fits.BinTableHDU.from_columns([ + fits.Column(name='wavelength', format='D', array=lam, unit='um'), + fits.Column(name='spectrum', format='D', array=spec, unit='ADU')]) + return det From 9c315d9ad886c66cf72600fe5a6505bb4d290753 Mon Sep 17 00:00:00 2001 From: oczoske Date: Wed, 17 Sep 2025 14:18:09 +0200 Subject: [PATCH 11/11] PR comments --- scopesim/effects/mosaic_trace_list.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/scopesim/effects/mosaic_trace_list.py b/scopesim/effects/mosaic_trace_list.py index a06a0b1d..3fbf48ac 100644 --- a/scopesim/effects/mosaic_trace_list.py +++ b/scopesim/effects/mosaic_trace_list.py @@ -127,6 +127,10 @@ def make_spectral_traces(self): self.catalog = Table(self._file[self.ext_cat].data) spec_traces = {} for row in self.catalog: + # image_plane_id = -1 marks rows that should not be read, + # e.g. the aperture list. Although not necessary if the catalogue + # is formatted in a way that only traces are listed, this provides + # a possibility to "mask" traces. if row["image_plane_id"] == -1: continue params = {col: row[col] for col in row.colnames} @@ -229,7 +233,6 @@ def __init__(self, **kwargs): def apply_to(self, det, **kwargs): """Apply to detector readout""" if not isinstance(det, Detector): - print("Type of det:", type(det)) return det image = det._hdu.data