diff --git a/scopesim/detector/detector_array.py b/scopesim/detector/detector_array.py index c073b27b..b53df817 100644 --- a/scopesim/detector/detector_array.py +++ b/scopesim/detector/detector_array.py @@ -19,7 +19,7 @@ def __init__(self, detector_list=None, **kwargs): self.detectors = [] self.latest_exposure = None - def readout(self, image_planes, array_effects=[], dtcr_effects=[], **kwargs): + def readout(self, image_planes, array_effects=None, dtcr_effects=None, **kwargs): """ Read out the detector array into a FITS file @@ -54,8 +54,8 @@ def readout(self, image_planes, array_effects=[], dtcr_effects=[], **kwargs): # - add ImageHDUs # - add ASCIITableHDU with Effects meta data in final table extension - self.array_effects = array_effects - self.dtcr_effects = dtcr_effects + self.array_effects = array_effects or [] + self.dtcr_effects = dtcr_effects or [] self.meta.update(kwargs) # 0. Get the image plane that corresponds to this detector array diff --git a/scopesim/optics/fov.py b/scopesim/optics/fov.py index c5e6c688..6a87d70f 100644 --- a/scopesim/optics/fov.py +++ b/scopesim/optics/fov.py @@ -1,5 +1,7 @@ """Defines FieldOfView class""" +import logging from copy import deepcopy +from itertools import chain import numpy as np from scipy.interpolate import interp1d @@ -83,6 +85,10 @@ def __init__(self, header, waverange, detector_header=None, **kwargs): self._wavelength = None self._volume = None + def _pixarea(self, hdr): + return (hdr["CDELT1"] * u.Unit(hdr["CUNIT1"]) * + hdr["CDELT2"] * u.Unit(hdr["CUNIT2"])).to(u.arcsec ** 2) + @property def trace_id(self): """Return the name of the trace""" @@ -91,11 +97,7 @@ def trace_id(self): @property def pixel_area(self): if self.meta["pixel_area"] is None: - hdr = self.header - pixarea = (hdr["CDELT1"] * u.Unit(hdr["CUNIT1"]) * - hdr["CDELT2"] * u.Unit(hdr["CUNIT2"])).to(u.arcsec ** 2) - self.meta["pixel_area"] = pixarea.value # [arcsec] - + self.meta["pixel_area"] = self._pixarea(self.header).value # [arcsec] return self.meta["pixel_area"] def extract_from(self, src): @@ -111,13 +113,16 @@ def extract_from(self, src): fields_in_fov = [field for field in src.fields if fu.is_field_in_fov(self.header, field)] + if not fields_in_fov: + logging.warning("No fields in FOV.") - spec_refs = [] + spec_refs = set() volume = self.volume() for ifld, fld in enumerate(fields_in_fov): if isinstance(fld, Table): - fields_in_fov[ifld] = fu.extract_area_from_table(fld, volume) - spec_refs += list(np.unique(fields_in_fov[ifld] ["ref"])) + extracted_field = fu.extract_area_from_table(fld, volume) + spec_refs.update(extracted_field["ref"]) + fields_in_fov[ifld] = extracted_field elif isinstance(fld, fits.ImageHDU): if fld.header["NAXIS"] in (2, 3): @@ -125,11 +130,11 @@ def extract_from(self, src): if fld.header["NAXIS"] == 2 or fld.header.get("BG_SRC"): ref = fld.header.get("SPEC_REF") if ref is not None: - spec_refs += [ref] + spec_refs.add(ref) waves = volume["waves"] * u.Unit(volume["wave_unit"]) spectra = {ref: fu.extract_range_from_spectrum(src.spectra[ref], waves) - for ref in np.unique(spec_refs)} + for ref in spec_refs} self.fields = fields_in_fov self.spectra = spectra @@ -153,6 +158,7 @@ def view(self, hdu_type="image", sub_pixel=None, use_photlam=None): self.meta["sub_pixel"] = sub_pixel if hdu_type == "image": + # FIXME: why not just make False the default value?? use_photlam = False if use_photlam is None else use_photlam self.hdu = self.make_image_hdu(use_photlam=use_photlam) elif hdu_type == "cube": @@ -163,82 +169,163 @@ def view(self, hdu_type="image", sub_pixel=None, use_photlam=None): return self.hdu def flatten(self): + """If cube, collapse along first axis.""" if self.hdu and self.hdu.header["NAXIS"] == 3: image = np.sum(self.hdu.data, axis=0) self.hdu.data = image - def make_spectrum(self): - """ - This is needed for when we do incoherent MOS instruments. - Each fibre doesn't care about the spatial information. - - 1. Make waveset and zeros flux - make dict of spectra evaluated at waveset + def _evaluate_spectrum_with_weight(self, ref, waveset, weight): + return self.spectra[ref](waveset).value * weight - 2. Find Cube fields - collapse cube along spatial dimensions --> spectrum vector - convert vector to PHOTLAM - interpolate at waveset - add to zeros flux vector - - 3. Find Image fields - sum image over both dimensions - evaluate SPEC_REF spectum at waveset - multiply by sum - add to zeros flux vector - - 4. Find Table fields - evaluate all spectra at waveset - for each unique ref, sum the weights - add each spectra * sum of weights to zeros flux vector - - Returns - ------- - spec : SourceSpectrum - [PHOTLAM] + def _calc_area_factor(self, field): + bg_solid_angle = u.Unit(field.header["SOLIDANG"]).to(u.arcsec**-2) + area_factor = self.pixel_area * bg_solid_angle # arcsec**2 * arcsec**-2 + return area_factor + def _make_spectrum_cubefields(self, fov_waveset): """ - fov_waveset = self.waveset - canvas_flux = np.zeros(len(fov_waveset)) - - specs = {ref: spec(fov_waveset) for ref, spec in self.spectra.items()} + Find Cube fields. + * collapse cube along spatial dimensions --> spectrum vector + * convert vector to PHOTLAM + * interpolate at waveset + * yield scaled flux to be added to canvas flux + """ for field in self.cube_fields: - hdu_waveset = fu.get_cube_waveset(field.header, return_quantity=True) + hdu_waveset = fu.get_cube_waveset(field.header, + return_quantity=True) fluxes = field.data.sum(axis=2).sum(axis=1) fov_waveset_fluxes = np.interp(fov_waveset, hdu_waveset, fluxes) field_unit = field.header.get("BUNIT", PHOTLAM) flux_scale_factor = u.Unit(field_unit).to(PHOTLAM) - canvas_flux += fov_waveset_fluxes * flux_scale_factor + yield fov_waveset_fluxes * flux_scale_factor + def _make_spectrum_imagefields(self, waveset): + """ + Find Image fields. + + * sum image over both dimensions + * evaluate SPEC_REF spectum at waveset + * yield spectrum multiply by sum to be added to canvas flux + """ for field in self.image_fields: - ref = field.header["SPEC_REF"] weight = np.sum(field.data) - canvas_flux += self.spectra[ref](fov_waveset).value * weight + yield self._evaluate_spectrum_with_weight(field.header["SPEC_REF"], + waveset, weight) + def _make_spectrum_tablefields(self, waveset): + """ + Find Table fields. + * evaluate all spectra at waveset + * for each unique ref, sum the weights + * yield each spectrum * sum of weights to be added to canvas flux + """ for field in self.table_fields: refs = np.array(field["ref"]) weights = np.array(field["weight"]) - weight_sums = {ref: np.sum(weights[refs == ref]) - for ref in np.unique(refs)} - for ref, weight in weight_sums.items(): - canvas_flux += self.spectra[ref](fov_waveset).value * weight + # TODO: could do grouping of table with both columns?? + for ref in set(refs): + weight = np.sum(weights, where=refs==ref) + yield self._evaluate_spectrum_with_weight(ref, waveset, weight) + def _make_spectrum_backfields(self, waveset): for field in self.background_fields: - bg_solid_angle = u.Unit(field.header["SOLIDANG"]).to(u.arcsec**-2) - area_factor = self.pixel_area * bg_solid_angle # arcsec**2 * arcsec**-2 + yield self._evaluate_spectrum_with_weight( + field.header["SPEC_REF"], waveset, + self._calc_area_factor(field)) + + def make_spectrum(self): + """ + This is needed for when we do incoherent MOS instruments. + Each fibre doesn't care about the spatial information. - ref = field.header["SPEC_REF"] - canvas_flux += self.spectra[ref](fov_waveset).value * area_factor + Returns + ------- + spec : SourceSpectrum + [PHOTLAM] + + """ + fov_waveset = self.waveset + # Start with zero flux no ensure correct array shape even if none of + # the sub-functions yield anything. + canvas_flux = sum(chain( + self._make_spectrum_cubefields(fov_waveset), + self._make_spectrum_imagefields(fov_waveset), + self._make_spectrum_tablefields(fov_waveset), + self._make_spectrum_backfields(fov_waveset), + ), start=np.zeros_like(fov_waveset.value)) spectrum = SourceSpectrum(Empirical1D, points=fov_waveset, lookup_table=canvas_flux) - return spectrum + def _make_image_cubefields(self, area): + """ + Find Cube fields. + + * collapse cube along wavelength axis + * rescale and reorient image + * yield cube image to be added to canvas image + """ + for field in self.cube_fields: + # cube_fields come in with units of photlam/arcsec2, need to convert to ph/s + # We need to the voxel volume (spectral and solid angle) for that. + # ..todo: implement branch for use_photlam is True + spectral_bin_width = (field.header["CDELT3"] * + u.Unit(field.header["CUNIT3"])).to(u.Angstrom) + # First collapse to image, then convert units + image = np.sum(field.data, axis=0) * PHOTLAM/u.arcsec**2 + image = (image * self._pixarea(field.header) * area * + spectral_bin_width).to(u.ph/u.s) + yield fits.ImageHDU(data=image, header=field.header) + + def _make_image_imagefields(self, fluxes): + """ + Find Image fields. + + * sum spectra between wavelength edges + * multiply image by summed spectrum + * yield image to be added to canvas image + """ + for field in self.image_fields: + image = deepcopy(field.data) + image *= fluxes[field.header["SPEC_REF"]] # ph / s + yield fits.ImageHDU(data=image, header=field.header) + + def _make_image_tablefields(self, fluxes): + """ + Find Table fields. + + * sum spectra between wavelength edges + * yield summed flux at x,y position to be added to canvas image + """ + for field in self.table_fields: + # x, y are ALWAYS in arcsec - crval is in deg + xpix, ypix = imp_utils.val2pix(self.header, + field["x"] / 3600, + field["y"] / 3600) + if utils.from_currsys(self.meta["sub_pixel"]): + for idx, row in enumerate(field): + xs, ys, fracs = imp_utils.sub_pixel_fractions(xpix[idx], + ypix[idx]) + for x, y, frac in zip(xs, ys, fracs): + yield fluxes[row["ref"]] * frac, row["weight"], x, y + else: + # Note: these had x/ypix+0.5 until a06ab75 + # TODO: could these be something more numpythonic grid-ish? + x = np.array(xpix).astype(int) + y = np.array(ypix).astype(int) # quickest way to round + flux = np.array([fluxes[ref] for ref in field["ref"]]) + yield flux, np.array(field["weight"]), x, y + + def _make_image_backfields(self, fluxes): + for field in self.background_fields: + flux = fluxes[field.header["SPEC_REF"]] + yield flux * self._calc_area_factor(field) + def make_image_hdu(self, use_photlam=False): """ Used for imaging @@ -249,23 +336,7 @@ def make_image_hdu(self, use_photlam=False): See make_cube for an explanation - 1. Make waveset and canvas image - make canvas image from NAXIS1,2 from fov.header - - 2. Find Cube fields - collapse cube along wavelength axis - rescale and reorient image - add cube image to canvas image - - 3. Find Image fields - rescale and reorient images - sum spectra between wavelength edges - multiply image by summed spectrum - add image to canvas image - - 4. Find Table fields - sum spectra between wavelength edges - add summed flux at x,y position in canvas image + Make canvas image from NAXIS1,2 from fov.header Parameters ---------- @@ -280,91 +351,131 @@ def make_image_hdu(self, use_photlam=False): """ spline_order = utils.from_currsys("!SIM.computing.spline_order") - # 1. Make waveset and canvas image - fov_waveset = np.unique(self.waveset) + # Make waveset and canvas image + fov_waveset = self.waveset bin_widths = np.diff(fov_waveset) # u.um bin_widths = 0.5 * (np.r_[0, bin_widths] + np.r_[bin_widths, 0]) area = utils.from_currsys(self.meta["area"]) # u.m2 # PHOTLAM * u.um * u.m2 --> ph / s - specs = {ref: spec(fov_waveset) for ref, spec in self.spectra.items()} - if use_photlam is False: - for key in specs: - specs[key] = (specs[key] * bin_widths * area).to(u.ph / u.s) - - fluxes = {ref: np.sum(spec).value for ref, spec in specs.items()} - naxis1, naxis2 = self.header["NAXIS1"], self.header["NAXIS2"] - canvas_image_hdu = fits.ImageHDU(data=np.zeros((naxis2, naxis1)), - header=self.header) - - # 2. Find Cube fields - for field in self.cube_fields: - # cube_fields come in with units of photlam/arcsec2, need to convert to ph/s - # We need to the voxel volume (spectral and solid angle) for that. - # ..todo: implement branch for use_photlam is True - spectral_bin_width = (field.header["CDELT3"] * - u.Unit(field.header["CUNIT3"])).to(u.Angstrom) - pixarea = (field.header["CDELT1"] * u.Unit(field.header["CUNIT1"]) * - field.header["CDELT2"] * u.Unit(field.header["CUNIT2"])).to(u.arcsec**2) - - # First collapse to image, then convert units - image = np.sum(field.data, axis=0) * PHOTLAM/u.arcsec**2 - image = (image * pixarea * area * spectral_bin_width).to(u.ph/u.s) + specs = {ref: spec(fov_waveset) if use_photlam + else (spec(fov_waveset) * bin_widths * area).to(u.ph / u.s) + for ref, spec in self.spectra.items()} + fluxes = {ref: np.sum(spec.value) for ref, spec in specs.items()} + canvas_image_hdu = fits.ImageHDU( + data=np.zeros((self.header["NAXIS2"], self.header["NAXIS1"])), + header=self.header) - tmp_hdu = fits.ImageHDU(data=image, header=field.header) + for tmp_hdu in chain(self._make_image_cubefields(area), + self._make_image_imagefields(fluxes)): canvas_image_hdu = imp_utils.add_imagehdu_to_imagehdu( tmp_hdu, canvas_image_hdu, conserve_flux=True, spline_order=spline_order) - # 2. Find Image fields + for flux, weight, x, y in self._make_image_tablefields(fluxes): + if utils.from_currsys(self.meta["sub_pixel"]): + canvas_image_hdu.data[y, x] += flux * weight + else: + # Mask out any stars that were pushed out of the fov by rounding + mask = ((x < canvas_image_hdu.data.shape[1]) * + (y < canvas_image_hdu.data.shape[0])) + canvas_image_hdu.data[y[mask], x[mask]] += flux[mask] * weight[mask] + + canvas_image_hdu.data = sum(self._make_image_backfields(fluxes), + start=canvas_image_hdu.data) + return canvas_image_hdu # [ph s-1 pixel-1] + + def _make_cube_cubefields(self, fov_waveset): + """ + Find Cube fields. + + * rescale and reorient cubes + * interp1d smaller cubes with waveset + * yield cubes to be added to cavas cube + """ + for field in self.cube_fields: + # Cube should be in PHOTLAM arcsec-2 for SpectralTrace mapping + # Assumption is that ImageHDUs have units of PHOTLAM arcsec-2 + field_waveset = fu.get_cube_waveset(field.header, + return_quantity=True) + + # ..todo: Deal with this bounds_error in a more elegant way + field_interp = interp1d(field_waveset.to(u.um).value, + field.data, axis=0, kind="linear", + bounds_error=False, fill_value=0) + + field_data = field_interp(fov_waveset.value) + + # Pixel scale conversion + field_data *= self._pixarea(field.header).value / self.pixel_area + field_hdu = fits.ImageHDU(data=field_data, header=field.header) + yield field_hdu + + def _make_cube_imagefields(self, specs, spline_order): + """ + Find Image fields. + + * rescale and reorient images + * evaluate spectra at waveset + * expand image by spectra to 3D form + * yield image cubes to be added to cavas cube + """ for field in self.image_fields: - image = deepcopy(field.data) - tmp_hdu = fits.ImageHDU(data=image, header=field.header) - tmp_hdu.data *= fluxes[field.header["SPEC_REF"]] # ph / s + # Cube should be in PHOTLAM arcsec-2 for SpectralTrace mapping + # Assumption is that ImageHDUs have units of PHOTLAM arcsec-2 + # ImageHDUs have photons/second/pixel. + # ..todo: Add a catch to get ImageHDU with BUNITs + canvas_image_hdu = fits.ImageHDU( + data=np.zeros((self.header["NAXIS2"], self.header["NAXIS1"])), + header=self.header) + # FIXME: Why here not "Pixel scale conversion" as above? + field.data /= self.pixel_area canvas_image_hdu = imp_utils.add_imagehdu_to_imagehdu( - tmp_hdu, + field, canvas_image_hdu, - conserve_flux=True, spline_order=spline_order) + spec = specs[field.header["SPEC_REF"]] + field_cube = canvas_image_hdu.data[None, :, :] * spec[:, None, None] # 2D * 1D -> 3D + yield field_cube.value + + def _make_cube_tablefields(self, specs): + """ + Find Table fields. - # 3. Find Table fields + * evaluate spectra at waveset + * yield spectrum at x,y position to be added cavas to cube + """ for field in self.table_fields: - # x, y are ALWAYS in arcsec - crval is in deg - xpix, ypix = imp_utils.val2pix(self.header, - field["x"] / 3600, - field["y"] / 3600) - if utils.from_currsys(self.meta["sub_pixel"]): - for i, row in enumerate(field): - xs, ys, fracs = imp_utils.sub_pixel_fractions(xpix[i], - ypix[i]) - ref, weight = row["ref"], row["weight"] - for x, y, f in zip(xs, ys, fracs): - canvas_image_hdu.data[y, x] += fluxes[ref] * weight * f - else: - # Note: these had x/ypix+0.5 until a06ab75 - x = np.array(xpix).astype(int) - y = np.array(ypix).astype(int) # quickest way to round - f = np.array([fluxes[ref] for ref in field["ref"]]) - weight = np.array(field["weight"]) + # Cube should be in PHOTLAM arcsec-2 for SpectralTrace mapping + # Point sources are in PHOTLAM per pixel + # Point sources need to be scaled up by inverse pixel_area + for row in field: + xsky, ysky = row["x"] / 3600, row["y"] / 3600 + # x, y are ALWAYS in arcsec - crval is in deg + xpix, ypix = imp_utils.val2pix(self.header, xsky, ysky) + flux_vector = specs[row["ref"]].value * row["weight"] / self.pixel_area - # Mask out any stars that were pushed out of the fov by rounding - m = (x < canvas_image_hdu.data.shape[1]) * \ - (y < canvas_image_hdu.data.shape[0]) - canvas_image_hdu.data[y[m], x[m]] += f[m] * weight[m] + if utils.from_currsys(self.meta["sub_pixel"]): + xs, ys, fracs = imp_utils.sub_pixel_fractions(xpix, ypix) + for i, j, k in zip(xs, ys, fracs): + yield flux_vector * k, i, j + else: + yield flux_vector, int(xpix), int(ypix) - # 4. Find Background fields + def _make_cube_backfields(self, specs): for field in self.background_fields: - bg_solid_angle = u.Unit(field.header["SOLIDANG"]).to(u.arcsec**-2) - area_factor = self.pixel_area * bg_solid_angle # arcsec**2 * arcsec**-2 - - flux = fluxes[field.header["SPEC_REF"]] * area_factor - canvas_image_hdu.data += flux - - image_hdu = canvas_image_hdu # [ph s-1 pixel-1] + # TODO: The following would have been identical to the other two + # make methods, but was commented out. Why? + # bg_solid_angle = u.Unit(field.header["SOLIDANG"]).to(u.arcsec**-2) # float [arcsec-2] + # pixel_area = utils.from_currsys(self.meta["pixel_scale"]) ** 2 # float [arcsec2] + # area_factor = pixel_area * bg_solid_angle # float [arcsec2 * arcsec-2] - return image_hdu + # Cube should be in PHOTLAM arcsec-2 for SpectralTrace mapping + # spec = specs[field.header["SPEC_REF"]] * area_factor + spec = specs[field.header["SPEC_REF"]] + yield spec[:, None, None].value def make_cube_hdu(self): """ @@ -391,23 +502,9 @@ def make_cube_hdu(self): make waveset from self.meta values make canvas cube based on waveset of largest cube and NAXIS1,2 from fov.header - 2. Find Cube fields:: - - rescale and reorient cubes - interp1d smaller cubes with waveset - add cubes to cavas cube - - 3. Find Image fields:: - - rescale and reorient images - evaluate spectra at waveset - expand image by spectra to 3D form - add image cubes to canvas cube - - 4. Find Table fields:: - - evaluate spectra at waveset - add spectrum at x,y position in canvas cube + 2. Find Cube fields (see ``FieldOfView._make_cube_cubefields()``). + 3. Find Image fields (see ``FieldOfView._make_cube_imagefields()``). + 4. Find Table fields (see ``FieldOfView._make_cube_tablefields()``). ``PHOTLAM = ph/s/m2/um``. Original source fields are in units of: @@ -428,15 +525,14 @@ def make_cube_hdu(self): # 1. Make waveset and canvas cube (area, bin_width are applied at end) - - wave_min = self.meta["wave_min"] # Quantity [um] - wave_max = self.meta["wave_max"] - + # TODO: Why is this not self.waveset? What's different? wave_unit = u.Unit(utils.from_currsys("!SIM.spectral.wave_unit")) - dwave = utils.from_currsys("!SIM.spectral.spectral_bin_width") # Not a quantity - fov_waveset = np.arange(wave_min.value, wave_max.value, dwave) * wave_unit + fov_waveset = np.arange( + self.meta["wave_min"].value, self.meta["wave_max"].value, + utils.from_currsys("!SIM.spectral.spectral_bin_width")) * wave_unit fov_waveset = fov_waveset.to(u.um) + # TODO: what's with this code?? # fov_waveset = self.waveset # wave_bin_n = len(fov_waveset) # if "lin" in self.meta["wave_bin_type"]: @@ -449,90 +545,30 @@ def make_cube_hdu(self): for ref, spec in self.spectra.items()} # make canvas cube based on waveset of largest cube and NAXIS1,2 from fov.header - naxis1, naxis2 = self.header["NAXIS1"], self.header["NAXIS2"] - naxis3 = len(fov_waveset) - canvas_cube_hdu = fits.ImageHDU(data=np.zeros((naxis3, naxis2, naxis1)), - header=self.header) + canvas_cube_hdu = fits.ImageHDU( + data=np.zeros((len(fov_waveset), + self.header["NAXIS2"], + self.header["NAXIS1"])), + header=self.header) canvas_cube_hdu.header["BUNIT"] = "ph s-1 cm-2 AA-1" - # 2. Add Cube fields - for field in self.cube_fields: - # Cube should be in PHOTLAM arcsec-2 for SpectralTrace mapping - # Assumption is that ImageHDUs have units of PHOTLAM arcsec-2 - field_waveset = fu.get_cube_waveset(field.header, - return_quantity=True) - - # ..todo: Deal with this bounds_error in a more elegant way - field_interp = interp1d(field_waveset.to(u.um).value, - field.data, axis=0, kind="linear", - bounds_error=False, fill_value=0) - - field_data = field_interp(fov_waveset.value) - - # Pixel scale conversion - field_pixarea = (field.header["CDELT1"] - * field.header["CDELT2"] - * u.Unit(field.header["CUNIT1"]) - * u.Unit(field.header["CUNIT2"])).to(u.arcsec**2) - field_pixarea = field_pixarea.value - field_data *= field_pixarea / self.pixel_area - field_hdu = fits.ImageHDU(data=field_data, header=field.header) + for field_hdu in self._make_cube_cubefields(fov_waveset): canvas_cube_hdu = imp_utils.add_imagehdu_to_imagehdu( field_hdu, canvas_cube_hdu, spline_order=spline_order) - # 3. Find Image fields - for field in self.image_fields: - # Cube should be in PHOTLAM arcsec-2 for SpectralTrace mapping - # Assumption is that ImageHDUs have units of PHOTLAM arcsec-2 - # ImageHDUs have photons/second/pixel. - # ..todo: Add a catch to get ImageHDU with BUNITs - canvas_image_hdu = fits.ImageHDU(data=np.zeros((naxis2, naxis1)), - header=self.header) - pixarea = (field.header["CDELT1"] * u.Unit(field.header["CUNIT1"]) * - field.header["CDELT2"] * u.Unit(field.header["CUNIT2"])).to(u.arcsec**2) - - field.data = field.data / self.pixel_area - canvas_image_hdu = imp_utils.add_imagehdu_to_imagehdu(field, - canvas_image_hdu, - spline_order=spline_order) - spec = specs[field.header["SPEC_REF"]] - field_cube = canvas_image_hdu.data[None, :, :] * spec[:, None, None] # 2D * 1D -> 3D - canvas_cube_hdu.data += field_cube.value + canvas_cube_hdu.data = sum(self._make_cube_imagefields(specs, + spline_order), + start=canvas_cube_hdu.data) - # 4. Find Table fields - for field in self.table_fields: - # Cube should be in PHOTLAM arcsec-2 for SpectralTrace mapping - # Point sources are in PHOTLAM per pixel - # Point sources need to be scaled up by inverse pixel_area - pixel_area = self.pixel_area - for row in field: - xsky, ysky = row["x"], row["y"] - ref, weight = row["ref"], row["weight"] - # x, y are ALWAYS in arcsec - crval is in deg - xpix, ypix = imp_utils.val2pix(self.header, xsky / 3600, ysky / 3600) - if utils.from_currsys(self.meta["sub_pixel"]): - xs, ys, fracs = imp_utils.sub_pixel_fractions(xpix, ypix) - for i, j, k in zip(xs, ys, fracs): - flux_vector = specs[ref].value * weight * k / pixel_area - canvas_cube_hdu.data[:, j, i] += flux_vector - else: - x, y = int(xpix), int(ypix) - flux_vector = specs[ref].value * weight / pixel_area - canvas_cube_hdu.data[:, y, x] += flux_vector + for flux, x, y in self._make_cube_tablefields(specs): + canvas_cube_hdu.data[:, y, x] += flux - # 5. Add Background fields - for field in self.background_fields: - # bg_solid_angle = u.Unit(field.header["SOLIDANG"]).to(u.arcsec**-2) # float [arcsec-2] - # pixel_area = utils.from_currsys(self.meta["pixel_scale"]) ** 2 # float [arcsec2] - # area_factor = pixel_area * bg_solid_angle # float [arcsec2 * arcsec-2] - - # Cube should be in PHOTLAM arcsec-2 for SpectralTrace mapping - # spec = specs[field.header["SPEC_REF"]] * area_factor - spec = specs[field.header["SPEC_REF"]] - canvas_cube_hdu.data += spec[:, None, None].value + canvas_cube_hdu.data = sum(self._make_cube_backfields(specs), + start=canvas_cube_hdu.data) + # TODO: what's with this code?? # 6. Convert from PHOTLAM to ph/s/voxel # PHOTLAM = ph/s/cm-2/AA # area = m2, fov_waveset = um @@ -541,6 +577,7 @@ def make_cube_hdu(self): canvas_cube_hdu.data *= area.to(u.cm ** 2).value canvas_cube_hdu.data *= 1e4 # ph/s/AA/arcsec2 --> ph/s/um/arcsec2 + # TODO: what's with this code?? # bin_widths = np.diff(fov_waveset).to(u.AA).value # bin_widths = 0.5 * (np.r_[0, bin_widths] + np.r_[bin_widths, 0]) # canvas_cube_hdu.data *= bin_widths[:, None, None] @@ -552,14 +589,12 @@ def make_cube_hdu(self): "CUNIT3": "um", "CTYPE3": "WAVE"}) # ..todo: Add the log wavelength keyword here, if log scale is needed - - cube_hdu = canvas_cube_hdu # [ph s-1 AA-1 (arcsec-2)] - - return cube_hdu + return canvas_cube_hdu # [ph s-1 AA-1 (arcsec-2)] def volume(self, wcs_prefix=""): xs, ys = imp_utils.calc_footprint(self.header, wcs_suffix=wcs_prefix) - wave_corners = self.waverange + # FIXME: This is unused!! + # wave_corners = self.waverange self._volume = {"xs": [min(xs), max(xs)], "ys": [min(ys), max(ys)], "waves": self.waverange, @@ -569,28 +604,27 @@ def volume(self, wcs_prefix=""): @property def data(self): + """Return either hdu.data, image, cube, spectrum or None.""" if self.hdu is not None: - data = self.hdu.data - elif self.image is not None: - data = self.image - elif self.cube is not None: - data = self.cube - elif self.spectrum is not None: - data = self.spectrum - else: - data = None - - return data + return self.hdu.data + if self.image is not None: + return self.image + if self.cube is not None: + return self.cube + if self.spectrum is not None: + return self.spectrum + return None @property def corners(self): + """Return sky footprint, image plane footprint.""" sky_corners = imp_utils.calc_footprint(self.header) imp_corners = imp_utils.calc_footprint(self.header, "D") return sky_corners, imp_corners @property def waverange(self): - """Returns wavelength range in um [wave_min, wave_max]""" + """Return wavelength range in um [wave_min, wave_max].""" if self._waverange is None: wave_min = utils.quantify(self.meta["wave_min"], u.um).value wave_max = utils.quantify(self.meta["wave_max"], u.um).value @@ -599,55 +633,63 @@ def waverange(self): @property def wavelength(self): - """Returns central wavelength in um""" + """Return central wavelength in um.""" if self._wavelength is None: self._wavelength = np.average(self.waverange) return utils.quantify(self._wavelength, u.um) @property def waveset(self): - """Returns a wavelength vector in um""" - - field_cubes = self.cube_fields - if len(field_cubes) > 0: - i = np.argmax([cube.header["NAXIS3"] for cube in field_cubes]) - _waveset = fu.get_cube_waveset(field_cubes[i].header, + """Return a wavelength vector in um.""" + if field_cubes := self.cube_fields: + naxis3_max = np.argmax([cube.header["NAXIS3"] + for cube in field_cubes]) + _waveset = fu.get_cube_waveset(field_cubes[naxis3_max].header, return_quantity=True) - elif len(self.spectra) > 0: + elif self.spectra: wavesets = [spec.waveset for spec in self.spectra.values()] _waveset = np.concatenate(wavesets) else: _waveset = self.waverange * u.um - _waveset = np.unique(_waveset) - - return _waveset.to(u.um) + # TODO: tie the round to a global precision setting (this is defined + # better in another TODO somewhere...) + # The rounding is necessary to not end up with things like: + # 0.7 um + # 0.7000000000000001 um + # 0.7000000000000002 um + # and yes, that actually happend... + _waveset = np.unique(_waveset.to(u.um).round(10)) + return _waveset @property def cube_fields(self): + """Return list of non-BG_SRC ImageHDU fields with NAXIS=3.""" return [field for field in self.fields if isinstance(field, fits.ImageHDU) and field.header["NAXIS"] == 3 - and field.header.get("BG_SRC", False) is False] + and not field.header.get("BG_SRC", False)] @property def image_fields(self): + """Return list of non-BG_SRC ImageHDU fields with NAXIS=2.""" return [field for field in self.fields if isinstance(field, fits.ImageHDU) and field.header["NAXIS"] == 2 - and field.header.get("BG_SRC", False) is False] + and not field.header.get("BG_SRC", False)] @property def table_fields(self): + """Return list of Table fields.""" return [field for field in self.fields if isinstance(field, Table)] - @property def background_fields(self): + """Return list of BG_SRC ImageHDU fields.""" return [field for field in self.fields if isinstance(field, fits.ImageHDU) - and field.header.get("BG_SRC", False) is True] + and field.header.get("BG_SRC", False)] def __repr__(self): waverange = [self.meta["wave_min"].value, self.meta["wave_max"].value] diff --git a/scopesim/optics/fov_manager.py b/scopesim/optics/fov_manager.py index 819385e9..6c361af9 100644 --- a/scopesim/optics/fov_manager.py +++ b/scopesim/optics/fov_manager.py @@ -43,10 +43,11 @@ # """ from copy import deepcopy -import numpy as np from typing import TextIO from io import StringIO +from collections.abc import Iterable, MutableSequence +import numpy as np from astropy import units as u from . import image_plane_utils as ipu @@ -72,7 +73,7 @@ class FOVManager: All observation parameters as passed from UserCommands """ - def __init__(self, effects=[], **kwargs): + def __init__(self, effects=None, **kwargs): self.meta = {"area": "!TEL.area", "pixel_scale": "!INST.pixel_scale", "plate_scale": "!INST.plate_scale", @@ -94,9 +95,9 @@ def __init__(self, effects=[], **kwargs): params["meta"] = from_currsys({key: self.meta[key] for key in fvl_meta}) self.volumes_list = FovVolumeList(initial_volume=params) - self.effects = effects + self.effects = effects or [] self._fovs_list = [] - self.is_spectroscope = eu.is_spectroscope(effects) + self.is_spectroscope = eu.is_spectroscope(self.effects) if from_currsys(self.meta["preload_fovs"]) is True: self._fovs_list = self.generate_fovs_list() @@ -173,11 +174,11 @@ def fovs(self): return self._fovs_list @property - def fov_footprints(self, which="both"): + def fov_footprints(self): return None -class FovVolumeList(FOVSetupBase): +class FovVolumeList(FOVSetupBase, MutableSequence): """ List of FOV volumes for FOVManager @@ -211,9 +212,9 @@ def __init__(self, initial_volume=None): "yd_min": 0, "yd_max": 0} - def split(self, axis, value, aperture_id=None): + def split(self, axis, value, aperture_id=None) -> None: """ - Splits the all volumes that include axis=value into two. + Split the all volumes that include axis=value into two. - Loop through all volume dict - Find any entries where min < value < max @@ -221,8 +222,8 @@ def split(self, axis, value, aperture_id=None): Parameters ---------- - axis : str, list of str - "wave", "x", "y" + axis : {"wave", "x", "y"}, or list thereof + Which axis (``str``) or axes (``list[str]``) to use. value : float, list of floats aperture_id : int, optional Default None. If ``None``, split all volumes. If ``int``, only split @@ -244,20 +245,25 @@ def split(self, axis, value, aperture_id=None): if isinstance(axis, (tuple, list)): for ax, val in zip(axis, value): self.split(ax, val) - elif isinstance(value, (tuple, list, np.ndarray)): + return + + if isinstance(value, Iterable): for val in value: self.split(axis, val) - else: - for i, vol_old in enumerate(self.volumes): - if aperture_id in (vol_old["meta"]["aperture_id"], None) \ - and vol_old[f"{axis}_min"] < value \ - and vol_old[f"{axis}_max"] > value: - vol_new = deepcopy(vol_old) - vol_new[f"{axis}_min"] = value - vol_old[f"{axis}_max"] = value - self.volumes.insert(i+1, vol_new) - - def shrink(self, axis, values, aperture_id=None): + return + + for vol in self: + if (aperture_id is not None and + aperture_id != vol["meta"]["aperture_id"]): + continue + if vol[f"{axis}_min"] >= value or vol[f"{axis}_max"] <= value: + continue + new_vol = deepcopy(vol) + new_vol[f"{axis}_min"] = value + vol[f"{axis}_max"] = value + self.insert(self.index(vol) + 1, new_vol) + + def shrink(self, axis, values, aperture_id=None) -> None: """ - Loop through all volume dict - Replace any entries where min < values.min @@ -265,8 +271,8 @@ def shrink(self, axis, values, aperture_id=None): Parameters ---------- - axis : str - "wave", "x", "y" + axis : {"wave", "x", "y"} or list thereof + Which axis (``str``) or axes (``list[str]``) to use. values : list of 2 floats [min, max], [min, None], [None, max] aperture_id : int, optional @@ -283,47 +289,52 @@ def shrink(self, axis, values, aperture_id=None): """ + # FIXME: Isn't this method just the same as setting self.volumes to the + # output list of self.extract()?? Except the None values. if isinstance(axis, (tuple, list)): for ax, val in zip(axis, values): self.shrink(ax, val) - else: - to_pop = [] + return + + to_pop = [] + + for vol in self: + if (aperture_id is not None and + aperture_id != vol["meta"]["aperture_id"]): + continue if values[0] is not None: - for i, vol in enumerate(self.volumes): - if aperture_id in (vol["meta"]["aperture_id"], None): - if vol[f"{axis}_max"] <= values[0]: - to_pop.append(i) - elif vol[f"{axis}_min"] < values[0]: - vol[f"{axis}_min"] = values[0] + if vol[f"{axis}_max"] <= values[0]: + to_pop.append(self.index(vol)) + continue + vol[f"{axis}_min"] = max(values[0], vol[f"{axis}_min"]) if values[1] is not None: - for i, vol in enumerate(self.volumes): - if aperture_id in (vol["meta"]["aperture_id"], None): - if vol[f"{axis}_min"] >= values[1]: - to_pop.append(i) - if vol[f"{axis}_max"] > values[1]: - vol[f"{axis}_max"] = values[1] + if vol[f"{axis}_min"] >= values[1]: + to_pop.append(self.index(vol)) + continue + vol[f"{axis}_max"] = min(values[1], vol[f"{axis}_max"]) - for i in sorted(to_pop)[::-1]: - self.volumes.pop(i) + for idx in reversed(sorted(to_pop)): + self.pop(idx) def extract(self, axes, edges, aperture_id=None): """ - Returns new volumes from within all existing volumes + Return new volumes from within all existing volumes. This method DOES NOT alter the existing self.volumes list To include the returned volumes, add them to the self.volumes list Parameters ---------- - axes : str, list of str - "wave", "x", "y" + axes : list of either {"wave", "x", "y"} + Which axis (``list`` of single ``str``) or axes (``list[str]``) + to use. Must be ``list`` in either case. edges : list, tuple of lists Edge points for each axes listed aperture_id : int, optional Default None. If ``None``, extract from all volumes. If ``int``, - only extract from volumes with this ``aperture_id`` in the meta dict + only extract from volumes with this `aperture_id` in the meta dict Examples -------- @@ -345,38 +356,36 @@ def extract(self, axes, edges, aperture_id=None): A list of all new volumes extracted from existing volumes """ - new_vols = [] - for old_vol in self.volumes: - if aperture_id in (old_vol["meta"]["aperture_id"], None): - add_flag = True - new_vol = deepcopy(old_vol) + def _get_new_vols(): + for vol in self: + if (aperture_id is not None and + aperture_id != vol["meta"]["aperture_id"]): + continue + if not all(_volume_in_range(vol, axis, edge) for axis, edge + in zip(axes, edges)): + continue + new_vol = deepcopy(vol) for axis, edge in zip(axes, edges): - if edge[0] <= old_vol[f"{axis}_max"] and \ - edge[1] >= old_vol[f"{axis}_min"]: - new_vol[f"{axis}_min"] = max(edge[0], old_vol[f"{axis}_min"]) - new_vol[f"{axis}_max"] = min(edge[1], old_vol[f"{axis}_max"]) - else: - add_flag = False + new_vol[f"{axis}_min"] = max(edge[0], vol[f"{axis}_min"]) + new_vol[f"{axis}_max"] = min(edge[1], vol[f"{axis}_max"]) + yield new_vol - if add_flag is True: - new_vols.append(new_vol) - - return new_vols + return list(_get_new_vols()) def __len__(self): return len(self.volumes) - def __iter__(self): - return iter(self.volumes) + def __getitem__(self, index): + return self.volumes[index] - def __getitem__(self, key): - return self.volumes[key] + def __setitem__(self, index, value): + self.volumes[index] = value - def __setitem__(self, key, value): - self.volumes[item] = value + def __delitem__(self, index): + del self.volumes[index] - def __delitem__(self, key): - del self.volumes[key] + def insert(self, index, value): + self.volumes.insert(index, value) def write_string(self, stream: TextIO) -> None: """Write formatted string representation to I/O stream""" @@ -403,21 +412,18 @@ def __str__(self) -> str: output = str_stream.getvalue() return output - def __iadd__(self, other): - if isinstance(other, list): - self.volumes += other - else: - raise ValueError(f"Can only add lists of volumes: {other}") - - return self - def __add__(self, other): + # TODO: Is this functionality actually used anywhere? new_self = deepcopy(self) new_self += other return new_self +def _volume_in_range(vol: dict, axis: str, edge) -> bool: + return edge[0] <= vol[f"{axis}_max"] and edge[1] >= vol[f"{axis}_min"] + + # Spectroscopy FOV setup # ====================== # diff --git a/scopesim/optics/fov_manager_utils.py b/scopesim/optics/fov_manager_utils.py index 05165939..f97db5aa 100644 --- a/scopesim/optics/fov_manager_utils.py +++ b/scopesim/optics/fov_manager_utils.py @@ -1,4 +1,7 @@ +import logging from copy import deepcopy +from itertools import product +from more_itertools import pairwise import numpy as np from astropy import units as u @@ -9,6 +12,8 @@ from ..effects.effects_utils import get_all_effects from ..utils import check_keys +# TODO: Where are all these functions used?? + def get_3d_shifts(effects, **kwargs): """ @@ -41,9 +46,11 @@ def get_3d_shifts(effects, **kwargs): shifts = [eff.fov_grid(which="shifts", **kwargs) for eff in effects] old_bin_edges = [shift[0] for shift in shifts if len(shift[0]) >= 2] + # TODO: could this use combine_wavesets? new_bin_edges = np.unique(np.sort(np.concatenate(old_bin_edges), kind="stable")) + # TODO: could this be zeros_like? x_shifts = np.zeros(len(new_bin_edges)) y_shifts = np.zeros(len(new_bin_edges)) # .. todo:: replace the 1e-7 with a variable in !SIM @@ -99,25 +106,27 @@ def get_imaging_waveset(effects_list, **kwargs): wave_bin_edges = [filt.fov_grid(which="waveset", **kwargs) for filt in filters] - if len(wave_bin_edges) > 0: - kwargs["wave_min"] = np.max([w[0].value for w in wave_bin_edges]) - kwargs["wave_max"] = np.min([w[1].value for w in wave_bin_edges]) + if wave_bin_edges: + kwargs["wave_min"] = max(wave[0].value for wave in wave_bin_edges) + kwargs["wave_max"] = min(wave[1].value for wave in wave_bin_edges) + # Bit confusing... wave_bin_edges = [[kwargs["wave_min"], kwargs["wave_max"]]] if kwargs["wave_min"] > kwargs["wave_max"]: - raise ValueError(f"Filter wavelength ranges do not overlap: {wave_bin_edges}") + raise ValueError("Filter wavelength ranges do not overlap: " + f"{wave_bin_edges}.") # ..todo: add in Atmospheric dispersion and ADC here for effect_class in [efs.PSF]: - effects = get_all_effects(effects_list, effect_class) - for eff in effects: + for eff in get_all_effects(effects_list, effect_class): waveset = eff.fov_grid(which="waveset", **kwargs) if waveset is not None: - wave_bin_edges += [waveset] + wave_bin_edges.append(waveset) - wave_bin_edges = combine_wavesets(wave_bin_edges) + wave_bin_edges = combine_wavesets(*wave_bin_edges) - if len(wave_bin_edges) == 0: + if not wave_bin_edges: + # This is already set at the top, why again here? wave_bin_edges = [kwargs["wave_min"], kwargs["wave_max"]] return wave_bin_edges @@ -125,7 +134,7 @@ def get_imaging_waveset(effects_list, **kwargs): def get_imaging_headers(effects, **kwargs): """ - Returns a list of Header objects for each of the FieldOfVIew objects + Return a generator of Header objects for each of the FieldOfVIew objects. Parameters ---------- @@ -135,7 +144,7 @@ def get_imaging_headers(effects, **kwargs): Returns ------- - hdrs : list of Header objects + hdrs : generator of Header objects Notes ----- @@ -163,33 +172,39 @@ def get_imaging_headers(effects, **kwargs): aperture_effects = get_all_effects(effects, (efs.ApertureMask, efs.SlitWheel, efs.ApertureList)) - if len(aperture_effects) == 0: + if not aperture_effects: detector_arrays = get_all_effects(effects, efs.DetectorList) - if len(detector_arrays) > 0: - aperture_effects += [detarr.fov_grid(which="edges", - pixel_scale=pixel_scale) - for detarr in detector_arrays] - else: - raise ValueError("No ApertureMask or DetectorList was provided. At " - "least one must be passed to make an ImagePlane: " - f"{effects}") - + if not detector_arrays: + raise ValueError("No ApertureMask or DetectorList was provided. " + "At least one must be passed to make an " + f"ImagePlane: {effects}.") + aperture_effects.extend( + detarr.fov_grid(which="edges", pixel_scale=pixel_scale) + for detarr in detector_arrays) + + # FIXME: all of this is a bit inconsistent; fov_grid(which="edges" is + # called afterwards, but when looking in detector_arrays, the same + # is called immediately; does that even work? is this all tested?? # get aperture headers from fov_grid() # - for-loop catches mutliple headers from ApertureList.fov_grid() - hdrs = [] - for ap_eff in aperture_effects: - # ..todo:: add this functionality to ApertureList effect - hdr = ap_eff.fov_grid(which="edges", pixel_scale=pixel_scale) - hdrs += hdr if isinstance(hdr, (list, tuple)) else [hdr] + def _get_hdrs(ap_effs): + for ap_eff in ap_effs: + # ..todo:: add this functionality to ApertureList effect + hdr = ap_eff.fov_grid(which="edges", pixel_scale=pixel_scale) + if isinstance(hdr, (list, tuple)): + yield from hdr + else: + yield hdr + headers = _get_hdrs(aperture_effects) # check size of aperture in pixels - split if necessary - sky_hdrs = [] - for hdr in hdrs: - if hdr["NAXIS1"] * hdr["NAXIS2"] > kwargs["max_segment_size"]: - sky_hdrs += imp_utils.split_header(hdr, kwargs["chunk_size"]) - else: - sky_hdrs += [hdr] - + def _get_sky_hdrs(hdrs): + for hdr in hdrs: + if hdr["NAXIS1"] * hdr["NAXIS2"] > kwargs["max_segment_size"]: + yield from imp_utils.split_header(hdr, kwargs["chunk_size"]) + else: + yield hdr + sky_hdrs = _get_sky_hdrs(headers) # ..todo:: Deal with the case that two or more ApertureMasks overlap # map the on-sky apertures directly to the image plane using plate_scale @@ -208,13 +223,12 @@ def get_imaging_headers(effects, **kwargs): dethdr = imp_utils.header_from_list_of_xy(x_det, y_det, pixel_size, "D") skyhdr.update(dethdr) - - return sky_hdrs + yield skyhdr def get_imaging_fovs(headers, waveset, shifts, **kwargs): """ - Returns a list of ``FieldOfView`` objects + Return a generator of ``FieldOfView`` objects. Parameters ---------- @@ -224,54 +238,50 @@ def get_imaging_fovs(headers, waveset, shifts, **kwargs): waveset : list of floats [um] N+1 wavelengths for N spectral layers - shifts : list of tuples + shifts : list of tuples (or actually arrays?) [deg] x,y shifts w.r.t to the optical axis plane. N shifts for N spectral layers Returns ------- - fovs : list of FieldOfView objects + fovs : generator of ``FieldOfView`` objects """ - - shift_waves = shifts["wavelengths"] # in [um] + # Ensure array for later indexing + shift_waves = np.array(shifts["wavelengths"]) # in [um] shift_dx = shifts["x_shifts"] # in [deg] shift_dy = shifts["y_shifts"] # combine the wavelength bins from 1D spectral effects and 3D shift effects - if len(shifts["wavelengths"]) > 0: - mask = (shift_waves > np.min(waveset)) * (shift_waves < np.max(waveset)) - waveset = combine_wavesets([waveset, shift_waves[mask]]) - - counter = 0 - fovs = [] + if shift_waves.size: + mask = (shift_waves > min(waveset)) * (shift_waves < max(waveset)) + waveset = combine_wavesets(waveset, shift_waves[mask]) - print(f"Preparing {(len(waveset)-1)*len(headers)} FieldOfViews", flush=True) + # Actually evaluating the generators here is only necessary for the log msg + waveset = list(waveset) + headers = list(headers) + logging.info("Preparing %d FieldOfViews", (len(waveset) - 1) * len(headers)) - for ii in range(len(waveset) - 1): - for hdr in headers: - # add any pre-instrument shifts to the FOV sky coords - wave_mid = 0.5 * (waveset[ii] + waveset[ii+1]) - x_shift = np.interp(wave_mid, shift_waves, shift_dx) - y_shift = np.interp(wave_mid, shift_waves, shift_dy) - - fov_hdr = deepcopy(hdr) - fov_hdr["CRVAL1"] += x_shift # headers are in [deg] - fov_hdr["CRVAL2"] += y_shift + combos = product(pairwise(waveset), headers) + for fov_id, ((wave_min, wave_max), hdr) in enumerate(combos): + # add any pre-instrument shifts to the FOV sky coords + wave_mid = 0.5 * (wave_min + wave_max) + x_shift = np.interp(wave_mid, shift_waves, shift_dx) + y_shift = np.interp(wave_mid, shift_waves, shift_dy) - # define the wavelength range for the FOV - waverange = [waveset[ii], waveset[ii + 1]] + fov_hdr = deepcopy(hdr) + fov_hdr["CRVAL1"] += x_shift # headers are in [deg] + fov_hdr["CRVAL2"] += y_shift - # Make the FOV - fov = FieldOfView(fov_hdr, waverange, id=counter, **kwargs) - fovs += [fov] - counter += 1 + # define the wavelength range for the FOV + waverange = [wave_min, wave_max] - return fovs + # Make the FOV + yield FieldOfView(fov_hdr, waverange, id=fov_id, **kwargs) def get_spectroscopy_headers(effects, **kwargs): - + """Return generator of Header objects.""" required_keys = ["pixel_scale", "plate_scale", "wave_min", "wave_max"] check_keys(kwargs, required_keys, action="error") @@ -285,15 +295,15 @@ def get_spectroscopy_headers(effects, **kwargs): efs.SlitWheel, efs.ApertureMask)) - if len(surface_list_effects) > 0: + if surface_list_effects: waves = surface_list_effects[0].fov_grid(which="waveset") if len(waves) == 2: kwargs["wave_min"] = np.max([waves[0].value, kwargs["wave_min"]]) kwargs["wave_max"] = np.min([waves[1].value, kwargs["wave_max"]]) - if len(detector_list_effects) > 0: + if detector_list_effects: implane_hdr = detector_list_effects[0].image_plane_header - elif len(spec_trace_effects) > 0: + elif spec_trace_effects: implane_hdr = spec_trace_effects[0].image_plane_header else: raise ValueError("Missing a way to determine the image plane size") @@ -304,6 +314,7 @@ def get_spectroscopy_headers(effects, **kwargs): f"{spec_trace_effects}") spec_trace = spec_trace_effects[0] + # TODO: The following is WET with the code in get_imaging_headers sky_hdrs = [] for ap_eff in aperture_effects: # if ApertureList, a list of ApertureMask headers is returned @@ -320,19 +331,23 @@ def get_spectroscopy_headers(effects, **kwargs): plate_scale=kwargs["plate_scale"] ) for sky_hdr in sky_hdrs] - fov_headers = [hdr for hdr_list in fov_headers for hdr in hdr_list] - # ..todo: check that each header is not larger than chunk_size - - return fov_headers + for hdr_list in fov_headers: + for hdr in hdr_list: + yield hdr + # TODO: check that each header is not larger than chunk_size + # that's already done in get_imaging_headers, isn't it? -def get_spectroscopy_fovs(headers, shifts, effects=[], **kwargs): +def get_spectroscopy_fovs(headers, shifts, effects=None, **kwargs): + """Return a generator of ``FieldOfView`` objects.""" + if effects is None: + effects = [] shift_waves = shifts["wavelengths"] # in [um] shift_dx = shifts["x_shifts"] # in [deg] shift_dy = shifts["y_shifts"] - print(f"Preparing {len(headers)} FieldOfViews", flush=True) + logging.info("Preparing %d FieldOfViews", len(headers)) apertures = get_all_effects(effects, (efs.ApertureList, efs.ApertureMask)) masks = [ap.fov_grid(which="masks") for ap in apertures] @@ -343,8 +358,7 @@ def get_spectroscopy_fovs(headers, shifts, effects=[], **kwargs): elif isinstance(mask, np.ndarray): mask_dict[len(mask_dict)] = mask - fovs = [] - for ii, hdr in enumerate(headers): + for fov_id, hdr in enumerate(headers): # add any pre-instrument shifts to the FOV sky coords wave_mid = hdr["WAVE_MID"] x_shift = np.interp(wave_mid, shift_waves, shift_dx) @@ -355,29 +369,30 @@ def get_spectroscopy_fovs(headers, shifts, effects=[], **kwargs): fov_hdr["CRVAL2"] += y_shift # Make the FOV - fov = FieldOfView(fov_hdr, waverange=[hdr["WAVE_MIN"], hdr["WAVE_MAX"]], - **kwargs) + waverange = [hdr["WAVE_MIN"], hdr["WAVE_MAX"]] + fov = FieldOfView(fov_hdr, waverange=waverange, **kwargs) fov.meta["distortion"]["rotation"] = hdr["ROTANGD"] fov.meta["distortion"]["shear"] = hdr["SKEWANGD"] fov.meta["conserve_image"] = hdr["IMG_CONS"] - fov.meta["fov_id"] = ii + # TODO: In the other function, the id is set via the contructor. + # What's the difference? + fov.meta["fov_id"] = fov_id fov.meta["aperture_id"] = hdr["APERTURE"] # .. todo: get these masks working # there needs to be fov_grid(which="mask") in ApertureList/Mask # fov.mask = mask_dict[hdr["APERTURE"]] - fovs += [fov] + yield fov - return fovs - -def combine_wavesets(waveset_list): +# FIXME: This functions doesn't seem to be covered by any separate unit test. +def combine_wavesets(*wavesets): """ - Joins and sorts several sets of wavelengths into a single 1D array + Join and sorts several sets of wavelengths into a single 1D array. Parameters ---------- - waveset_list : list + wavesets : one or more iterables A group of wavelength arrays or lists Returns @@ -385,12 +400,22 @@ def combine_wavesets(waveset_list): wave_set : np.ndarray Combined set of wavelengths + Note + ---- + This assumes that all wavesets are given in the same unit! """ - wave_set = [] - for wbe in waveset_list: # wbe = waveset bin edges - wbe = wbe.value if isinstance(wbe, u.Quantity) else wbe - wave_set += list(wbe) - # ..todo:: set variable in !SIM.computing for rounding to the 7th decimal - wave_set = np.unique(np.round(np.sort(wave_set, kind="stable"), 7)) - + # TODO: set variable in !SIM.computing for rounding to the 7th decimal + decimals = 7 + + def _get_waves(waves): + for wave in waves: + if isinstance(wave, u.Quantity): + round_wave = wave.round(decimals).value + else: + round_wave = np.round(wave, decimals) + yield from round_wave + + # NOTE: This function previously used np.sort(wave_set, kind="stable"). + # If any issues occur with the buitin sorted, go back to that! + wave_set = sorted(set(_get_waves(wavesets))) return wave_set diff --git a/scopesim/optics/fov_utils.py b/scopesim/optics/fov_utils.py index 8c554b49..86747eb4 100644 --- a/scopesim/optics/fov_utils.py +++ b/scopesim/optics/fov_utils.py @@ -329,7 +329,11 @@ def extract_area_from_imagehdu(imagehdu, fov_volume): # OC [2021-12-14] if fov range is not covered by the source return nothing if not np.any(mask): - print(f"FOV {fov_waves[0]} um - {fov_waves[1]} um: not covered by Source") + logging.warning("FOV %s um - %s um: not covered by Source", + fov_waves[0], fov_waves[1]) + # FIXME: returning None here breaks the principle that a function + # should always return the same type. Maybe this should + # instead raise an exception that's caught higher up... return None i0p, i1p = np.where(mask)[0][0], np.where(mask)[0][-1] @@ -360,6 +364,9 @@ def extract_area_from_imagehdu(imagehdu, fov_volume): data = imagehdu.data[y0p:y1p, x0p:x1p] new_hdr["SPEC_REF"] = hdr.get("SPEC_REF") + if not data.size: + logging.warning("Empty image HDU.") + new_imagehdu = fits.ImageHDU(data=data) new_imagehdu.header.update(new_hdr) diff --git a/scopesim/optics/image_plane_utils.py b/scopesim/optics/image_plane_utils.py index 116c6c18..6944175f 100644 --- a/scopesim/optics/image_plane_utils.py +++ b/scopesim/optics/image_plane_utils.py @@ -756,14 +756,15 @@ def val2pix(header, a, b, wcs_suffix=""): if isinstance(header, fits.ImageHDU): header = header.header - if "PC1_1"+s in header: - pc11 = header["PC1_1"+s] - pc12 = header["PC1_2"+s] - pc21 = header["PC2_1"+s] - pc22 = header["PC2_2"+s] + pckeys = [key + s for key in ["PC1_1", "PC1_2", "PC2_1", "PC2_2"]] + if all(key in header for key in pckeys): + pc11, pc12, pc21, pc22 = (header[key] for key in pckeys) else: pc11, pc12, pc21, pc22 = 1, 0, 0, 1 + if (pc11 * pc22 + pc12 * pc21) != 1.0: + logging.error("PC matrix det != 1.0") + da = float(header["CDELT1"+s]) db = float(header["CDELT2"+s]) x0 = float(header["CRPIX1"+s]) diff --git a/scopesim/tests/tests_optics/test_FieldOfView.py b/scopesim/tests/tests_optics/test_FieldOfView.py index 99ff7152..45e571f4 100644 --- a/scopesim/tests/tests_optics/test_FieldOfView.py +++ b/scopesim/tests/tests_optics/test_FieldOfView.py @@ -140,8 +140,9 @@ def test_contains_all_fields_inside_fov(self, basic_fov_header, assert isinstance(the_fov.fields[2], Table) -@pytest.mark.xfail(reason="apply make_cube's fov.waveset available to the outside ") +# @pytest.mark.xfail(reason="apply make_cube's fov.waveset available to the outside ") class TestMakeCube: + @pytest.mark.xfail(reason="apply make_cube's fov.waveset available to the outside ") def test_makes_cube_from_table(self): src_table = so._table_source() # 10x10" @ 0.2"/pix, [0.5, 2.5]m @ 0.02µm fov = _fov_190_210_um() @@ -163,6 +164,7 @@ def test_makes_cube_from_table(self): plt.imshow(cube.data[0, :, :], origin="lower") plt.show() + @pytest.mark.xfail(reason="apply make_cube's fov.waveset available to the outside ") def test_makes_cube_from_imagehdu(self): src_image = so._image_source() # 10x10" @ 0.2"/pix, [0.5, 2.5]m @ 0.02µm fov = _fov_190_210_um() @@ -181,6 +183,7 @@ def test_makes_cube_from_imagehdu(self): plt.imshow(cube.data[0, :, :], origin="lower") plt.show() + @pytest.mark.xfail(reason="apply make_cube's fov.waveset available to the outside ") def test_makes_cube_from_other_cube_imagehdu(self): import scopesim as sim sim.rc.__currsys__["!SIM.spectral.spectral_bin_width"] = 0.01 @@ -201,6 +204,7 @@ def test_makes_cube_from_other_cube_imagehdu(self): plt.imshow(cube.data[0, :, :], origin="lower") plt.show() + @pytest.mark.xfail(reason="apply make_cube's fov.waveset available to the outside ") def test_makes_cube_from_two_similar_cube_imagehdus(self): src_cube = so._cube_source() + so._cube_source(dx=1) # 2 cubes 10x10" @ 0.2"/pix, [0.5, 2.5]m @ 0.02µm fov = _fov_197_202_um() @@ -220,6 +224,7 @@ def test_makes_cube_from_two_similar_cube_imagehdus(self): plt.imshow(cube.data[0, :, :], origin="lower") plt.show() + @pytest.mark.xfail(reason="apply make_cube's fov.waveset available to the outside ") def test_makes_cube_from_all_types_of_source_object(self): src_all = so._table_source() + \ so._image_source(dx=-4, dy=-4) + \ @@ -271,7 +276,7 @@ def test_cube_has_full_wcs(self, src): assert "CTYPE3" in cube.header -@pytest.mark.xfail(reason="revisit fov.waveset e.g. use make_cube waveset") +# @pytest.mark.xfail(reason="revisit fov.waveset e.g. use make_cube waveset") class TestMakeImage: def test_makes_image_from_table(self): src_table = so._table_source() # 10x10" @ 0.2"/pix, [0.5, 2.5]m @ 0.02µm @@ -317,6 +322,7 @@ def test_makes_image_from_image(self): plt.imshow(img.data, origin="lower") plt.show() + @pytest.mark.xfail(reason="revisit fov.waveset e.g. use make_cube waveset") def test_makes_image_from_cube(self): src_cube = so._cube_source() # 10x10" @ 0.2"/pix, [0.5, 2.5]m @ 0.02µm fov = _fov_197_202_um() @@ -334,6 +340,7 @@ def test_makes_image_from_cube(self): assert out_sum == approx(in_sum) + @pytest.mark.xfail(reason="revisit fov.waveset e.g. use make_cube waveset") def test_makes_image_from_all_types_of_source_object(self): src_all = so._table_source() + \ so._image_source(dx=-4, dy=-4) + \ @@ -370,7 +377,7 @@ def test_makes_image_from_all_types_of_source_object(self): plt.show() -@pytest.mark.xfail(reason="revisit fov.waveset e.g. use make_cube waveset") +# @pytest.mark.xfail(reason="revisit fov.waveset e.g. use make_cube waveset") class TestMakeSpectrum: def test_make_spectrum_from_table(self): src_table = so._table_source() # 10x10" @ 0.2"/pix, [0.5, 2.5]m @ 0.02µm diff --git a/scopesim/tests/tests_optics/test_fov_manager_utils.py b/scopesim/tests/tests_optics/test_fov_manager_utils.py index 1007af7e..63468858 100644 --- a/scopesim/tests/tests_optics/test_fov_manager_utils.py +++ b/scopesim/tests/tests_optics/test_fov_manager_utils.py @@ -13,11 +13,31 @@ from scopesim.tests.mocks.py_objects import trace_list_objects as tlo +@pytest.fixture(scope="function") +def wave_kwargs(): + return {"wave_min": 0.5, "wave_max": 2.5} + + +@pytest.fixture(scope="function") +def th_filt(): + return eo._filter_tophat_curve() + + @pytest.fixture(scope="function") def full_trace_list(): return tlo.make_trace_hdulist() +@pytest.fixture(scope="function") +def spec_hdrs(full_trace_list): + params = {"pixel_scale": 0.1, "plate_scale": 0.1, + "wave_min": 0.7, "wave_max": 2.5} + spt = SpectralTraceList(hdulist=full_trace_list, **params) + apm = apo._basic_aperture() + hdrs = fm_utils.get_spectroscopy_headers(effects=[spt, apm], **params) + return list(hdrs) + + PLOTS = False @@ -73,35 +93,35 @@ def test_combined_shifts_reduced_to_usable_number(self, sub_pix_frac): class TestGetImagingWaveset: - def test_returns_default_wave_range_when_passed_no_effects(self): - kwargs = {"wave_min": 0.5, "wave_max": 2.5} - wave_bin_edges = fm_utils.get_imaging_waveset([], **kwargs) + @pytest.mark.usefixtures("wave_kwargs") + def test_returns_default_wave_range_when_passed_no_effects(self, wave_kwargs): + wave_bin_edges = fm_utils.get_imaging_waveset([], **wave_kwargs) assert len(wave_bin_edges) == 2 - def test_returns_waveset_of_filter(self): - filt = eo._filter_tophat_curve() - kwargs = {"wave_min": 0.5, "wave_max": 2.5} - wave_bin_edges = fm_utils.get_imaging_waveset([filt], **kwargs) + @pytest.mark.usefixtures("wave_kwargs", "th_filt") + def test_returns_waveset_of_filter(self, wave_kwargs, th_filt): + wave_bin_edges = fm_utils.get_imaging_waveset([th_filt], **wave_kwargs) assert len(wave_bin_edges) == 2 - def test_returns_waveset_of_psf(self): + @pytest.mark.usefixtures("wave_kwargs") + def test_returns_waveset_of_psf(self, wave_kwargs): psf = eo._const_psf() - kwargs = {"wave_min": 0.5, "wave_max": 2.5} - wave_bin_edges = fm_utils.get_imaging_waveset([psf], **kwargs) + wave_bin_edges = fm_utils.get_imaging_waveset([psf], **wave_kwargs) assert len(wave_bin_edges) == 4 - def test_returns_waveset_of_psf_and_filter(self): - filt = eo._filter_tophat_curve() + @pytest.mark.usefixtures("wave_kwargs", "th_filt") + def test_returns_waveset_of_psf_and_filter(self, wave_kwargs, th_filt): psf = eo._const_psf() - kwargs = {"wave_min": 0.5, "wave_max": 2.5} - wave_bin_edges = fm_utils.get_imaging_waveset([filt, psf], **kwargs) + wave_bin_edges = fm_utils.get_imaging_waveset([th_filt, psf], + **wave_kwargs) assert len(wave_bin_edges) == 4 - def test_returns_waveset_of_ncpa_psf_inside_filter_edges(self): - filt = eo._filter_tophat_curve() + @pytest.mark.usefixtures("wave_kwargs", "th_filt") + def test_returns_waveset_of_ncpa_psf_inside_filter_edges(self, wave_kwargs, + th_filt): psf = eo._ncpa_psf() - kwargs = {"wave_min": 0.5, "wave_max": 2.5} - wave_bin_edges = fm_utils.get_imaging_waveset([psf, filt], **kwargs) + wave_bin_edges = fm_utils.get_imaging_waveset([psf, th_filt], + **wave_kwargs) assert min(wave_bin_edges) == 1. assert max(wave_bin_edges) == 2. assert len(wave_bin_edges) == 9 @@ -111,7 +131,7 @@ class TestGetImagingHeaders: def test_throws_error_if_not_all_keywords_are_passed(self): apm = eo._img_aperture_mask() with pytest.raises(ValueError): - fm_utils.get_imaging_headers([apm]) + list(fm_utils.get_imaging_headers([apm])) def test_returns_set_of_headers_from_aperture_effects(self): apm = eo._img_aperture_mask(array_dict={"x": [-1.28, 1., 1., -1.28], @@ -129,6 +149,8 @@ def test_returns_set_of_headers_from_detector_list_effect(self): kwargs = {"pixel_scale": 0.004, "plate_scale": 0.26666666666, "max_segment_size": 2048 ** 2, "chunk_size": 1024} hdrs = fm_utils.get_imaging_headers([det], **kwargs) + # Evaluate generator for testing + hdrs = list(hdrs) area_sum = np.sum([hdr["NAXIS1"] * hdr["NAXIS2"] for hdr in hdrs]) assert area_sum == 4096**2 @@ -161,6 +183,8 @@ def test_returns_fov_objects_for_basic_input(self): "max_segment_size": 100 ** 2, "chunk_size": 100} hdrs = fm_utils.get_imaging_headers([apm], **kwargs) + # Evaluate generator for testing + hdrs = list(hdrs) waveset = np.linspace(1, 2, 6) shifts = {"wavelengths": np.array([1, 2]), "x_shifts": np.zeros(2), @@ -168,7 +192,10 @@ def test_returns_fov_objects_for_basic_input(self): fovs = fm_utils.get_imaging_fovs(headers=hdrs, waveset=waveset, shifts=shifts) + # Evaluate generator for testing + fovs = list(fovs) assert len(fovs) == (len(waveset)-1) * len(hdrs) + assert fovs if PLOTS: from scopesim.optics.image_plane_utils import calc_footprint @@ -189,19 +216,13 @@ def test_returns_fov_objects_for_basic_input(self): plt.show() -@pytest.mark.usefixtures("full_trace_list") +@pytest.mark.usefixtures("spec_hdrs") class TestGetSpectroscopyHeaders: - def test_returns_headers(self, full_trace_list): - params = {"pixel_scale": 0.1, "plate_scale": 0.1, - "wave_min": 0.7, "wave_max": 2.5} - spt = SpectralTraceList(hdulist=full_trace_list, **params) - apm = apo._basic_aperture() - - hdrs = fm_utils.get_spectroscopy_headers(effects=[spt, apm], **params) - assert all([isinstance(hdr, PoorMansHeader) for hdr in hdrs]) + def test_returns_headers(self, spec_hdrs): + assert all([isinstance(hdr, PoorMansHeader) for hdr in spec_hdrs]) if PLOTS: - for hdr in hdrs: + for hdr in spec_hdrs: x = np.array([0, hdr["NAXIS1"], hdr["NAXIS1"], 0]) y = np.array([0, 0, hdr["NAXIS2"], hdr["NAXIS2"]]) xw, yw = pix2val(hdr, x, y, "D") @@ -210,19 +231,15 @@ def test_returns_headers(self, full_trace_list): plt.show() -@pytest.mark.usefixtures("full_trace_list") +@pytest.mark.usefixtures("spec_hdrs") class TestGetSpectroscopyFOVs: - def test_returns_fovs(self, full_trace_list): - params = {"pixel_scale": 0.1, "plate_scale": 0.1, - "wave_min": 0.7, "wave_max": 2.5} - spt = SpectralTraceList(hdulist=full_trace_list, **params) - apm = apo._basic_aperture() - + def test_returns_fovs(self, spec_hdrs): shifts = {"wavelengths": np.array([0.7, 2.5]), "x_shifts": np.array([0, 0]), "y_shifts": np.array([0, 1/3600.])} - hdrs = fm_utils.get_spectroscopy_headers(effects=[spt, apm], **params) - fovs = fm_utils.get_spectroscopy_fovs(hdrs, shifts) + fovs = fm_utils.get_spectroscopy_fovs(spec_hdrs, shifts) + # Evaluate generator for testing + fovs = list(fovs) assert all([isinstance(fov, FieldOfView) for fov in fovs])