diff --git a/scopesim/effects/__init__.py b/scopesim/effects/__init__.py index da03ffc0..170b4192 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 .micado_ifu_trace_list import * from .mosaic_trace_list import * from .surface_list import * from .ter_curves import * diff --git a/scopesim/effects/apertures.py b/scopesim/effects/apertures.py index fd6984a5..b8ca2ce2 100644 --- a/scopesim/effects/apertures.py +++ b/scopesim/effects/apertures.py @@ -118,6 +118,7 @@ def apply_to(self, obj, **kwargs): """See parent docstring.""" if isinstance(obj, FovVolumeList): logger.debug("Executing %s, FoV setup", self.meta['name']) + x = quantity_from_table("x", self.table, u.arcsec).to_value(u.arcsec) y = quantity_from_table("y", self.table, @@ -129,6 +130,12 @@ def apply_to(self, obj, **kwargs): vol["meta"]["xi_min"] = min(x) * u.arcsec vol["meta"]["xi_max"] = max(x) * u.arcsec + if "offset_type" in self.meta: + vol["meta"]["offset_type"] = self.meta["offset_type"] + vol["meta"]["offset_x"] = self.meta["offset_x"] + vol["meta"]["offset_y"] = self.meta["offset_y"] + else: + vol["meta"]["offset_type"] = "none" return obj # Outdated. Remove when removing all old FOVManager code from effects diff --git a/scopesim/effects/metis_lms_trace_list.py b/scopesim/effects/metis_lms_trace_list.py index bf8c0d88..e6770a90 100644 --- a/scopesim/effects/metis_lms_trace_list.py +++ b/scopesim/effects/metis_lms_trace_list.py @@ -53,93 +53,98 @@ def __init__(self, **kwargs): self.slicelist["top"].max() - self.slicelist["bottom"].min()]) - def apply_to(self, obj, **kwargs): - """See parent docstring.""" - if isinstance(obj, FovVolumeList): - # Create a single volume that covers the aperture and - # the maximum wavelength range of LMS - 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 - - 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: - obj.cube = obj.make_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) - ny_slice = self.meta["slice_samples"] - - # Spatial pixel coordinates - xslice, yslice = np.meshgrid(np.arange(n_x), - np.arange(ny_slice)) - - fovimage = np.zeros((obj.detector_header["NAXIS2"], - obj.detector_header["NAXIS1"]), - dtype=np.float32) - - for sptid, spt in tqdm(self.spectral_traces.items(), - desc=" Spectral Traces", position=2): - 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[1] = (ny_slice + 1) / 2 - slicewcs.wcs.crval[1] = (ymin + ymax) / 2 / 3600 - slicewcs.wcs.cdelt[1] = (ymax - ymin) / ny_slice / 3600 - slicewcs_spat = slicewcs.sub(2) - - # World coordinates for the slice - xworld, yworld = slicewcs_spat.all_pix2world(xslice, - yslice, 0) - # FOV pixel coordinates for the slice - xfov, yfov = fovwcs_spat.all_world2pix(xworld, yworld, 0) - - slicecube = np.zeros((n_z, ny_slice, n_x)) - for islice in range(n_z): - ifov = RectBivariateSpline(np.arange(n_y), - np.arange(n_x), - fovcube[islice], kx=1, ky=1) - slicecube[islice] = ifov(yfov, xfov, grid=False) - - slicefov = FieldOfView3D(obj.header, - [obj.meta["wave_min"], - obj.meta["wave_max"]]) - slicefov.detector_header = obj.detector_header - slicefov.meta["xi_min"] = obj.meta["xi_min"] - slicefov.meta["xi_max"] = obj.meta["xi_max"] - slicefov.meta["trace_id"] = sptid - slicefov.cube = fits.ImageHDU(header=slicewcs.to_header(), - data=slicecube) - # slicefov.cube.writeto(f"slicefov_{sptid}.fits", overwrite=True) - slicefov.hdu = spt.map_spectra_to_focal_plane(slicefov) - if slicefov.hdu is not None: - sxmin = slicefov.hdu.header["XMIN"] - sxmax = slicefov.hdu.header["XMAX"] - symin = slicefov.hdu.header["YMIN"] - symax = slicefov.hdu.header["YMAX"] - fovimage[symin:symax, sxmin:sxmax] += slicefov.hdu.data - - obj.hdu = fits.ImageHDU(data=fovimage, header=obj.detector_header) + # Overrides for the two apply_to() methods + def apply_to_fovvolumelist(self, obj): + """Setup + Create a single volume that covers the aperture and + the maximum wavelength range of LMS. + """ + logger.debug("Executing %s: Setup", self.meta['name']) + 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 + return obj + + def apply_to_fov(self, obj): + """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: + obj.cube = obj.make_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) + ny_slice = self.meta["slice_samples"] + + # Spatial pixel coordinates + xslice, yslice = np.meshgrid(np.arange(n_x), + np.arange(ny_slice)) + + fovimage = np.zeros((obj.detector_header["NAXIS2"], + obj.detector_header["NAXIS1"]), + dtype=np.float32) + + for sptid, spt in tqdm(self.spectral_traces.items(), + desc=" Spectral Traces", position=2): + 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[1] = (ny_slice + 1) / 2 + slicewcs.wcs.crval[1] = (ymin + ymax) / 2 / 3600 + slicewcs.wcs.cdelt[1] = (ymax - ymin) / ny_slice / 3600 + slicewcs_spat = slicewcs.sub(2) + + # World coordinates for the slice + xworld, yworld = slicewcs_spat.all_pix2world(xslice, + yslice, 0) + # FOV pixel coordinates for the slice + xfov, yfov = fovwcs_spat.all_world2pix(xworld, yworld, 0) + + slicecube = np.zeros((n_z, ny_slice, n_x)) + for islice in range(n_z): + ifov = RectBivariateSpline(np.arange(n_y), + np.arange(n_x), + fovcube[islice], kx=1, ky=1) + slicecube[islice] = ifov(yfov, xfov, grid=False) + + slicefov = FieldOfView3D(obj.header, + [obj.meta["wave_min"], + obj.meta["wave_max"]]) + slicefov.detector_header = obj.detector_header + slicefov.meta["xi_min"] = obj.meta["xi_min"] + slicefov.meta["xi_max"] = obj.meta["xi_max"] + slicefov.meta["trace_id"] = sptid + slicefov.cube = fits.ImageHDU(header=slicewcs.to_header(), + data=slicecube) + # slicefov.cube.writeto(f"slicefov_{sptid}.fits", overwrite=True) + slicefov.hdu = spt.map_spectra_to_focal_plane(slicefov) + if slicefov.hdu is not None: + sxmin = slicefov.hdu.header["XMIN"] + sxmax = slicefov.hdu.header["XMAX"] + symin = slicefov.hdu.header["YMIN"] + symax = slicefov.hdu.header["YMAX"] + fovimage[symin:symax, sxmin:sxmax] += slicefov.hdu.data + + obj.hdu = fits.ImageHDU(data=fovimage, header=obj.detector_header) return obj + def make_spectral_traces(self): """Compute the transformations by interpolation.""" # nslice = len(self._file["Aperture List"].data) diff --git a/scopesim/effects/micado_ifu_trace_list.py b/scopesim/effects/micado_ifu_trace_list.py new file mode 100644 index 00000000..cd2595b4 --- /dev/null +++ b/scopesim/effects/micado_ifu_trace_list.py @@ -0,0 +1,217 @@ +# -*- coding: utf-8 -*- +"""SpectralTraceList and SpectralTrace for the MICADO IFU.""" + +import copy +import warnings + +from tqdm.auto import tqdm +import numpy as np +from scipy.interpolate import RectBivariateSpline + +from astropy.io import fits +from astropy.io import ascii as ioascii +from astropy.table import Table +from astropy.wcs import WCS +from astropy import units as u + +from ..utils import from_currsys, find_file, quantify, get_logger +from .spectral_trace_list import SpectralTraceList +from .spectral_trace_list_utils import SpectralTrace +from .spectral_trace_list_utils import Transform2D +from .spectral_trace_list_utils import make_image_interpolations, det_offset +from .apertures import ApertureMask +from .ter_curves import TERCurve +from ..optics.fov import FieldOfView, FieldOfView3D +from ..optics.fov_volume_list import FovVolumeList + +logger = get_logger(__name__) + + +class MicadoIFUSpectralTraceList(SpectralTraceList): + """SpectralTraceList for the MICADO IFU.""" + + _class_params = { + "naxis1": 112, + "nslice": 32, + "slicewidth": 0.012, # arcsec + "wave_min": 1.5, + "wave_max": 2.5, # todo: do not hardcode + } + + + def __init__(self, **kwargs): + super().__init__(**kwargs) + self.wave_min = self.meta["wave_min"] + self.wave_max = self.meta["wave_max"] + + def apply_to_fovvolumelist(self, obj): + """ + Setup + + Create a single volume that covers the aperture and the maximum + wavelength range of the MICADO IFU. + """ + logger.debug("Executing %s: Setup", self.meta['name']) + 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 + self.wave_min = wave_min + self.wave_max = wave_max + return obj + + def apply_to_fov(self, obj): + """Apply to FieldOfView""" + + 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: + obj.cube = obj.make_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) + ny_slice = self.meta["slice_samples"] + + # Spatial pixel coordinates + xslice, yslice = np.meshgrid(np.arange(n_x), + np.arange(ny_slice)) + + fovimage = np.zeros((obj.detector_header["NAXIS2"], + obj.detector_header["NAXIS1"]), + dtype=np.float32) + + for sptid, spt in tqdm(self.spectral_traces.items(), + desc=" Spectral Traces", position=2): + 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[1] = (ny_slice + 1) / 2 + slicewcs.wcs.crval[1] = (ymin + ymax) / 2 / 3600 + slicewcs.wcs.cdelt[1] = (ymax - ymin) / ny_slice / 3600 + slicewcs_spat = slicewcs.sub(2) + + # World coordinates for the slice + xworld, yworld = slicewcs_spat.all_pix2world(xslice, + yslice, 0) + # FOV pixel coordinates for the slice + xfov, yfov = fovwcs_spat.all_world2pix(xworld, yworld, 0) + + slicecube = np.zeros((n_z, ny_slice, n_x)) + for islice in range(n_z): + ifov = RectBivariateSpline(np.arange(n_y), + np.arange(n_x), + fovcube[islice], kx=1, ky=1) + slicecube[islice] = ifov(yfov, xfov, grid=False) + + slicefov = FieldOfView3D(obj.header, + [obj.meta["wave_min"], + obj.meta["wave_max"]]) + slicefov.detector_header = obj.detector_header + slicefov.meta["xi_min"] = spt.meta["fov"]["xi_min"] + slicefov.meta["xi_max"] = spt.meta["fov"]["xi_max"] + slicefov.meta["trace_id"] = sptid + slicefov.cube = fits.ImageHDU(header=slicewcs.to_header(), + data=slicecube) + # slicefov.cube.writeto(f"slicefov_{sptid}.fits", overwrite=True) + slicefov.hdu = spt.map_spectra_to_focal_plane(slicefov) + if slicefov.hdu is not None: + sxmin = slicefov.hdu.header["XMIN"] + sxmax = slicefov.hdu.header["XMAX"] + symin = slicefov.hdu.header["YMIN"] + symax = slicefov.hdu.header["YMAX"] + fovimage[symin:symax, sxmin:sxmax] += slicefov.hdu.data + + obj.hdu = fits.ImageHDU(data=fovimage, header=obj.detector_header) + return obj + + + + + def make_spectral_traces(self): + """Make a spectral trace for each combination of order and aperture""" + 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) + self.slicelist = self._file["Aperture List"].data + spec_traces = {} + for row in self.catalog: + if row["image_plane_id"] == -99: + continue + params = {col: row[col] for col in row.colnames} + params.update(self.meta) + trace_hdu = self._file[row["extension_id"]] + + for spslice, apid in enumerate(self.slicelist["id"]): + specid = f"{row['description']}_{apid:02d}" + thistrace = MicadoIFUSpectralTrace(trace_hdu, self.slicelist, spslice, + **params) + thistrace.meta["trace_id"] = specid + spec_traces[specid] = thistrace + + self.spectral_traces = spec_traces + + + +class MicadoIFUSpectralTrace(SpectralTrace): + """SpectralTrace for the MICADO IFU""" + + _class_params = { + "naxis1": 112, + "nslice": 32, + "slicewidth": 0.012, # arcsec + "pixsize": 0.015, # mm + } + + + def __init__(self, trace_tbl, aplist, spslice, cmds=None, **kwargs): + super().__init__(trace_tbl, **kwargs) + + self._set_dispersion(self.wave_min, self.wave_max) + # Provisional + self.aplist = aplist + self.meta["slice"] = spslice + self.meta["fov"] = self.fov_grid() + + lam = np.linspace(self.wave_min, self.wave_max, 1001) * u.um + off_y = aplist['offset'][spslice] + dlam = np.median(self.dlam_per_pix(lam)) / self.meta["pixsize"] * off_y + self.xy2lam.posttransform = (det_offset, {"offset": -dlam}) + self.xilam2x.pretransform_y = (det_offset, {"offset": -dlam}) + self.xilam2y.pretransform_y = (det_offset, {"offset": -dlam}) + + def fov_grid(self): + """ + Provide information on the source space volume required by the effect. + """ + aperture = self.aplist[self.meta["slice"]] + x_min = aperture["left"] + x_max = aperture["right"] + y_min = aperture["bottom"] + y_max = aperture["top"] + xi_min = aperture["xi1"] + xi_max = aperture["xi2"] + + # ..todo: just a hack - xi and x are the same except xi is a quantity + xi_min = quantify(xi_min, u.arcsec) + xi_max = quantify(xi_max, u.arcsec) + + return {"x_min": x_min, "x_max": x_max, + "y_min": y_min, "y_max": y_max, + "xi_min": xi_min, "xi_max": xi_max, + "wave_min": self.meta["wave_min"], + "wave_max": self.meta["wave_max"], + "trace_id": self.trace_id} diff --git a/scopesim/effects/spectral_trace_list.py b/scopesim/effects/spectral_trace_list.py index d455fc08..4422e183 100644 --- a/scopesim/effects/spectral_trace_list.py +++ b/scopesim/effects/spectral_trace_list.py @@ -16,7 +16,8 @@ from .effects import Effect from .ter_curves import FilterCurve -from .spectral_trace_list_utils import SpectralTrace, make_image_interpolations +from .spectral_trace_list_utils import (SpectralTrace, make_image_interpolations, + det_offset) from ..optics.image_plane_utils import header_from_list_of_xy from ..optics.fov import FieldOfView from ..optics.fov_volume_list import FovVolumeList @@ -125,8 +126,6 @@ def __init__(self, **kwargs): "x_colname": "x", "y_colname": "y", "s_colname": "s", - "offset_x": 0, # [mm] in detector plane - "offset_y": 0, # [mm] in detector plane "wave_colname": "wavelength", "center_on_wave_mid": False, "dwave": 0.002, # [um] for finding the best fit dispersion @@ -149,6 +148,8 @@ def make_spectral_traces(self): self.catalog = Table(self._file[self.ext_cat].data) spec_traces = {} for row in self.catalog: + if row["image_plane_id"] == -99: + continue params = {col: row[col] for col in row.colnames} params.update(self.meta) hdu = self._file[row["extension_id"]] @@ -183,6 +184,7 @@ def update_meta(self): self.meta["y_min"] = min(ylim) self.meta["y_max"] = max(ylim) + def apply_to(self, obj, **kwargs): """ Interface between ``FieldOfView`` and ``SpectralTraceList``. @@ -197,71 +199,118 @@ def apply_to(self, obj, **kwargs): list, identified by meta["trace_id"]. """ if isinstance(obj, FovVolumeList): - logger.debug("%s applied to %s", self.display_name, - obj.__class__.__name__) - # Setup of FieldOfView object - # volumes = [spectral_trace.fov_grid() - # for spectral_trace in self.spectral_traces.values()] - - new_vols_list = [] - - # for vol in volumes: - for spt in self.spectral_traces.values(): - vol = spt.fov_grid() - wave_edges = [vol["wave_min"], vol["wave_max"]] - if "x_min" in vol: - x_edges = [vol["x_min"], vol["x_max"]] - y_edges = [vol["y_min"], vol["y_max"]] - extracted_vols = obj.extract( - axes=["wave", "x", "y"], - edges=(wave_edges, x_edges, y_edges), - aperture_id=vol["aperture_id"]) - else: - extracted_vols = obj.extract( - axes=["wave"], - edges=(wave_edges, ), - aperture_id=vol["aperture_id"]) - - for ex_vol in extracted_vols: - ex_vol["meta"].update(vol) - ex_vol["meta"].pop("wave_min") - ex_vol["meta"].pop("wave_max") - new_vols_list.extend(extracted_vols) - - obj.volumes = new_vols_list + obj = self.apply_to_fovvolumelist(obj) if isinstance(obj, FieldOfView): - logger.debug("%s applied to %s", self.display_name, - obj.__class__.__name__) - # Application to field of view - if obj.hdu is not None and obj.hdu.header["NAXIS"] == 3: - obj.cube = obj.hdu - elif obj.hdu is not None and obj.hdu.header["NAXIS"] == 2: - # todo: catch the case of obj.hdu.header["NAXIS"] == 2 - # for MAAT - pass - elif obj.hdu is None and obj.cube is None: - logger.info("Making cube") - obj.cube = obj.make_hdu() - - # Check whether an offset slit is used. If so, recompute spectral traces. - offset_x = obj.cube.header["CRVAL1D"] - offset_y = obj.cube.header["CRVAL2D"] - if (offset_x != self.meta["offset_x"] or - offset_y != self.meta["offset_y"]): - logger.debug("Recomputing spectral traces for offset (%.1g, %.1g)", - offset_x, offset_y) - self.meta["offset_x"] = offset_x - self.meta["offset_y"] = offset_y - self.make_spectral_traces() - self.update_meta() - - spt = self.spectral_traces[obj.trace_id] - obj.hdu = spt.map_spectra_to_focal_plane(obj) + obj = self.apply_to_fov(obj) logger.debug("%s done", self.display_name) return obj + def apply_to_fovvolumelist(self, obj): + """ + Setup of SpectralTraceList + + During setup of the required FieldOfView objects, the + SpectralTraceList is asked for the source space volumes that + it requires (spatial limits and wavelength limits). + """ + logger.debug("%s applied to %s", self.display_name, + obj.__class__.__name__) + + new_vols_list = [] + + for spt in self.spectral_traces.values(): + vol = spt.fov_grid() + wave_edges = [vol["wave_min"], vol["wave_max"]] + if "x_min" in vol: + x_edges = [vol["x_min"], vol["x_max"]] + y_edges = [vol["y_min"], vol["y_max"]] + extracted_vols = obj.extract( + axes=["wave", "x", "y"], + edges=(wave_edges, x_edges, y_edges), + aperture_id=vol["aperture_id"]) + else: + extracted_vols = obj.extract( + axes=["wave"], + edges=(wave_edges, ), + aperture_id=vol["aperture_id"]) + + for ex_vol in extracted_vols: + ex_vol["meta"].update(vol) + ex_vol["meta"].pop("wave_min") + ex_vol["meta"].pop("wave_max") + new_vols_list.extend(extracted_vols) + + obj.volumes = new_vols_list + return obj + + + def apply_to_fov(self, obj): + """Apply to FieldOfView + + During "observation" the method is passed a single FieldOfView + object and applies the mapping to the image plane to it. + The FieldOfView object is associated to one SpectralTrace from the + list, identified by meta["trace_id"]. + """ + logger.debug("%s applied to %s %s", self.display_name, + obj.__class__.__name__, obj.trace_id) + + if obj.hdu is not None and obj.hdu.header["NAXIS"] == 3: + obj.cube = obj.hdu + elif obj.hdu is not None and obj.hdu.header["NAXIS"] == 2: + # todo: catch the case of obj.hdu.header["NAXIS"] == 2 + # for MAAT + pass + elif obj.hdu is None and obj.cube is None: + logger.info("Making cube") + obj.cube = obj.make_hdu() + + # Check whether an offset slit is used. If so, recompute spectral traces. + # TODO: this needs to be removed. If transforms are available that + # can be applied to the actual slit position, they ought be used + # without modification. If they are not available but need to be + # hacked, then that hack needs to be defined via the yaml and only + # be used in a particular case, not in general. + #offset_x = obj.cube.header["CRVAL1D"] + #offset_y = obj.cube.header["CRVAL2D"] + #if (offset_x != self.meta["offset_x"] or + # offset_y != self.meta["offset_y"]): + # logger.debug("Recomputing spectral traces for offset (%.1g, %.1g)", + # offset_x, offset_y) + # self.meta["offset_x"] = offset_x + # self.meta["offset_y"] = offset_y + # self.make_spectral_traces() + # self.update_meta() + + spt = self.spectral_traces[obj.trace_id] + remove_pre = False + if ("offset_type" in obj.meta and + obj.meta["offset_type"] == "detector"): + self.meta["offset_type"] = "detector" + self.meta["offset_x"] = obj.meta["offset_x"] + self.meta["offset_y"] = obj.meta["offset_y"] + spt.xy2lam.pretransform_x = (det_offset, + {"offset": obj.meta["offset_x"]}) + spt.xy2lam.pretransform_y = (det_offset, + {"offset": obj.meta["offset_y"]}) + spt.xilam2x.posttransform = (det_offset, + {"offset": -obj.meta["offset_x"]}) + spt.xilam2y.posttransform = (det_offset, + {"offset": -obj.meta["offset_y"]}) + remove_pre = True + + obj.hdu = spt.map_spectra_to_focal_plane(obj) + + if remove_pre: + spt.xy2lam.pretransform_y = None + spt.xy2xi.pretransform_x = None + spt.xilam2x.posttransform = None + spt.xilam2y.posttransform = None + + return obj + @property def footprint(self): """Return the footprint of the entire SpectralTraceList.""" @@ -361,15 +410,36 @@ def rectify_traces(self, hdulist, xi_min=None, xi_max=None, interps=None, for i, trace_id in tqdm(enumerate(self.spectral_traces, start=1), desc=" Traces", total=len(self.spectral_traces)): - hdu = self[trace_id].rectify(hdulist, - interps=interps, - bin_width=bin_width, - xi_min=xi_min, xi_max=xi_max, - wave_min=wave_min, wave_max=wave_max) + spt = self[trace_id] + remove_pre = False + if ("offset_type" in self.meta and + self.meta["offset_type"] == "detector"): + spt.xy2lam.pretransform_x = (det_offset, + {"offset": self.meta["offset_x"]}) + spt.xy2lam.pretransform_y = (det_offset, + {"offset": self.meta["offset_y"]}) + spt.xilam2x.posttransform = (det_offset, + {"offset": -self.meta["offset_x"]}) + spt.xilam2y.posttransform = (det_offset, + {"offset": -self.meta["offset_y"]}) + remove_pre = True + + + hdu = spt.rectify(hdulist, + interps=interps, + bin_width=bin_width, + xi_min=xi_min, xi_max=xi_max, + wave_min=wave_min, wave_max=wave_max) if hdu is not None: # ..todo: rectify does not do that yet outhdul.append(hdu) outhdul[0].header[f"EXTNAME{i}"] = trace_id + if remove_pre: + spt.xy2lam.pretransform_y = None + spt.xy2xi.pretransform_x = None + spt.xilam2x.posttransform = None + spt.xilam2y.posttransform = None + outhdul[0].header.update(inhdul[0].header) return outhdul diff --git a/scopesim/effects/spectral_trace_list_utils.py b/scopesim/effects/spectral_trace_list_utils.py index 7541f9cb..1fc8e999 100644 --- a/scopesim/effects/spectral_trace_list_utils.py +++ b/scopesim/effects/spectral_trace_list_utils.py @@ -49,8 +49,6 @@ class SpectralTrace: "x_colname": "x", "y_colname": "y", "s_colname": "s", - "offset_x": 0, - "offset_y": 0, "wave_colname": "wavelength", "dwave": 0.002, "aperture_id": 0, @@ -114,6 +112,7 @@ def fov_grid(self): return {"wave_min": wave_min, "wave_max": wave_max, "trace_id": self.trace_id, "aperture_id": aperture_id} + def compute_interpolation_functions(self): """ Compute various interpolation functions between slit and focal plane. @@ -122,17 +121,21 @@ def compute_interpolation_functions(self): `xi` (spatial coordinate along the slit, in arcsec) and `lam` (wavelength, in um). - The interpolation functions include a shift in the focal-plane - coordinates, determined from the CRVAL of the source FOV. + #The interpolation functions include a shift in the focal-plane + #coordinates, determined from the CRVAL of the source FOV. """ - x_arr = self.table[self.meta["x_colname"]] + self.meta["offset_x"] - y_arr = self.table[self.meta["y_colname"]] + self.meta["offset_y"] + 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 + + # TODO There should be an option to include pre- and + # posttransforms. But how should they be defined? self.xy2xi = Transform2D.fit(x_arr, y_arr, xi_arr) self.xy2lam = Transform2D.fit(x_arr, y_arr, lam_arr) self.xilam2x = Transform2D.fit(xi_arr, lam_arr, x_arr) @@ -1042,3 +1045,7 @@ def get_affine_parameters(coords): shears = (np.average(shears, axis=0) * rad2deg) - (90 + rotations) return rotations, shears + +def det_offset(x, offset): + """Apply offset to x""" + return x + offset diff --git a/scopesim/optics/image_plane.py b/scopesim/optics/image_plane.py index 794c03fb..b58ac50c 100644 --- a/scopesim/optics/image_plane.py +++ b/scopesim/optics/image_plane.py @@ -68,7 +68,7 @@ def __init__(self, header, cmds=None, **kwargs): raise ValueError(f"header must have a valid image-plane WCS: " f"{dict(header)}") - image = zeros_from_header(header) + image = zeros_from_header(header, dtype=np.float32) self.hdu = fits.ImageHDU(data=image, header=header) self.hdu.header["BUNIT"] = "ph s-1" # photons per second (per pixel)