From 029692ba8830ecc7beb40c7609e4df0a41ca20be Mon Sep 17 00:00:00 2001 From: Anthony Horton Date: Tue, 26 Nov 2019 10:57:57 +1100 Subject: [PATCH] Image simulation (#16) * Added dark frame generation to Camera, moved other astroimsim code * Initial implementations of image simulation methods * Docs build fix * Remove now redundant astroimsim.py * Memory use reduction (no longer creates huge temporary arrays) * Fix off-by-one, array vs Quantity & floating point stopping bugs * Fix axis 18 (#20) * Axis order changed * Camera.py docstring improved * clipping below zero pixel values * Correcting test_psf.py * camera.py docstring improved * below zero assert added * below zero and infinit asserts added * combining pairs * `@pytest.mark.parameterize` implemented * Some problems with ... * something that works for now * Crval not set as per Issue #17 (#25) * use centre to always set CRVAL first thing See issue: https://github.com/AstroHuntsman/gunagala/issues/17 * created get_pixel_coords method for Imager * remove centre for get_pixel_coords & check for CRVAL * add :skip: InvalidTransformError to docs/index.rst * added test of get_pixel_coords * fix unit=deg problems in SkyCoord * also unit=deg fix * another unit=deg fix * split test and make scope function * make_noiseless_image unit=deg * kwargs docstrings & split star_ & centre_ kwargs * @bazkiaei trying to parametrize as it has been advised * @bazkiaei trying to parametrize as it has been advised * Revert "Merge remote-tracking branch 'upstream/image_sim' into improve-pytest-#26" This reverts commit 308168f633c2e4fade9a6e0848056b7ca4c0f674. * Revert "Merge remote-tracking branch 'upstream/image_sim' into improve-pytest-#26" This reverts commit 4b4bc04d7f671fd57da316618f17535ccceda22c. * Revert "Revert "Merge remote-tracking branch 'upstream/image_sim' into improve-pytest-#26"" This reverts commit 25a4c047ff85a2c2f166b04b6c11d870ded1bf91. * Revert "@bazkiaei" This reverts commit 94273a26e3ce4d81d5c4c4d616135727793ee9c8. * Revert "@bazkiaei" This reverts commit c811d411ddbad9feb44e395b3a29d972fab466d8. * Revert "Merge remote-tracking branch 'upstream/image_sim' into improve-pytest-#26" This reverts commit 308168f633c2e4fade9a6e0848056b7ca4c0f674. * Revert "Revert "Merge remote-tracking branch 'upstream/image_sim' into improve-pytest-#26"" This reverts commit 9eb810c2854c296a95232da4d940c7fd535d30ce. * Improving psf_pytest.py #26 (#28) The `test_psf.py` parametrised and decorated. * auto pep8 * fix pixellated call * add notimplemented error for analytical psf * fix real image generator bug * rename ZWO QE --- .travis.yml | 4 +- docs/index.rst | 4 +- gunagala/astroimsim.py | 213 ---------------- gunagala/camera.py | 70 +++++- gunagala/data/performance.yaml | 2 +- .../{ZWO_QE.csv => ZWO_QE_asi1600.csv} | 0 gunagala/imager.py | 232 +++++++++++++++--- gunagala/psf.py | 118 ++++++--- gunagala/tests/test_camera.py | 93 +++++++ gunagala/tests/test_imager.py | 23 ++ gunagala/tests/test_psf.py | 130 +++++++--- 11 files changed, 571 insertions(+), 318 deletions(-) delete mode 100644 gunagala/astroimsim.py rename gunagala/data/performance_data/{ZWO_QE.csv => ZWO_QE_asi1600.csv} (100%) diff --git a/.travis.yml b/.travis.yml index 6e7f4a3..253138b 100644 --- a/.travis.yml +++ b/.travis.yml @@ -88,8 +88,8 @@ matrix: EVENT_TYPE='pull_request push cron' #- os: linux # env: PYTHON_VERSION=2.7 ASTROPY_VERSION=lts - - os: linux - env: ASTROPY_VERSION=lts + # - os: linux + # env: ASTROPY_VERSION=lts # Try all python versions and Numpy versions. Since we can assume that # the Numpy developers have taken care of testing Numpy with different diff --git a/docs/index.rst b/docs/index.rst index afb7736..ce76b44 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -123,11 +123,13 @@ Reference/API :skip: PSF :skip: Sky :skip: Simple - :skip: Table :skip: MoffatPSF :skip: ZodiacalLight :skip: WCS + :skip: SkyCoord + :skip: CCDData :skip: interp1d + :skip: InvalidTransformError :no-inheritance-diagram: .. automodapi:: gunagala.optic diff --git a/gunagala/astroimsim.py b/gunagala/astroimsim.py deleted file mode 100644 index 315b523..0000000 --- a/gunagala/astroimsim.py +++ /dev/null @@ -1,213 +0,0 @@ -import numpy as np -from scipy.interpolate import RectSphereBivariateSpline, SmoothBivariateSpline, InterpolatedUnivariateSpline -from scipy.stats import poisson, norm, lognorm -from scipy.special import eval_chebyt -import astropy.io.fits as fits -import astropy.units as u -import astropy.constants as c -from astropy.coordinates import SkyCoord, GeocentricTrueEcliptic, get_sun, Angle -from astropy.time import Time -from astropy.wcs import WCS -from astropy.table import Table -import ccdproc - - -class Imager: - """ - Class representing an imaging instrument. - """ - def __init__(self, npix_x, npix_y, pixel_scale, aperture_area, throughput, filters, QE, gain, read_noise, temperature, zl): - - self.pixel_scale = pixel_scale - - # Construct a simple template WCS to store the focal plane configuration parameters - self.wcs = WCS(naxis=2) - self.wcs._naxis1 = npix_x - self.wcs._naxis2 = npix_y - self.wcs.wcs.crpix = [(npix_x + 1)/2, (npix_y + 1)/2] - self.wcs.wcs.cdelt = [pixel_scale.to(u.degree).value, pixel_scale.to(u.degree).value] - self.wcs.wcs.ctype = ['RA---TAN', 'DEC--TAN'] - - # Store throughput related parameters - self.aperture_area = aperture_area - self.throughput = throughput - self.filters = filters - self.QE = QE - - self.gain = gain - self.read_noise = read_noise - - self.zl = zl - - # Pre-calculate effective aperture areas. pivot wavelengths and sensitivity integrals - self._eff_areas = self._effective_areas() - self._pivot_waves = self._pivot_wavelengths() - self._sensitivities = self._sensitivity_integral() - - # Pre-calculate normalisation for observed ZL - self._zl_ep = self._zl_obs_ep() - - # Precalculate dark frame - self.dark_current, self.dark_frame = self._make_dark_frame(temperature) - - def _make_dark_frame(self, temperature, alpha=0.0488/u.Kelvin, beta=-12.772, shape=0.4, seed=None): - """ - Function to create a dark current 'image' in electrons per second per pixel given - an image sensor temperature and a set of coefficients for a simple dark current model. - Modal dark current for the image sensor as a whole is modelled as D.C. = 10**(alpha * T + beta) where - T is the temperature in Kelvin. Individual pixel dark currents are random uncorrelated values from - a log normal distribution so that there is a semi-realistic tail of 'hot pixels'. - For reproducible dark frames the random number generator seed can optionally be specified. - """ - temperature = temperature.to(u.Kelvin, equivalencies=u.equivalencies.temperature()) - mode = 10**(alpha * temperature + beta) * u.electron / (u.second) - scale = mode * np.exp(shape**2) - if seed: - np.random.seed(seed) - dark_frame = lognorm.rvs(shape, scale=scale, size=(self.wcs._naxis2, self.wcs._naxis1)) - - return mode, dark_frame - - def _effective_areas(self): - """ - Utility function to calculate the effective aperture area of for each filter as a function of - wavelength, i.e. aperture area * optical throughput * image sensor QE - """ - eff_areas = {} - - for (f_name, f_data) in self.filters.items(): - # Interpolate throughput data at same wavelengths as filter transmission data - t = np.interp(f_data['Wavelength'], self.throughput['Wavelength'], \ - self.throughput['Throughput']) * self.throughput['Throughput'].unit - # Interpolate QE data at same wavelengths as filter transmission data - q = np.interp(f_data['Wavelength'], self.QE['Wavelength'], self.QE['QE']) * self.QE['QE'].unit - eff_areas[f_name] = Table(names = ('Wavelength', 'Effective Area'), \ - data = (f_data['Wavelength'], self.aperture_area * t * f_data['Transmission'] * q)) - - return eff_areas - - def _pivot_wavelengths(self): - """ - Utility function to calculate the pivot wavelengths for each of the filters using the - effective area data (must be pre-calculated before calling this function). - """ - pivot_waves = {} - - # Generally this is definied in terms of system efficiency instead of aperture - # effective aperture area but they're equivalent as the aperture area factor - # cancels. - for (f_name, eff_data) in self._eff_areas.items(): - p1 = np.trapz(eff_data['Effective Area'] * eff_data['Wavelength'], x=eff_data['Wavelength']) - p2 = np.trapz(eff_data['Effective Area'] / eff_data['Wavelength'], x=eff_data['Wavelength']) - pivot_waves[f_name] = (p1/p2)**0.5 * eff_data['Wavelength'].unit - - return pivot_waves - - def _sensitivity_integral(self): - """ - Utility function to calculate the sensitivity integral for each of the filters, - i.e. the factor to convert a constant spectral flux density in F_lambda units - to a count rate in electrons per second. - """ - # Need to make sure units get preserved here. - sensitivities = {} - - for (f_name, eff_data) in self._eff_areas.items(): - s = np.trapz(eff_data['Wavelength'] * eff_data['Effective Area'], x=eff_data['Wavelength']) - s = s * eff_data['Effective Area'].unit * eff_data['Wavelength'].unit**2 * u.photon / (c.h * c.c) - sensitivities[f_name] = s.to((u.electron/u.second) / (u.Watt / (u.m**2 * u.micron))) - - return sensitivities - - def _zl_obs_ep(self): - """ - Utility function to pre-calculate observed ecliptic pole zodiacal light count rates - """ - # Integrate product of zodiacal light photon SFD and effective aperture area - # over wavelength to get observed ecliptic pole surface brightness for each filter. - # Note, these are constant with time so can and should precalculate once. - zl_ep = {} - - for f in self.filters.keys(): - electrons = np.zeros((self.wcs._naxis2, self.wcs._naxis1)) * u.electron / u.second - eff_area_interp = np.interp(self.zl.waves, self._eff_areas[f]['Wavelength'], \ - self._eff_areas[f]['Effective Area']) * \ - self._eff_areas[f]['Effective Area'].unit - zl_ep[f] = (np.trapz(self.zl.photon_sfd * eff_area_interp, x=self.zl.waves)) - - return zl_ep - - def get_pixel_coords(self, centre): - """ - Utility function to return a SkyCoord array containing the on sky position - of the centre of all the pixels in the image centre, given a SkyCoord for - the field centre - """ - # Ensure centre is a SkyCoord (this allows entering centre as a string) - if not isinstance(centre, SkyCoord): - centre = SkyCoord(centre) - - # Set field centre coordinates in internal WCS - self.wcs.wcs.crval = [centre.icrs.ra.value, centre.icrs.dec.value] - - # Arrays of pixel coordinates - XY = np.meshgrid(np.arange(self.wcs._naxis1), np.arange(self.wcs._naxis2)) - - # Convert to arrays of RA, dec (ICRS, decimal degrees) - RAdec = self.wcs.all_pix2world(XY[0], XY[1], 0) - - return SkyCoord(RAdec[0], RAdec[1], unit='deg') - - def make_noiseless_image(self, centre, time, f): - """ - Function to create a noiseless simulated image for a given image centre and observation time. - """ - electrons = np.zeros((self.wcs._naxis2, self.wcs._naxis1)) * u.electron / u.second - - # Calculate observed zodiacal light background. - # Get relative zodical light brightness for each pixel - # Note, side effect of this is setting centre of self.wcs - pixel_coords = self.get_pixel_coords(centre) - zl_rel = zl.relative_brightness(pixel_coords, time) - - # TODO: calculate area of each pixel, for now use nominal pixel scale^2 - # Finally multiply to get an observed zodical light image - zl_obs = self.zl_obs_ep * zl_rel * self.pixel_scale**2 - - electrons += zl_obs - - noiseless = ccdproc.CCDData(electrons, wcs=self.wcs) - - return noiseless - - def make_image_real(self, noiseless, exp_time, subtract_dark = False): - """ - Given a noiseless simulated image in electrons per pixel add dark current, - Poisson noise and read noise, and convert to ADU using the predefined gain. - """ - # Scale photoelectron rates by exposure time - data = noiseless.data * noiseless.unit * exp_time - # Add dark current - data += self.dark_frame * exp_time - # Force to electron units - data = data.to(u.electron) - # Apply Poisson noise. This is not unit-aware, need to restore them manually - data = (poisson.rvs(data/u.electron)).astype('float64') * u.electron - # Apply read noise. Again need to restore units manually - data += norm.rvs(scale=self.read_noise/u.electron, size=data.shape) * u.electron - # Optionally subtract a Perfect Dark - if subtract_dark: - data -= (self.dark_frame * exp_time).to(u.electron) - # Convert to ADU - data /= self.gain - # Force to adu (just here to catch unit problems) - data = data.to(u.adu) - # 'Analogue to digital conversion' - data = np.where(data < 2**16 * u.adu, data, (2**16 - 1) * u.adu) - data = data.astype('uint16') - # Data type conversion strips units so need to put them back manually - image = ccdproc.CCDData(data, wcs=noiseless.wcs, unit=u.adu) - image.header['EXPTIME'] = exp_time - image.header['DARKSUB'] = subtract_dark - - return image diff --git a/gunagala/camera.py b/gunagala/camera.py index a1964d3..ccee3e8 100644 --- a/gunagala/camera.py +++ b/gunagala/camera.py @@ -2,6 +2,7 @@ Cameras (stricly the image sensor subsystem, not including optics, optical filters, etc) """ import os +import numpy as np from astropy import units as u from astropy.table import Table @@ -37,7 +38,8 @@ class Camera: Pixel pitch. Square pixels are assumed. Resolution : astropy.units.Quantity Two element Quantity containing the number of pixels across - the image sensor in both horizontal & vertical directions. + the image sensor in both vertical & horizontal directions. + (y, x) read_noise astropy.units.Quantity Intrinsic noise of image sensor and readout electronics, in electrons/pixel units. @@ -55,6 +57,13 @@ class Camera: minimum_exposure : astropy.units.Quantity Length of the shortest exposure that the camera is able to take. + dark_current_dist : scipy.stats.rv_continuous, optional + A 'frozen' continuous random variable object that describes the distribution + of dark currents for the pixels in the image sensor. Used to create a `dark frame` + of uncorrelated dark current values. If not given no dark frame is created. + dark_current_seed: int, optional + Seed used to initialise the random number generator before creating the dark + frame. Set to a fixed value if you need to repeatedly generate the same dark frame. Attributes ---------- @@ -86,9 +95,23 @@ class Camera: Sequence of wavelengths from the QE data QE : astropy.units.Quantity Sequence of quantum efficiency values from the QE data. + dark_frame: astropy.units.Quantity or None + Array containing the dark current values for the pixels of the image sensor """ - def __init__(self, bit_depth, full_well, gain, bias, readout_time, pixel_size, resolution, - read_noise, dark_current, QE, minimum_exposure): + def __init__(self, + bit_depth, + full_well, + gain, + bias, + readout_time, + pixel_size, + resolution, + read_noise, + dark_current, + QE, + minimum_exposure, + dark_current_dist=None, + dark_current_seed=None): self.bit_depth = int(bit_depth) self.full_well = ensure_unit(full_well, u.electron / u.pixel) @@ -112,3 +135,44 @@ def __init__(self, bit_depth, full_well, gain, bias, readout_time, pixel_size, r self.wavelengths, self.QE = get_table_data(QE, column_names=('Wavelength', 'QE'), column_units=(u.nm, u.electron / u.photon)) + + # Generate dark frame + if dark_current_dist is not None: + try: + dark_current_dist.rvs + except AttributeError: + raise ValueError("dark_current_dist ({}) has no rvs() method!".format(dark_current_dist)) + self.dark_frame = self._make_dark_frame(dark_current_dist, dark_current_seed) + else: + self.dark_frame = None + + def _make_dark_frame(self, distribution, seed=None): + """ + Function to create a dark current 'image' in electrons per second per pixel. + + Creates an array of random, uncorrelated dark current values drawn from the + statistical distribution defined by the `distribution` parameter and returns + it as an astropy.units.Quantity. + + Parameters + ---------- + distribution : scipy.stats.rv_continuous + A 'frozen' continuous random variable object that describes the distribution + of dark current values for the pixels in the image sensor. + seed: int + Seed used to initialise the random number generator before creating the dark + frame. Set to a fixed value if you need to repeatedly generate the same dark + frame. + """ + if seed is not None: + # Initialise RNG + np.random.seed(seed) + + dark_frame = distribution.rvs(size=self.resolution.value.astype(int)) + dark_frame = dark_frame * u.electron / (u.second * u.pixel)\ + + if seed is not None: + # Re-initialise RNG with random seed + np.random.seed() + + return dark_frame diff --git a/gunagala/data/performance.yaml b/gunagala/data/performance.yaml index 5f1327f..3cf7519 100644 --- a/gunagala/data/performance.yaml +++ b/gunagala/data/performance.yaml @@ -57,7 +57,7 @@ cameras: resolution: 4656, 3520 # pixel read_noise: 2.5 # electron / pixel dark_current: 0.008 # electron / (second * pixel) - QE: ZWO_QE.csv + QE: ZWO_QE_asi1600.csv minimum_exposure: 0.000032 # second cis115: # e2v CIS115 'Sirius' back side illuminated CMOS image sensor diff --git a/gunagala/data/performance_data/ZWO_QE.csv b/gunagala/data/performance_data/ZWO_QE_asi1600.csv similarity index 100% rename from gunagala/data/performance_data/ZWO_QE.csv rename to gunagala/data/performance_data/ZWO_QE_asi1600.csv diff --git a/gunagala/imager.py b/gunagala/imager.py index 5a8143d..5eed8f1 100644 --- a/gunagala/imager.py +++ b/gunagala/imager.py @@ -6,23 +6,26 @@ import numpy as np from scipy.interpolate import interp1d +from scipy.stats import poisson, norm import matplotlib matplotlib.use('Agg') from matplotlib import pyplot as plt from astropy import constants as c from astropy import units as u -from astropy.table import Table -from astropy.wcs import WCS +from astropy.wcs import WCS, InvalidTransformError +from astropy.coordinates import SkyCoord +from astropy.nddata import CCDData from gunagala.optic import Optic from gunagala.optical_filter import Filter from gunagala.camera import Camera -from gunagala.psf import PSF, MoffatPSF +from gunagala.psf import PSF, MoffatPSF, FittablePSF from gunagala.sky import Sky, Simple, ZodiacalLight from gunagala.config import load_config from gunagala.utils import ensure_unit + def create_imagers(config=None): """ Parse config and create a corresponding dictionary of Imager objects. @@ -205,6 +208,7 @@ class Imager: Detected electrons/s/pixel due to the sky background for each filter bandpass. """ + def __init__(self, optic, camera, filters, psf, sky, num_imagers=1, num_per_computer=1): if not isinstance(optic, Optic): @@ -246,13 +250,13 @@ def __init__(self, optic, camera, filters, psf, sky, num_imagers=1, num_per_comp # Construct a simple template WCS to store the focal plane configuration parameters self.wcs = WCS(naxis=2) - self.wcs._naxis1 = self.camera.resolution[0].value - self.wcs._naxis2 = self.camera.resolution[1].value - self.wcs.wcs.crpix = [(self.camera.resolution[0].value + 1)/2, - (self.camera.resolution[1].value + 1)/2] + self.wcs._naxis2, self.wcs._naxis1 = self.camera.resolution.value.astype(int) + self.wcs.wcs.crpix = [(self.camera.resolution[1].value - 1) / 2, + (self.camera.resolution[0].value - 1) / 2] self.wcs.wcs.cdelt = [self.pixel_scale.to(u.degree / u.pixel).value, self.pixel_scale.to(u.degree / u.pixel).value] self.wcs.wcs.ctype = ['RA---TAN', 'DEC--TAN'] + self.wcs.wcs.crval = [None, None] # Calculate end to end efficiencies, etc. self._efficiencies() @@ -268,7 +272,8 @@ def __init__(self, optic, camera, filters, psf, sky, num_imagers=1, num_per_comp # Work out what *sort* of surface brightness we got and do something appropriate try: - surface_brightness = surface_brightness.to(u.photon / (u.second * u.m**2 * u.arcsecond**2 * u.nm)) + surface_brightness = surface_brightness.to( + u.photon / (u.second * u.m**2 * u.arcsecond**2 * u.nm)) except u.UnitConversionError: raise ValueError("I don't know what to do with this!") else: @@ -282,6 +287,27 @@ def __init__(self, optic, camera, filters, psf, sky, num_imagers=1, num_per_comp # Not a callable, should be a Simple sky model which just returns AB magnitudes per square arcsecond self.sky_rate[filter_name] = self.SB_to_rate(sb, filter_name) + def set_WCS_centre(self, centre, *args, **kwargs): + """ + Set the WCS CRVALs of the Imager instance to centre. + + Parameters + ---------- + centre : astropy.coordinates.SkyCoord or str + Sky coordinates of the image centre. Must be either a SkyCoord object or convertible + to one by the constructor of SkyCoord. + unit : str (optional) + Unit is the same as for the astropy.coordinates.SkyCoord(), e.g. + SkyCoord(RAdec[0], RAdec[1], unit='deg') + """ + + # Ensure centre is a SkyCoord (this allows entering centre as a string) + if not isinstance(centre, SkyCoord): + centre = SkyCoord(centre, **kwargs) + + # Set field centre coordinates in internal WCS + self.wcs.wcs.crval = [centre.icrs.ra.value, centre.icrs.dec.value] + def extended_source_signal_noise(self, surface_brightness, filter_name, total_exp_time, sub_exp_time, calc_type='per pixel', saturation_check=True, binning=1): """ @@ -334,7 +360,8 @@ def extended_source_signal_noise(self, surface_brightness, filter_name, total_ex raise ValueError("Invalid calculation type '{}'!".format(calc_type)) if calc_type == 'per arcsecond squared' and binning != 1: - raise ValueError("Cannot specify pixel binning with calculation type 'per arcsecond squared'!") + raise ValueError( + "Cannot specify pixel binning with calculation type 'per arcsecond squared'!") if surface_brightness: # Given a source brightness @@ -378,11 +405,13 @@ def extended_source_signal_noise(self, surface_brightness, filter_name, total_ex # Noise sources (per pixel for single imager) signal = (rate * total_exp_time).to(u.electron / u.pixel) # If calculating the signal & noise for the sky itself need to avoid double counting it here - sky_counts = self.sky_rate[filter_name] * total_exp_time if surface_brightness else 0 * u.electron / u.pixel + sky_counts = self.sky_rate[filter_name] * \ + total_exp_time if surface_brightness else 0 * u.electron / u.pixel dark_counts = self.camera.dark_current * total_exp_time total_read_noise = number_subs**0.5 * self.camera.read_noise - noise = ((signal + sky_counts + dark_counts) * (u.electron / u.pixel) + total_read_noise**2)**0.5 + noise = ((signal + sky_counts + dark_counts) * + (u.electron / u.pixel) + total_read_noise**2)**0.5 noise = noise.to(u.electron / u.pixel) # Saturation check @@ -391,7 +420,8 @@ def extended_source_signal_noise(self, surface_brightness, filter_name, total_ex saturated = self._is_saturated(rate, sub_exp_time, filter_name) else: # Sky counts already included in _is_saturated, need to avoid counting them twice - saturated = self._is_saturated(0 * u.electron / (u.pixel * u.second), sub_exp_time, filter_name) + saturated = self._is_saturated( + 0 * u.electron / (u.pixel * u.second), sub_exp_time, filter_name) # np.where strips units, need to manually put them back. signal = np.where(saturated, 0, signal) * u.electron / u.pixel noise = np.where(saturated, 0, noise) * u.electron / u.pixel @@ -509,7 +539,8 @@ def extended_source_etc(self, surface_brightness, filter_name, snr_target, sub_e raise ValueError("Invalid calculation type '{}'!".format(calc_type)) if calc_type == 'per arcsecond squared' and binning != 1: - raise ValueError("Cannot specify pixel binning with calculation type 'per arcsecond squared'!") + raise ValueError( + "Cannot specify pixel binning with calculation type 'per arcsecond squared'!") # Convert target SNR per array combined, binned pixel to SNR per unbinned pixel snr_target = ensure_unit(snr_target, u.dimensionless_unscaled) @@ -543,12 +574,12 @@ def extended_source_etc(self, surface_brightness, filter_name, snr_target, sub_e noise_squared_rate = ((rate + self.sky_rate[filter_name] + self.camera.dark_current) * (u.electron / u.pixel) + - self.camera.read_noise**2 / sub_exp_time) + self.camera.read_noise**2 / sub_exp_time) else: # Avoiding counting sky noise twice when the target is the sky background itself noise_squared_rate = ((rate + self.camera.dark_current) * (u.electron / u.pixel) + - self.camera.read_noise**2 / sub_exp_time) + self.camera.read_noise**2 / sub_exp_time) noise_squared_rate = noise_squared_rate.to(u.electron**2 / (u.pixel**2 * u.second)) total_exp_time = (snr_target**2 * noise_squared_rate / rate**2).to(u.second) @@ -564,7 +595,8 @@ def extended_source_etc(self, surface_brightness, filter_name, snr_target, sub_e saturated = self._is_saturated(rate, sub_exp_time, filter_name) else: # Sky counts already included in _is_saturated, need to avoid counting them twice - saturated = self._is_saturated(0 * u.electron / (u.pixel * u.second), sub_exp_time, filter_name) + saturated = self._is_saturated( + 0 * u.electron / (u.pixel * u.second), sub_exp_time, filter_name) number_subs = np.where(saturated, 0, number_subs) @@ -615,7 +647,8 @@ def extended_source_limit(self, total_exp_time, filter_name, snr_target, sub_exp raise ValueError("Invalid calculation type '{}'!".format(calc_type)) if calc_type == 'per arcsecond squared' and binning != 1: - raise ValueError("Cannot specify pixel binning with calculation type 'per arcsecond squared'!") + raise ValueError( + "Cannot specify pixel binning with calculation type 'per arcsecond squared'!") # Convert target SNR per array combined, binned pixel to SNR per unbinned pixel snr_target = ensure_unit(snr_target, u.dimensionless_unscaled) @@ -635,8 +668,10 @@ def extended_source_limit(self, total_exp_time, filter_name, snr_target, sub_exp total_exp_time = number_subs * sub_exp_time # Noise sources - sky_counts = self.sky_rate[filter_name] * total_exp_time if enable_sky_noise else 0.0 * u.electron / u.pixel - dark_counts = self.camera.dark_current * total_exp_time if enable_dark_noise else 0.0 * u.electron / u.pixel + sky_counts = self.sky_rate[filter_name] * \ + total_exp_time if enable_sky_noise else 0.0 * u.electron / u.pixel + dark_counts = self.camera.dark_current * \ + total_exp_time if enable_dark_noise else 0.0 * u.electron / u.pixel total_read_noise = number_subs**0.5 * \ self.camera.read_noise if enable_read_noise else 0.0 * u.electron / u.pixel @@ -645,7 +680,8 @@ def extended_source_limit(self, total_exp_time, filter_name, snr_target, sub_exp # Calculate science count rate for target signal to noise ratio a = total_exp_time**2 - b = -snr_target**2 * total_exp_time * u.electron / u.pixel # Units come from converting signal counts to noise + # Units come from converting signal counts to noise + b = -snr_target**2 * total_exp_time * u.electron / u.pixel c = -snr_target**2 * noise_squared rate = (-b + (b**2 - 4 * a * c)**0.5) / (2 * a) @@ -837,7 +873,7 @@ def flux_to_ABmag(self, flux, filter_name): # First convert from total flux to spectral flux densitity per # unit wavelength, using the pre-caluclated integral. - f_nu = flux * u.electron / (self._iminus2[filter_name] * c.c * u.photon) + f_nu = flux * u.electron / (self._iminus2[filter_name] * c.c * u.photon) # Then convert spectral flux density to magnitudes return f_nu.to(u.ABmag, @@ -866,7 +902,8 @@ def total_exposure_time(self, total_elapsed_time, sub_exp_time): total_elapsed_time = ensure_unit(total_elapsed_time, u.second) sub_exp_time = ensure_unit(sub_exp_time, u.second) - num_of_subs = np.floor(total_elapsed_time / (sub_exp_time + self.camera.readout_time * self.num_per_computer)) + num_of_subs = np.floor(total_elapsed_time / (sub_exp_time + + self.camera.readout_time * self.num_per_computer)) total_exposure_time = num_of_subs * sub_exp_time return total_exposure_time @@ -1315,14 +1352,16 @@ def exp_time_sequence(self, raise ValueError("This Imager has no filter '{}'!".format(filter_name)) if bool(bright_limit) == bool(shortest_exp_time): - raise ValueError("One and only one of bright_limit and shortest_exp_time must be specified!") + raise ValueError( + "One and only one of bright_limit and shortest_exp_time must be specified!") if bool(faint_limit) == bool(num_long_exp): raise ValueError("one and only one of faint_limit and num_long_exp must be specified!") longest_exp_time = ensure_unit(longest_exp_time, u.second) if longest_exp_time < self.camera.minimum_exposure: - raise ValueError("Longest exposure time shorter than minimum exposure time of the camera!") + raise ValueError( + "Longest exposure time shorter than minimum exposure time of the camera!") if bright_limit: # First calculate exposure time that will just saturate on the brightest sources. @@ -1346,7 +1385,8 @@ def exp_time_sequence(self, # Round down the shortest exposure time so that it is a exp_time_ratio^integer multiple of the longest # exposure time - num_exp_times = int(math.ceil(math.log(longest_exp_time / shortest_exp_time, exp_time_ratio))) + num_exp_times = int( + math.ceil(math.log(longest_exp_time / shortest_exp_time, exp_time_ratio))) shortest_exp_time = (longest_exp_time / (exp_time_ratio ** num_exp_times)) # Ensuring the shortest exposure time is not lower than the minimum exposure time of the cameras @@ -1357,7 +1397,8 @@ def exp_time_sequence(self, # Creating a list of exposure times from the shortest exposure time to the one directly below the # longest exposure time - exp_times = [shortest_exp_time.to(u.second) * exp_time_ratio**i for i in range(num_exp_times)] + exp_times = [shortest_exp_time.to(u.second) * exp_time_ratio ** + i for i in range(num_exp_times)] if faint_limit: num_long_exp = 0 @@ -1510,15 +1551,141 @@ def snr_vs_ABmag(self, exp_times, filter_name, magnitude_interval=0.02 * u.ABmag return magnitudes.to(u.ABmag), snrs.to(u.dimensionless_unscaled) + def get_pixel_coords(self): + """ + Utility function to return a SkyCoord array containing the on sky position + of the centre of all the pixels in the image. + """ + # Arrays of pixel coordinates + XY = np.meshgrid(np.arange(self.wcs._naxis1), np.arange(self.wcs._naxis2)) + + # Convert to arrays of RA, dec (ICRS, decimal degrees) + try: + RAdec = self.wcs.all_pix2world(XY[0], XY[1], 0) + except InvalidTransformError: + raise ValueError("CRVAL not set! Must call set_WCS_centre before get_pixel_coords") + + return SkyCoord(RAdec[0], RAdec[1], unit='deg') + + def make_noiseless_image(self, + centre, + obs_time, + filter_name, + centre_kwargs={}, + stars=None, + star_kwargs={}): + """ + Creates a noiseless simulated image for a given image centre and + observation time. + + Parameters + ---------- + centre : astropy.coordinates.SkyCoord or str + Sky coordinates of the image centre. Must be either a SkyCoord + object or convertible to one by the constructor of SkyCoord. + obs_time : astropy.time.Time or str + Time of the obseration. This can be relevant when calculating the + sky background and source positions. Must be either a Time object + or convertible to one by the constructor of Time. + filter_name : str + Name of the optical filter to use. + stars : sequence, optional + Sequence containing + centre_kwargs : dict (optional) + kwargs for centre kwargs to send to astropy.coordinates.SkyCoord(), + e.g. SkyCoord(centre, unit='deg') + star_kwargs : dict (optional) + kwargs for star kwargs to send to astropy.coordinates.SkyCoord(), + e.g. SkyCoord(coords, unit='deg') + + Returns + ------- + noiseless : astropy.nddata.CDDData + Noiseless image in the form of a CCDData object. + """ + + if isinstance(self.psf, FittablePSF): + raise NotImplementedError("Analytical PSFs currently don't work.\n They will take all system memory. See: https://github.com/AstroHuntsman/gunagala/pull/16#issuecomment-426844974 ") + + electrons = np.zeros((self.wcs._naxis2, + self.wcs._naxis1)) * u.electron / (u.second * u.pixel) + self.set_WCS_centre(centre, **centre_kwargs) + + # Calculate observed sky background + sky_rate = self.sky_rate[filter_name] + if hasattr(self.sky, 'relative_brightness'): + pixel_coords = self.get_pixel_coords() + relative_sky = self.sky.relative_brightness(pixel_coords, obs_time) + sky_rate = sky_rate * relative_sky + electrons = electrons + sky_rate + + if stars is not None: + for (coords, magnitude) in stars: + coords = SkyCoord(coords, **star_kwargs) + pixel_coords = self.wcs.all_world2pix(((coords.ra.degree, coords.dec.degree),), 0) \ + - self.wcs.wcs.crpix + star_rate = self.ABmag_to_rate(magnitude, filter_name) + star_image = star_rate * self.psf.pixellated(size=self.camera.resolution / u.pixel, + offsets=pixel_coords[0]) / u.pixel + electrons = electrons + star_image + + noiseless = CCDData(electrons, wcs=self.wcs) + + return noiseless + + def make_image_real(self, noiseless, exp_time, subtract_dark=False): + """ + Given a noiseless simulated image in electrons per pixel add dark current, + Poisson noise and read noise, and convert to ADU using the predefined gain. + """ + # Scale photoelectron rates by exposure time + data = noiseless.data * noiseless.unit * exp_time + # Add dark current + if self.camera.dark_frame is None: + data += self.camera.dark_current * exp_time + else: + data += self.camera.dark_frame * exp_time + # Force to electron units + data = (data * u.pixel).to(u.electron) + # Apply Poisson noise. This is not unit-aware, need to restore them manually + data = (poisson.rvs(data / u.electron)).astype('float64') * u.electron + # Apply read noise. Again need to restore units manually + data += norm.rvs(scale=self.camera.read_noise / (u.electron * u.pixel), + size=data.shape) * u.electron + # Optionally subtract a Perfect Dark + if subtract_dark: + if self.camera.dark_frame is None: + data -= self.camera.dark_current * exp_time * u.pixel + else: + data -= self.camera.dark_frame * exp_time * u.pixel + # Convert to ADU + data /= self.camera.gain + # Force to adu (just here to catch unit problems) + data = data.to(u.adu) + # 'Analogue to digital conversion' + data += self.camera.bias * u.pixel + data = np.where(data < 2**self.camera.bit_depth * u.adu, + data, + (2**self.camera.bit_depth - 1) * u.adu) + data = data.astype('uint16') + # Data type conversion strips units so need to put them back manually + image = CCDData(data, wcs=noiseless.wcs, unit=u.adu) + image.header['EXPTIME'] = exp_time.value + image.header['DARKSUB'] = subtract_dark + + return image + def _is_saturated(self, rate, sub_exp_time, filter_name, n_sigma=3.0): # Total electrons per pixel from source, sky and dark current - electrons_per_pixel = (rate + self.sky_rate[filter_name] + self.camera.dark_current) * sub_exp_time + electrons_per_pixel = ( + rate + self.sky_rate[filter_name] + self.camera.dark_current) * sub_exp_time # Consider saturated if electrons per pixel is closer than n sigmas of noise to the saturation level return electrons_per_pixel > self.camera.saturation_level - n_sigma * self.camera.max_noise def _efficiencies(self): # Fine wavelength grid spaning maximum range of instrument response - waves = np.arange(self.camera.wavelengths.value.min(), self.camera.wavelengths.value.max(), 0.05) * u.nm + waves = np.arange(self.camera.wavelengths.value.min(), + self.camera.wavelengths.value.max(), 0.05) * u.nm self.wavelengths = waves # Sensitivity integrals for each filter bandpass @@ -1536,8 +1703,10 @@ def _efficiencies(self): self.bandwidth = {} # Interpolators for throughput and QE. Will move these into the Optics and Camera classes later. - tau = interp1d(self.optic.wavelengths, self.optic.throughput, kind='linear', fill_value='extrapolate') - qe = interp1d(self.camera.wavelengths, self.camera.QE, kind='linear', fill_value='extrapolate') + tau = interp1d(self.optic.wavelengths, self.optic.throughput, + kind='linear', fill_value='extrapolate') + qe = interp1d(self.camera.wavelengths, self.camera.QE, + kind='linear', fill_value='extrapolate') for name, band in self.filters.items(): # End-to-end efficiency. Need to put units back after interpolation @@ -1568,6 +1737,7 @@ def _gamma0(self): # Average photon energy energy = c.h * c.c / (self.pivot_wave * u.photon) # Divide by photon energy & multiply by aperture area, pixel area and bandwidth to get photons/s/pixel - photon_flux = (sfd_sb_0 / energy) * self.optic.aperture_area * self.pixel_area * self.bandwidth + photon_flux = (sfd_sb_0 / energy) * self.optic.aperture_area * \ + self.pixel_area * self.bandwidth self.gamma0 = photon_flux.to(u.photon / (u.s * u.pixel)) diff --git a/gunagala/psf.py b/gunagala/psf.py index 34fe72c..acbe974 100644 --- a/gunagala/psf.py +++ b/gunagala/psf.py @@ -12,7 +12,7 @@ from astropy.modeling.functional_models import Moffat2D from gunagala import utils - +from warnings import warn class PSF(): """ @@ -97,17 +97,17 @@ def _get_peak(self): saturation limits for point sources. """ # Odd number of pixels (1) so offsets = (0, 0) is centred on a pixel - centred_psf = self.pixellated(size=1, offsets=(0, 0)) + centred_psf = self.pixellated(size=(1, 1), offsets=(0, 0)) return centred_psf[0, 0] / u.pixel - def _get_n_pix(self, size=20): + def _get_n_pix(self, size=(20, 20)): """ Calculate the effective number of pixels for PSF fitting photometry with this PSF, in the worse case where the PSF is centred on the corner of a pixel. """ # Want a even number of pixels. - size = size + size % 2 + size = tuple(s + s % 2 for s in size) # Even number of pixels so offsets = (0, 0) is centred on pixel corners corner_psf = self.pixellated(size, offsets=(0, 0)) return 1 / ((corner_psf**2).sum()) * u.pixel @@ -166,7 +166,7 @@ def FWHM(self, FWHM): if self.pixel_scale: self._update_model() - def pixellated(self, size=21, offsets=(0.0, 0.0)): + def pixellated(self, size=(21, 21), offsets=(0.0, 0.0)): """ Calculates a pixellated version of the PSF. @@ -176,9 +176,8 @@ def pixellated(self, size=21, offsets=(0.0, 0.0)): Parameters ---------- - size : int, optional - Size of the pixellated PSF to calculate, the returned image - will have `size` x `size` pixels. Default value 21. + size : (int, int) optional + y, x size of the pixellated PSF to calculate. Default value (21, 21). offset : tuple of floats, optional y and x axis offsets of the centre of the PSF from the centre of the returned image, in pixels. @@ -191,16 +190,16 @@ def pixellated(self, size=21, offsets=(0.0, 0.0)): pixellated PSF will be somewhat less due to truncation of the PSF wings by the edge of the image. """ - size = int(size) - if size <= 0: + size = tuple(int(s) for s in size) + if size[0] <= 0 or size[1] <=0: raise ValueError("`size` must be > 0, got {}!".format(size)) # Update PSF centre coordinates - self.x_0 = offsets[0] - self.y_0 = offsets[1] + self.x_0 = offsets[1] + self.y_0 = offsets[0] - xrange = (-(size - 1) / 2, (size + 1) / 2) - yrange = (-(size - 1) / 2, (size + 1) / 2) + xrange = (-(size[1] - 1) / 2, (size[1] + 1) / 2) + yrange = (-(size[0] - 1) / 2, (size[0] + 1) / 2) return discretize_model(self, xrange, yrange, mode='oversample', factor=10) @@ -326,7 +325,7 @@ def __init__(self, self._psf_data = psf_data self._psf_sampling = utils.ensure_unit(psf_sampling, u.arcsecond / u.pixel) if psf_centre is None: - self._psf_centre = np.array(psf_data.shape) / 2 + self._psf_centre = (np.array(psf_data.shape) - 1) / 2 else: self._psf_centre = np.array(psf_centre) self._oversampling = int(oversampling) @@ -336,7 +335,7 @@ def __init__(self, self.pixel_scale = pixel_scale - def pixellated(self, size=21, offsets=(0.0, 0.0)): + def pixellated(self, size=(21, 21), offsets=(0.0, 0.0)): """ Calculates a pixellated version of the PSF. @@ -345,9 +344,8 @@ def pixellated(self, size=21, offsets=(0.0, 0.0)): Parameters ---------- - size : int, optional - Size of the pixellated PSF to calculate, the returned image - will have `size` x `size` pixels. Default value 21. + size : (int, int), optional + y, x size of the pixellated PSF to calculate. Default value (21, 21). offsets : tuple of floats, optional y and x axis offsets of the centre of the PSF from the centre of the returned image, in pixels. @@ -355,37 +353,87 @@ def pixellated(self, size=21, offsets=(0.0, 0.0)): Returns ------- pixellated : numpy.array - Pixellated PSF image with `size` by `size` pixels. The PSF + Pixellated PSF image with `size[0]` by `size[1]` pixels. The PSF is normalised to an integral of 1 however the sum of the pixellated PSF will be somewhat less due to truncation of the PSF wings by the edge of the image. """ - resampled_size = int(size) * self._oversampling - # Arrays of pixel coordinates relative to array centre - resampled_coordinates = np.mgrid[-(resampled_size - 1) / 2:resampled_size / 2, - -(resampled_size - 1) / 2:resampled_size / 2] - # Apply offsets so origin is at desired PSF centre - resampled_offsets = np.reshape(np.array(offsets) * self._oversampling, (2, 1, 1)) - resampled_coordinates = resampled_coordinates - resampled_offsets - # Convert from resampled PSF pixel units to original PSF pixel units - resampled_coordinates = resampled_coordinates / self._resampling_factor + size = np.array(size, dtype=np.int) + offsets = np.array(offsets) + # Only want to caclulate resampled PSF for positions that fall within the PSF data, + # otherwise end up filling the RAM with lots of double precision zeros. + + # Initialise output array + pixellated = np.zeros(size) + + # Limits of psf_data footprint in its own pixel coordinates + limits = np.array(((-0.5, -0.5), + (self._psf_data.shape[0] - 0.5, self._psf_data.shape[1] - 0.5))) + # Subtract position of PSF centre so origin is at the PSF centre + limits = limits - self._psf_centre + # Convert from psf_data pixel units to output array pixel units + limits = limits * (self._psf_sampling / self.pixel_scale).to(u.dimensionless_unscaled).value + # Reverse offsets so origin is at output array centre + limits = limits + np.array(offsets) + # Move origin to output array origin + limits = limits + (size - 1) / 2 + # Round limits to the centres of the pixels containing the boundary + limits = np.rint(limits).astype(np.int) + # Crop to output array edges + limits = np.array((np.where(limits[0] >= 0, limits[0], 0), + np.where(limits[1] < size, limits[1], size - 1))) + # Store output array limits for later + y0 = limits[0, 0] + y1 = limits[1, 0] + 1 + x0 = limits[0, 1] + x1 = limits[1, 1] + 1 + # Origin back to output array centre + limits = limits - (size - 1) / 2 + # Expand to pixel edges + limits = np.array((limits[0] - 0.5, limits[1] + 0.5)) + # Convert from output array pixels to oversampled array pixels + limits = limits * self._oversampling + # Contract by half a pixel to align with oversampled pixel centres + limits = np.array((limits[0] + 0.5, limits[1] - 0.5)) + # Apply offset so origin is at desired PSF centre + limits = limits - offsets * self._oversampling + # Convert from resampled PSF pixel units to psf_data pixel units + limits = limits / self._resampling_factor # Add position of PSF centre in psf_data so origin is the same as the origin of psf_data - resampled_coordinates = resampled_coordinates + np.reshape(self._psf_centre, (2, 1, 1)) + limits = limits + self._psf_centre + # Arrays of coordinates relative to output array centre. The half steps are to avoid + # problems with floating point precision & the stopping condition. + step = 1 / self._resampling_factor + resampled_coordinates = np.mgrid[limits[0, 0]:limits[1, 0] + step / 2:step, + limits[0, 1]:limits[1, 1] + step / 2:step] # Calculate resampled PSF using cubic spline interpolation resampled_psf = ndimage.map_coordinates(self._psf_data, resampled_coordinates) - # Rebin to the correct pixel scale - pixellated = utils.bin_array(resampled_psf, self._oversampling) + # Rebin to the output array pixel scale + resampled_psf = utils.bin_array(resampled_psf, self._oversampling) # Renormalise to correct for the effect of resampling - return pixellated / self._resampling_factor**2 + resampled_psf = resampled_psf / self._resampling_factor**2 + # Check and if there are some below zero pixel values, clip them to zero and renormalize resampled_psf. + if (resampled_psf < 0).any(): + warn("Warning: below zero values in resampled PSF. Clipping to zero.") + previous_total = resampled_psf.sum() + resampled_psf[resampled_psf < 0] = 0 + resampled_psf *= previous_total / resampled_psf.sum() + # Insert into output array in the correct place. + pixellated[y0:y1,x0:x1] = resampled_psf + return pixellated def _get_n_pix(self): # For accurate results want the calculation to include the whole PSF. - size = int(math.ceil(max(self._psf_data.shape) * self._psf_sampling / self.pixel_scale)) + psf_data_size = np.array(self._psf_data.shape) + size = psf_data_size + np.abs(self._psf_centre - psf_data_size / 2) + size = size * self._psf_sampling / self.pixel_scale + size = np.ceil(size + 0.5) return super()._get_n_pix(size=size) def _update_model(self): self._resampled_scale = self.pixel_scale / self._oversampling - self._resampling_factor = self._psf_sampling / self._resampled_scale + self._resampling_factor = (self._psf_sampling / self._resampled_scale) + self._resampling_factor = self._resampling_factor.to(u.dimensionless_unscaled).value self._n_pix = self._get_n_pix() self._peak = self._get_peak() diff --git a/gunagala/tests/test_camera.py b/gunagala/tests/test_camera.py index 0fd9b8a..90a00ed 100644 --- a/gunagala/tests/test_camera.py +++ b/gunagala/tests/test_camera.py @@ -1,4 +1,5 @@ import pytest +from scipy import stats import astropy.units as u from gunagala.camera import Camera @@ -31,3 +32,95 @@ def test_saturation_level(ccd): def test_max_noise(ccd): assert ccd.max_noise == (ccd.saturation_level * u.electron / u.pixel + (9.3 * u.electron / u.pixel)**2)**0.5 + + +def test_dark_frame(): + shape = 0.5 + loc = 0.4 + scale = 0.02 + dist = stats.lognorm(s=shape, loc=loc, scale=scale) + ccd = Camera(bit_depth=16, + full_well=25500 * u.electron / u.pixel, + gain=0.37 * u.electron / u.adu, + bias=1100 * u.adu / u.pixel, + readout_time=0.9 * u.second, + pixel_size=5.4 * u.micron / u.pixel, + resolution=(3326, 2504) * u.pixel, + read_noise=9.3 * u.electron / u.pixel, + dark_current=0.04 * u.electron / (u.pixel * u.second), + minimum_exposure=0.1 * u.second, + QE='ML8300M_QE.csv', + dark_current_dist=dist, + dark_current_seed=42) + assert isinstance(ccd.dark_frame, u.Quantity) + assert (ccd.dark_frame.shape * u.pixel == ccd.resolution).all() + fitted_params = stats.lognorm.fit(ccd.dark_frame.value[0:100,0:100]) + assert fitted_params[0] == pytest.approx(shape, rel=0.02) + assert fitted_params[1] == pytest.approx(loc, rel=0.02) + assert fitted_params[2] == pytest.approx(scale, rel=0.02) + + +def test_dark_seed(): + shape = 0.5 + loc = 0.4 + scale = 0.02 + dist = stats.lognorm(s=shape, loc=loc, scale=scale) + ccd1 = Camera(bit_depth=16, + full_well=25500 * u.electron / u.pixel, + gain=0.37 * u.electron / u.adu, + bias=1100 * u.adu / u.pixel, + readout_time=0.9 * u.second, + pixel_size=5.4 * u.micron / u.pixel, + resolution=(3326, 2504) * u.pixel, + read_noise=9.3 * u.electron / u.pixel, + dark_current=0.04 * u.electron / (u.pixel * u.second), + minimum_exposure=0.1 * u.second, + QE='ML8300M_QE.csv', + dark_current_dist=dist, + dark_current_seed=42) + ccd2 = Camera(bit_depth=16, + full_well=25500 * u.electron / u.pixel, + gain=0.37 * u.electron / u.adu, + bias=1100 * u.adu / u.pixel, + readout_time=0.9 * u.second, + pixel_size=5.4 * u.micron / u.pixel, + resolution=(3326, 2504) * u.pixel, + read_noise=9.3 * u.electron / u.pixel, + dark_current=0.04 * u.electron / (u.pixel * u.second), + minimum_exposure=0.1 * u.second, + QE='ML8300M_QE.csv', + dark_current_dist=dist, + dark_current_seed=42) + assert (ccd1.dark_frame == ccd2.dark_frame).all() + + +def test_dark_no_seed(): + shape = 0.5 + loc = 0.4 + scale = 0.02 + dist = stats.lognorm(s=shape, loc=loc, scale=scale) + ccd1 = Camera(bit_depth=16, + full_well=25500 * u.electron / u.pixel, + gain=0.37 * u.electron / u.adu, + bias=1100 * u.adu / u.pixel, + readout_time=0.9 * u.second, + pixel_size=5.4 * u.micron / u.pixel, + resolution=(3326, 2504) * u.pixel, + read_noise=9.3 * u.electron / u.pixel, + dark_current=0.04 * u.electron / (u.pixel * u.second), + minimum_exposure=0.1 * u.second, + QE='ML8300M_QE.csv', + dark_current_dist=dist) + ccd2 = Camera(bit_depth=16, + full_well=25500 * u.electron / u.pixel, + gain=0.37 * u.electron / u.adu, + bias=1100 * u.adu / u.pixel, + readout_time=0.9 * u.second, + pixel_size=5.4 * u.micron / u.pixel, + resolution=(3326, 2504) * u.pixel, + read_noise=9.3 * u.electron / u.pixel, + dark_current=0.04 * u.electron / (u.pixel * u.second), + minimum_exposure=0.1 * u.second, + QE='ML8300M_QE.csv', + dark_current_dist=dist) + assert (ccd1.dark_frame != ccd2.dark_frame).all() diff --git a/gunagala/tests/test_imager.py b/gunagala/tests/test_imager.py index 5d9cb81..ccfe25e 100644 --- a/gunagala/tests/test_imager.py +++ b/gunagala/tests/test_imager.py @@ -3,6 +3,7 @@ import numpy as np import astropy.units as u import astropy.constants as c +from astropy.coordinates import SkyCoord from gunagala.optic import Optic from gunagala.optical_filter import Filter @@ -77,6 +78,12 @@ def imager(lens, ccd, filters, psf, sky): return imager +@pytest.fixture(scope='function') +def imager_function_scope(lens, ccd, filters, psf, sky): + imager = Imager(optic=lens, camera=ccd, filters=filters, psf=psf, sky=sky, num_imagers=5, num_per_computer=5) + return imager + + def test_init(imager): assert isinstance(imager, Imager) assert imager.pixel_scale == (5.4 * u.micron / (391 * u.mm * u.pixel)).to(u.arcsecond / u.pixel, @@ -735,6 +742,22 @@ def test_snr_vs_mag(imager, filter_name, tmpdir): assert mags3[-1].value == pytest.approx(mags[-1].value - 2.5, abs=0.1) +def test_get_pixel_coords_no_WCS_call_first(imager_function_scope): + # test if get_pixel_coords is called without first running set_WCS_centre + with pytest.raises(ValueError): + imager_function_scope.get_pixel_coords() + + +def test_get_pixel_coords_with_WCS_call_first(imager_function_scope): + test_coord_string = "189.9976325 -11.6230544" + + # now set WCS centre first, then try and get_pixel_coords + imager_function_scope.set_WCS_centre(test_coord_string, unit='deg') + centre_field_pixels = imager_function_scope.get_pixel_coords() + assert imager_function_scope.wcs._naxis1 == centre_field_pixels.shape[1] + assert imager_function_scope.wcs._naxis2 == centre_field_pixels.shape[0] + + def test_create_imagers(): imagers = create_imagers() assert isinstance(imagers, dict) diff --git a/gunagala/tests/test_psf.py b/gunagala/tests/test_psf.py index efc82e1..e932c54 100644 --- a/gunagala/tests/test_psf.py +++ b/gunagala/tests/test_psf.py @@ -6,13 +6,13 @@ @pytest.fixture(scope='module') -def psf(): +def psf_moffat(): psf = MoffatPSF(FWHM=1 / 30 * u.arcminute, shape=4.7) return psf @pytest.fixture(scope='module') -def pix_psf(): +def psf_pixellated(): psf_data = np.array([[0.0, 0.0, 0.1, 0.0, 0.0], [0.0, 0.3, 0.7, 0.4, 0.0], [0.1, 0.8, 1.0, 0.6, 0.1], @@ -32,8 +32,13 @@ def test_base(): psf_base = PSF(FWHM=1 / 30 * u.arcminute) -def test_moffat(psf): - assert isinstance(psf, MoffatPSF) +@pytest.mark.parametrize("psf, type", [ + (psf_moffat(), MoffatPSF), + (psf_pixellated(), PixellatedPSF)], + ids=["moffat", "pixellated"] +) +def test_instance(psf, type): + assert isinstance(psf, type) assert isinstance(psf, PSF) @@ -83,33 +88,94 @@ def test_shape(psf): psf.shape = 2.5 assert psf.shape == 2.5 with pytest.raises(ValueError): - psf.shape = 0.5 - psf.shape = 4.7 - - -def test_pixellated(psf): - pixellated = psf.pixellated() - assert isinstance(pixellated, np.ndarray) - assert pixellated.shape == (21, 21) - pixellated = psf.pixellated(size=7.2) - assert isinstance(pixellated, np.ndarray) - assert pixellated.shape == (7, 7) - pixellated = psf.pixellated(offsets=(0.3, -0.7)) - assert isinstance(pixellated, np.ndarray) - assert pixellated.shape == (21, 21) + psf_moffat.FWHM = -1 * u.degree + psf_moffat.FWHM = 2 * u.arcsecond + + +@pytest.mark.parametrize("psf, t_pixel_scale, pixel_scale", [ + (psf_moffat(), 2.85, 2.85), + (psf_pixellated(), (1 / 3), (2 / 3))], + ids=["moffat", "pixellated"] +) +def test_pixel_scale(psf, t_pixel_scale, pixel_scale): + psf.pixel_scale = t_pixel_scale * u.arcsecond / u.pixel + assert psf.pixel_scale == t_pixel_scale * u.arcsecond / u.pixel + psf.pixel_scale = pixel_scale * u.arcsecond / u.pixel + + +@pytest.mark.parametrize("psf, expected_n_pix, pixel_scale", [ + (psf_moffat(), 4.25754067000986, 2.85), + (psf_pixellated(), 21.06994544, (2 / 3))], + ids=["moffat", "pixellated"] +) +def test_n_pix(psf, expected_n_pix, pixel_scale): + psf.pixel_scale = pixel_scale * u.arcsecond / u.pixel + assert psf.n_pix.to(u.pixel).value == pytest.approx(expected_n_pix) + + +@pytest.mark.parametrize("psf, expected_peak, pixel_scale", [ + (psf_moffat(), 0.7134084656751443, 2.85), + (psf_pixellated(), 0.08073066, (2 / 3))], + ids=["moffat", "pixellated"] +) +def test_peak(psf, expected_peak, pixel_scale): + psf.pixel_scale = pixel_scale * u.arcsecond / u.pixel + assert psf.peak.to(1 / (u.pixel)).value == pytest.approx(expected_peak) + + +def test_shape(psf_moffat): + assert psf_moffat.shape == 4.7 + psf_moffat.shape = 2.5 + assert psf_moffat.shape == 2.5 with pytest.raises(ValueError): - psf.pixellated(size=-1.3) - - -def test_pixellated(pix_psf): - pixellated = pix_psf.pixellated() - assert isinstance(pixellated, np.ndarray) - assert pixellated.shape == (21, 21) - pixellated = pix_psf.pixellated(size=7.2) - assert isinstance(pixellated, np.ndarray) - assert pixellated.shape == (7, 7) - pixellated = pix_psf.pixellated(offsets=(0.3, -0.7)) - assert isinstance(pixellated, np.ndarray) - assert pixellated.shape == (21, 21) + psf_moffat.shape = 0.5 + psf_moffat.shape = 4.7 + + +@pytest.mark.parametrize("psf, image_size", [ + (psf_moffat(), (21, 21)), + (psf_pixellated(), (21, 21)), + (psf_moffat(), (7, 9)), + (psf_pixellated(), (7, 9))], + ids=["moffat_square", + "pixellated_square", + "moffat_rectangle", + "pixellated_rectangle"] +) +def test_pixellated_dimension(psf, image_size): + assert isinstance(psf.pixellated(), np.ndarray) + assert isinstance(psf.pixellated(size=( + image_size[0] + 0.2, image_size[1] + 0.2)), np.ndarray) + assert psf.pixellated(size=( + image_size[0] + 0.2, image_size[1] + 0.2)).shape == image_size + assert (psf.pixellated(size=( + image_size[0] + 0.2, image_size[1] + 0.2)) >= 0).all() + assert np.isfinite(psf.pixellated(size=( + image_size[0] + 0.2, image_size[1] + 0.2))).all() + + +@pytest.mark.parametrize("psf, offset", [ + (psf_moffat(), (0.0, 0.0)), + (psf_pixellated(), (0.0, 0.0)), + (psf_moffat(), (0.3, -0.7)), + (psf_pixellated(), (0.3, -0.7))], + ids=["moffat_centre_offsets", + "pixellated_centre_offsets", + "moffat_noncentre_offsets", + "pixellated_noncentre_offsets"] +) +def test_offsets(psf, offset): + assert isinstance(psf.pixellated(offsets=offset), np.ndarray) + assert psf.pixellated().shape == (21, 21) + assert (psf.pixellated() >= 0).all() + assert np.isfinite(psf.pixellated()).all() + + +@pytest.mark.parametrize("psf, test_size", [ + (psf_moffat(), (1.3, -1.3)), + (psf_pixellated(), (-1.3, 1.3))], + ids=["moffat", "pixellated"] +) +def test_pixellated_invalid_size(psf, test_size): with pytest.raises(ValueError): - pix_psf.pixellated(size=-1.3) + psf.pixellated(size=test_size)