Skip to content

[07/29] Use calibrated traces and physical slit and detector coordinates #27

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
baileyji wants to merge 1 commit into
base: review/xs-06-echelle-geometry
Choose a base branch
from
Open

[07/29] Use calibrated traces and physical slit and detector coordinates #27

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 4 additions & 0 deletions scopesim/effects/apertures.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -128,6 +128,8 @@ def apply_to(self, obj, **kwargs):
for vol in obj.volumes:
vol["meta"]["xi_min"] = min(x) * u.arcsec
vol["meta"]["xi_max"] = max(x) * u.arcsec
vol["meta"]["eta_min"] = min(y) * u.arcsec
vol["meta"]["eta_max"] = max(y) * u.arcsec
vol["meta"]["slit_name"] = self.meta["name"]

if "offset_type" in self.meta:
Expand Down Expand Up @@ -305,6 +307,8 @@ def apply_to(self, obj, **kwargs):
# ..todo: HUGE HACK - Get rid of this!
vol["meta"]["xi_min"] = row["left"] * u.arcsec
vol["meta"]["xi_max"] = row["right"] * u.arcsec
vol["meta"]["eta_min"] = row["bottom"] * u.arcsec
vol["meta"]["eta_max"] = row["top"] * u.arcsec
vol["conserve_image"] = row["conserve_image"]
vol["shape"] = row["shape"]
vol["angle"] = row["angle"]
Expand Down
11 changes: 9 additions & 2 deletions scopesim/effects/illumination.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -538,6 +538,7 @@ class ImagePlaneBackground(Effect):
trace list. The supplied value or model should already include downstream
throughput and detector QE. For tapered QE maps, use an average or otherwise
representative positional QE when deriving this non-dispersed background.
A FITS ``filename`` and ``extname`` can supply an already registered map.
"""

z_order: ClassVar[tuple[int, ...]] = (760,)
Expand All @@ -562,15 +563,21 @@ def apply_to(self, obj, **kwargs):
return obj

shape = obj.hdu.data.shape
if self._map is None or shape != self._map_shape:
if self._file is not None or self._map is None or shape != self._map_shape:
self._map = self._make_map(shape)
self._map_shape = shape

obj.hdu.data = obj.hdu.data + self._map
return obj

def _make_map(self, shape):
if self._model is None:
if self._file is not None:
# The calibration image must be registered to this image plane.
extname = str(from_currsys(self.meta["extname"], self.cmds))
background = self._file[extname].data
if background.shape != shape:
raise ValueError("Background calibration shape differs from the image plane")
elif self._model is None:
value = from_currsys(self.meta["value"], self.cmds)
background = np.full(shape, value, dtype=np.float32)
else:
Expand Down
5 changes: 4 additions & 1 deletion scopesim/effects/metis_lms_trace_list.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -144,6 +144,8 @@ def apply_to_fov(self, obj):
slicefov.detector_header = obj.detector_header
slicefov.meta["xi_min"] = obj.meta["xi_min"]
slicefov.meta["xi_max"] = obj.meta["xi_max"]
slicefov.meta["eta_min"] = ymin * u.arcsec
slicefov.meta["eta_max"] = ymax * u.arcsec
slicefov.meta["trace_id"] = sptid
slicefov.cube = fits.ImageHDU(header=slicewcs.to_header(),
data=slicecube)
Expand Down Expand Up @@ -322,6 +324,7 @@ class MetisLMSSpectralTrace(SpectralTrace):
"""SpectralTrace for the METIS LM spectrograph."""

_class_params = {
**SpectralTrace._class_params,
"naxis1": 122,
"nslice": 28,
"slicewidth": 0.0207, # arcsec
Expand Down Expand Up @@ -393,7 +396,7 @@ def get_waverange(self, det_mm_lims):

lam0 = from_currsys(self.meta["wavelen"], self.cmds)
xi0 = 0.
ymid = self.xilam2y(xi0, lam0)[0] # estimate y level of trace
ymid = self.xilam2y(xi0, lam0) # scalar input gives a scalar trace coordinate
waverange = self.xy2lam(np.array([xmin, xmax]), np.array([ymid, ymid]),
grid=False)

Expand Down
180 changes: 107 additions & 73 deletions scopesim/effects/spectral_trace_list_utils.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,8 +14,8 @@

import numpy as np

from scipy.interpolate import RectBivariateSpline
from scipy.interpolate import interp1d
from scipy.interpolate import interp1d, RectBivariateSpline
from scipy.ndimage import gaussian_filter1d, map_coordinates

from astropy.table import Table, vstack
from astropy.modeling import fitting
Expand Down Expand Up @@ -103,6 +103,7 @@ class SpectralTrace:
"pixel_size": None,
"description": "<no description>",
"trace_flux_jacobian": "!SIM.spectral.trace_flux_jacobian",
"fit_degree": 4,
}

def __init__(self, trace_tbl, cmds=None, **kwargs):
Expand All @@ -128,6 +129,9 @@ def __init__(self, trace_tbl, cmds=None, **kwargs):
("PLTSCALE", "plate_scale"),
("DETPAD", "detector_pad"),
("DETANG", "detector_angle"),
("FITDEG", "fit_degree"),
("LSFRES", "optical_resolving_power"),
("SLITCOL", "slit_dispersion_colname"),
):
if header_key in trace_tbl.header:
self.meta[meta_key] = trace_tbl.header[header_key]
Expand Down Expand Up @@ -194,12 +198,19 @@ 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

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)
self.xilam2y = Transform2D.fit(xi_arr, lam_arr, y_arr)
self._xiy2x = Transform2D.fit(xi_arr, y_arr, x_arr)
self._xiy2lam = Transform2D.fit(xi_arr, y_arr, lam_arr)
degree = self.meta["fit_degree"]
self.xy2xi = Transform2D.fit(x_arr, y_arr, xi_arr, degree=degree)
self.xy2lam = Transform2D.fit(x_arr, y_arr, lam_arr, degree=degree)
self.xilam2x = Transform2D.fit(xi_arr, lam_arr, x_arr, degree=degree)
self.xilam2y = Transform2D.fit(xi_arr, lam_arr, y_arr, degree=degree)
self._xiy2x = Transform2D.fit(xi_arr, y_arr, x_arr, degree=degree)
self._xiy2lam = Transform2D.fit(xi_arr, y_arr, lam_arr, degree=degree)
if "slit_dispersion_colname" in self.meta:
center = xi_arr == 0
derivative = self.table[self.meta["slit_dispersion_colname"]].quantity.to_value(u.um / u.arcsec)
# FOVs are restricted to this trace; allow unit-conversion roundoff
# at the wavelength endpoints when interpolating the derivative.
self.slit_dispersion = interp1d(lam_arr[center], derivative[center], bounds_error=False, fill_value="extrapolate")

if self.dispersion_axis == "unknown":
dlam_dx, dlam_dy = self.xy2lam.gradient()
Expand Down Expand Up @@ -261,6 +272,24 @@ def _focal_plane_grid(self, fov):
y_mm = np.linspace(ymin_mm, ymax_mm, sub_naxis2, endpoint=False, dtype=np.float32)
return xmin, ymin, x_mm, y_mm, det_wcs

def spectral_padding(self, fov):
"""Input wavelength support for slit displacement and optical blur."""
lower, upper = fov.meta['wave_min'].to_value(u.um), fov.meta['wave_max'].to_value(u.um)
wave = np.unique(u.Quantity(self.table[self.meta['wave_colname']]).to_value(u.um))
wave = np.r_[lower, wave[(wave > lower) & (wave < upper)], upper]
width = (fov.meta['eta_max'] - fov.meta['eta_min']).to_value(u.arcsec)
if 'slit_dispersion_colname' in self.meta:
derivative = self.slit_dispersion(wave)
else:
pixsize = (fov.detector_header['CDELT1D']*u.Unit(fov.detector_header['CUNIT1D'])).to_value(u.mm)
self._set_dispersion(lower, upper, pixsize=pixsize)
step = abs(fov.header['CDELT2'])*3600
derivative = self.dlam_per_pix(wave)/step
padding = np.max(abs(derivative))*width/2
if 'optical_resolving_power' in self.meta:
padding += 5*upper/(self.meta['optical_resolving_power']*np.sqrt(8*np.log(2)))
return padding*u.um

def map_spectra_to_focal_plane(self, fov):
"""
Apply the spectral trace mapping to a spectral cube.
Expand Down Expand Up @@ -294,14 +323,13 @@ def map_spectra_to_focal_plane(self, fov):
image = np.zeros((sub_naxis2, sub_naxis1), dtype=np.float32)

self._set_dispersion(wave_min, wave_max, pixsize=pixsize)
try:
xilam = XiLamImage(fov, self.dlam_per_pix)
self._xilamimg = xilam # ..todo: remove or make available with a debug flag?
except ValueError:
logger.warning(" ---> %s gave ValueError", self.trace_id)

npix_xi, npix_lam = xilam.npix_xi, xilam.npix_lam
xilam_wcs = xilam.wcs
slit_dispersion = None
if "slit_dispersion_colname" in self.meta:
slit_dispersion = self.slit_dispersion
xilam = XiLamImage(fov, self.dlam_per_pix, slit_dispersion=slit_dispersion)
if "optical_resolving_power" in self.meta:
xilam.convolve_lsf(self.meta["optical_resolving_power"])
self._xilamimg = xilam

# focal-plane coordinate images
ximg_fpa, yimg_fpa = np.meshgrid(x_mm, y_mm)
Expand All @@ -310,33 +338,11 @@ def map_spectra_to_focal_plane(self, fov):
xi_fpa = self.xy2xi(ximg_fpa, yimg_fpa).astype(np.float32)
lam_fpa = self.xy2lam(ximg_fpa, yimg_fpa).astype(np.float32)

# mask everything outside the wavelength range
mask = (xi_fpa >= xi_min) & (xi_fpa <= xi_max)
xi_fpa *= mask
lam_fpa *= mask

# mask everything outside the wavelength range
# mask = (lam_fpa >= wave_min) & (lam_fpa <= wave_max)
# xi_fpa *= mask
# lam_fpa *= mask

# Convert to pixel images
# These are the pixel coordinates in the image corresponding to
# xi, lambda
# It is much quicker to do the linear transformation by hand
# than to use the astropy.wcs functions for conversion.
i_img = ((lam_fpa - xilam_wcs.wcs.crval[0])
/ xilam_wcs.wcs.cdelt[0]).astype(int)
j_img = ((xi_fpa - xilam_wcs.wcs.crval[1])
/ xilam_wcs.wcs.cdelt[1]).astype(int)

# truncate images to remove pixel coordinates outside the image
ijmask = ((i_img >= 0) * (i_img < npix_lam)
* (j_img >= 0) * (j_img < npix_xi))

# do the actual interpolation
# image is in [ph/s/um/arcsec]
image = xilam.interp(xi_fpa, lam_fpa, grid=False) * ijmask
# The physical slit ends include the outer half-cells. A truncated
# integer index incorrectly admits pixels below the lower boundary.
inside = (xi_fpa >= xi_min) & (xi_fpa < xi_max)
inside &= (lam_fpa >= wave_min) & (lam_fpa < wave_max)
image = xilam.interp(xi_fpa, lam_fpa, grid=False)*inside

# Scale to ph / s / pixel
image *= self._trace_flux_scale(
Expand Down Expand Up @@ -380,11 +386,7 @@ def map_spectra_to_focal_plane(self, fov):

def _trace_flux_scale(self, x_mm, y_mm, pixsize, pixscale, det_header):
"""Return the local [arcsec um] per detector pixel flux scale."""
try:
use_jacobian = bool(from_currsys(
self.meta["trace_flux_jacobian"], self.cmds))
except ValueError:
use_jacobian = False
use_jacobian = from_currsys(self.meta["trace_flux_jacobian"], self.cmds)

if not use_jacobian:
dlam_by_dx, dlam_by_dy = self.xy2lam.gradient()
Expand All @@ -401,13 +403,10 @@ def _trace_flux_scale(self, x_mm, y_mm, pixsize, pixscale, det_header):
@staticmethod
def _detector_pixel_area(det_header):
"""Return detector WCS pixel area in mm2."""
try:
pixel_scale_matrix = WCS(det_header, key="D").pixel_scale_matrix
return np.abs(np.linalg.det(pixel_scale_matrix))
except Exception: # pragma: no cover - FITS fallback for odd headers
cdelt1 = det_header["CDELT1D"] * u.Unit(det_header["CUNIT1D"])
cdelt2 = det_header["CDELT2D"] * u.Unit(det_header["CUNIT2D"])
return np.abs((cdelt1 * cdelt2).to_value(u.mm**2))
detector_wcs = WCS(det_header, key="D")
pixel_area = np.abs(np.linalg.det(detector_wcs.pixel_scale_matrix))
pixel_area *= u.Unit(detector_wcs.wcs.cunit[0]) * u.Unit(detector_wcs.wcs.cunit[1])
return pixel_area.to_value(u.mm**2)

def rectify(
self,
Expand Down Expand Up @@ -851,9 +850,13 @@ class XiLamImage():
dlam_per_pix : a 1-D interpolation function from wavelength (in um) to
dispersion (in um/pixel); alternatively a number giving an average
dispersion
slit_dispersion : callable, optional
Physical cross-slit derivative d(lambda)/d(eta) in um/arcsec as a
function of wavelength in um. When supplied, this replaces the
legacy coupling between sky-grid sampling and detector dispersion.
"""

def __init__(self, fov, dlam_per_pix):
def __init__(self, fov, dlam_per_pix, slit_dispersion=None):
# ..todo: we assume that we always have a cube. We use SpecCADO's
# add_cube_layer method
cube_wcs = WCS(fov.cube.header, key=" ")
Expand All @@ -871,8 +874,16 @@ def __init__(self, fov, dlam_per_pix):
(n_lam, n_eta, n_xi) = fov.cube.data.shape

# arrays of cube coordinates
cube_xi = d_xi * np.arange(n_xi) + fov.meta["xi_min"].value
cube_eta = d_eta * (np.arange(n_eta) - (n_eta - 1) / 2)
cube_xi = cube_wcs.all_pix2world(np.arange(n_xi), np.zeros(n_xi), np.zeros(n_xi), 0)[0]
cube_eta = cube_wcs.all_pix2world(np.zeros(n_eta), np.arange(n_eta), np.zeros(n_eta), 0)[1]
if cube_wcs.has_celestial:
# Slit coordinates are offsets about zero; RA wraps at 360 degrees.
cube_xi = (cube_xi + 180) % 360 - 180
cube_xi *= u.Unit(cube_wcs.wcs.cunit[0]).to(u.arcsec)
cube_eta *= u.Unit(cube_wcs.wcs.cunit[1]).to(u.arcsec)
# Each IFU slice has its own optical axis, independent of sky position.
eta_center = (fov.meta["eta_min"] + fov.meta["eta_max"]) / 2
cube_eta -= eta_center.to_value(u.arcsec)
cube_lam = wcs_lam.all_pix2world(np.arange(n_lam), 0)[0]
cube_lam *= u.Unit(wcs_lam.wcs.cunit[0]).to(u.um)

Expand All @@ -886,24 +897,31 @@ def __init__(self, fov, dlam_per_pix):
logger.warning("Using scalar dlam_per_pix = %.2g",
dlam_per_pix_val)

dlam_per_arcsec = dlam_per_pix_val / d_eta if slit_dispersion is None else slit_dispersion(self.lam)
for i, eta in enumerate(cube_eta):
lam0 = self.lam + dlam_per_pix_val * eta / d_eta
# A positive calibrated derivative moves an input line to a larger
# apparent wavelength: sample the source at lambda_out - shift.
lam0 = self.lam + dlam_per_arcsec * eta if slit_dispersion is None else self.lam - dlam_per_arcsec * eta

# lam0 is the target wavelength. We need to check that this
# overlaps with the wavelength range covered by the cube
if lam0.min() < cube_lam.max() and lam0.max() > cube_lam.min():
plane = fov.cube.data[:, i, :].T
plane_interp = RectBivariateSpline(cube_xi, cube_lam, plane,
kx=1, ky=1)
self.image += plane_interp(cube_xi, lam0)
# Only wavelength changes here; xi is already on its native
# grid. Building a 2D spline for every slit cell is redundant.
index = np.searchsorted(cube_lam, lam0, side='right') - 1
index = np.clip(index, 0, n_lam - 2)
fraction = (lam0 - cube_lam[index])/d_lam
valid = (lam0 >= cube_lam[0]) & (lam0 <= cube_lam[-1])
plane = fov.cube.data[:, i, :]
self.image += ((1 - fraction[:, None])*plane[index] + fraction[:, None]*plane[index + 1]).T*valid

self.image *= d_eta # ph/s/um/arcsec2 --> ph/s/um/arcsec
self.image *= abs(d_eta) # ph/s/um/arcsec2 --> ph/s/um/arcsec

# WCS for the xi-lambda image, i.e. the rectified 2D spectrum
# Default WCS with xi in arcsec
self.wcs = WCS(naxis=2)
self.wcs.wcs.crpix = [1, 1]
self.wcs.wcs.crval = [self.lam[0], fov.meta["xi_min"].value]
self.wcs.wcs.crval = [self.lam[0], cube_xi[0]]
self.wcs.wcs.pc = [[1, 0], [0, 1]]
self.wcs.wcs.cdelt = [d_lam, d_xi]
self.wcs.wcs.ctype = ["LINEAR", "LINEAR"]
Expand All @@ -923,14 +941,30 @@ def __init__(self, fov, dlam_per_pix):
self.xi = self.wcs.all_pix2world(self.lam[0], np.arange(n_xi), 0)[1]
self.npix_xi = n_xi
self.npix_lam = n_lam
# ..todo: cubic spline introduces negative values, linear does not.
# Alternative might be to cubic-spline interpolate on sqrt(image),
# with subsequent squaring of the result. This would require
# wrapping RectBivariateSpline in a new (sub)class.
spline_order = (1, 1)
self.interp = RectBivariateSpline(self.xi, self.lam, self.image,
kx=spline_order[0],
ky=spline_order[1])

def interp(self, xi, wavelength, grid=True):
"""Bilinear density interpolation on the existing uniform WCS grid."""
xi, wavelength = np.meshgrid(xi, wavelength, indexing='ij') if grid else np.broadcast_arrays(xi, wavelength)
row = (xi - self.xi[0])/self.wcs.wcs.cdelt[1]
column = (wavelength - self.lam[0])/self.wcs.wcs.cdelt[0]
return map_coordinates(self.image, [row, column], order=1, mode='nearest', prefilter=False)

def convolve_lsf(self, resolving_power):
"""Apply an optical Gaussian LSF at constant R, preserving photon density.

This is the optical blur alone; slit illumination is integrated above.
Convolving lambda * I_lambda on a uniform log(lambda) grid preserves
integrated flux rather than treating I_lambda as a log-bin count.
"""
if resolving_power <= 0:
raise ValueError("Optical resolving power must be positive")
log_wave = np.log(self.lam)
n = int(np.ceil((log_wave[-1] - log_wave[0]) / np.min(np.diff(log_wave)))) + 1
log_grid = np.linspace(log_wave[0], log_wave[-1], n)
sigma = 1 / (resolving_power * np.sqrt(8 * np.log(2)) * (log_grid[1] - log_grid[0]))
density = interp1d(log_wave, self.image * self.lam, axis=1)(log_grid)
blurred = gaussian_filter1d(density, sigma, axis=1, mode="nearest")
self.image = interp1d(log_grid, blurred, axis=1)(log_wave) / self.lam


class Transform2D():
Expand Down
7 changes: 6 additions & 1 deletion scopesim/effects/ter_curves.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -547,9 +547,14 @@ def load_skycalc_table(self):

else:
path = find_file(use_local_file)
fits_tbl = fits.getdata(path, ext=1)
with fits.open(path) as sky_file:
fits_tbl = sky_file[1].data
# Older cached files require an explicit unit if TUNIT is absent.
unit = self.meta["local_skycalc_wavelength_unit"] if "local_skycalc_wavelength_unit" in self.meta else sky_file[1].columns["lam"].unit
wave_unit = u.Unit(unit)
fits_hdr = fits.getheader(path, ext=0)
tbl = Table(fits_tbl)
tbl["lam"] = tbl["lam"] * wave_unit.to(u.um)
tbl["lam"].unit = u.um
for colname in tbl.colnames:
if "flux" in colname:
Expand Down
9 changes: 9 additions & 0 deletions scopesim/tests/test_utils_functions.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -199,6 +199,15 @@ def test_converts_list(self):
def test_converts_numpy_array(self):
assert utils.from_currsys(np.array(["!SIM.random.seed"]*2))[1] is None

@pytest.mark.parametrize("dtype", [np.float32, np.float64, np.int32])
def test_numeric_calibration_array_keeps_values_and_precision(self, dtype):
values = np.arange(12, dtype=dtype).reshape(3, 4)
resolved = utils.from_currsys(values)
np.testing.assert_array_equal(resolved, values)
assert resolved.dtype == dtype
resolved[0, 0] = 20
assert values[0, 0] == 0

def test_converts_dict(self):
assert utils.from_currsys({"seed": "!SIM.random.seed"})["seed"] is None

Expand Down
Loading