diff --git a/stips/astro_image/astro_image.py b/stips/astro_image/astro_image.py index e1d8939..c207237 100644 --- a/stips/astro_image/astro_image.py +++ b/stips/astro_image/astro_image.py @@ -611,13 +611,14 @@ def makePoint(id, x, y, rate, psf_params, overrides, logger): """ ix, iy = int(np.floor(x)), int(np.floor(y)) if psf_params['has_psf']: - psf_img = getReducedPsf(psf_params, rate, x, y) - residual_ipc = SelectParameter("residual_ipc", overrides) - psf_type = SelectParameter("psf_type", overrides) - if residual_ipc and psf_type == "orig": - psf_img = apply_ipc(psf_img, 1) + psf_img = PSFGrid.reduced_psf_from_fits(psf_params['psf_file'], + psf_params['psf_ext'], + rate, + x, + y) else: - psf_img = np.array([[0., 0., 0.], [0., rate, 0.], [0., 0., 0.]], dtype=np.float64) + psf_img = np.array([[0., 0., 0.], [0., rate, 0.], [0., 0., 0.]], + dtype=np.float64) return psf_img,ix,iy @@ -630,8 +631,7 @@ def addPoints(self, ids, xs, ys, rates, *args, **kwargs): self._log("info", msg.format(len(xs),self.name)) psf_file = os.path.join(self.psf_dir, self.psf_name) - psf_params = {'psf_file': psf_file, 'psf_oversample': self.psf_oversample, - 'has_psf': self.has_psf} + psf_params = {'psf_file': psf_file, 'has_psf': self.has_psf} psf_type = SelectParameter('psf_type', self.overrides) if SelectParameter('residual_ipc', self.overrides): if 'ipc' not in psf_type: @@ -842,7 +842,11 @@ def galsimSersic(id, gal_params, psf_params, xsize, ysize, dir, overrides, logge central_flux = img[iyc, ixc] if psf_params['has_psf']: - psf_img = getReducedPsf(psf_params, 1., gal_params['x'], gal_params['y']) + psf_img = PSFGrid.reduced_psf_from_fits(psf_params['psf_file'], + psf_params['psf_ext'], + 1., + gal_params['x'], + gal_params['y']) result = convolve_fft(img, psf_img, allow_huge=True) else: result = img @@ -951,7 +955,11 @@ def astropySersic(id, gal_params, psf_params, xsize, ysize, dir, overrides, logg img *= flux / np.sum(img) if psf_params['has_psf']: - psf_img = getReducedPsf(psf_params, 1., gal_params['x'], gal_params['y']) + psf_img = PSFGrid.reduced_psf_from_fits(psf_params['psf_file'], + psf_params['psf_ext'], + 1., + gal_params['x'], + gal_params['y']) result = convolve_fft(img, psf_img, allow_huge=True) else: result = img @@ -1057,7 +1065,11 @@ def pandeiaSersic(galnum, gal_params, psf_params, xsize, ysize, dir, overrides, if psf_params['has_psf']: # print("\tGenerating PSF ({})".format(time.ctime())) - psf_img = getReducedPsf(psf_params, 1., gal_params['x'], gal_params['y']) + psf_img = PSFGrid.reduced_psf_from_fits(psf_params['psf_file'], + psf_params['psf_ext'], + 1., + gal_params['x'], + gal_params['y']) # print("\tPSF image has size {} dtype {}".format(psf_img.shape, psf_img.dtype)) # print("\tConvolving by PSF ({})".format(time.ctime())) result = convolve_fft(img, psf_img) @@ -1129,8 +1141,8 @@ def addSersics(self, ids, xs, ys, fluxes, ns, res, phis, ratios, *args, **kwargs central_fluxes = [] psf_file = os.path.join(self.psf_dir, self.psf_name) - psf_params = {'psf_file': psf_file, 'psf_oversample': self.psf_oversample, - 'has_psf': self.has_psf, 'psf_ext': self.psf_ext} + psf_params = {'psf_file': psf_file, 'has_psf': self.has_psf, + 'psf_ext': self.psf_ext} if self.parallel_enable: self._log("info", "Adding sersic profiles in parallel.") diff --git a/stips/utilities/PSFGrid.py b/stips/utilities/PSFGrid.py index fe6daf5..2458984 100644 --- a/stips/utilities/PSFGrid.py +++ b/stips/utilities/PSFGrid.py @@ -15,9 +15,11 @@ """ +import logging import numpy as np import os import shutil +import sys import webbpsf from collections import defaultdict @@ -50,7 +52,7 @@ def make_epsf_kernel(oversample): if oversample%2 == 0: kernel_size += 1 - kernel = np.ones((kernel_size, kernel_size), dtype=np.int64) + kernel = np.ones((kernel_size, kernel_size), dtype=np.float64) if oversample%2 == 0: kernel[0,:] = kernel[0,:] * 0.5 @@ -116,6 +118,15 @@ def __init__(self, data, **kwargs): Slight override to the parent class __init__ method, adding in some PSF information and parameters to the metadata attribute. """ + # if self.logger.level == logging.DEBUG: + self.logger = logging.getLogger('__stips__') + log_level = SelectParameter("log_level") + self.logger.setLevel(getattr(logging, log_level)) + if not len(self.logger.handlers): + stream_handler = logging.StreamHandler(sys.stderr) + format = '%(asctime)s %(levelname)s: %(message)s' + stream_handler.setFormatter(logging.Formatter(format)) + self.logger.addHandler(stream_handler) self.data = {} if 'meta' not in kwargs: meta = { @@ -148,7 +159,6 @@ def create_from_parameters(cls, **kwargs): - oversample - grid size """ - logger = kwargs.get('logger', None) data = {} meta = {} meta['psf_fov_pixels'] = SelectParameter('psf_fov_pixels', kwargs) @@ -158,7 +168,6 @@ def create_from_parameters(cls, **kwargs): meta['instrument'] = kwargs.get('instrument', 'wfi').lower() meta['detector'] = kwargs.get('detector', 'sca01').lower() meta['filter'] = kwargs.get('filter', 'f067').lower() - oversample_by_scale = kwargs.get('oversample_by_scale', True) pixel_scale = kwargs.get('pixel_scale', 0.11) psf_grid_size = SelectParameter('psf_grid_default_size', kwargs) num_psfs = psf_grid_size*psf_grid_size @@ -174,28 +183,21 @@ def create_from_parameters(cls, **kwargs): psf_xgrid = np.linspace(pix_limit['low'], pix_limit['high'], psf_grid_size) psf_ygrid = np.linspace(pix_limit['low'], pix_limit['high'], psf_grid_size) i = 0 - for xp in psf_xgrid: - for yp in psf_ygrid: - if logger is not None: - logger.info("Calculating PSF at ({},{})".format(xp, yp)) + for yp in psf_ygrid: + for xp in psf_xgrid: + self.logger.info("Calculating PSF at ({},{})".format(xp, yp)) meta['grid_xypos'].append((xp, yp)) ins = cls.get_instrument(meta['telescope'], meta['instrument']) ins.options['parity'] = 'odd' ins.filter = meta['filter'] ins.detector = meta['detector'] ins.detector_position = (xp, yp) - if oversample_by_scale: - ins.pixelscale = pixel_scale/meta['oversampling'] - psf_fits = ins.calc_psf(fov_pixels=fov, - nlambda=meta['psf_nlambda'], - oversample=1) - psf = psf_fits["DET_SAMP"].data - else: - ins.pixelscale = pixel_scale - psf_fits = ins.calc_psf(fov_pixels=meta['psf_fov_pixels'], - nlambda=meta['psf_nlambda'], - oversample=meta['oversampling']) - psf = psf_fits['OVERSAMP'].data + ins.pixelscale = pixel_scale/meta['oversampling'] + psf_fits = ins.calc_psf(fov_pixels=fov, + nlambda=meta['psf_nlambda'], + oversample=1) + psf = psf_fits["DET_SAMP"].data + self.logger.info("\tPSF size is {}".format(psf.shape)) data['orig'][i][:] = deepcopy(psf) * meta['oversampling']**2 data['orig_ipc'][i][:] = deepcopy(data['orig'][i][:]) data['epsf'][i][:] = make_epsf(deepcopy(psf), meta['oversampling']) @@ -243,6 +245,16 @@ def create_from_fits(cls, file_name): i += 1 data[psf_type] = GriddedPSFModel(NDData(psf_data, meta=meta)) return cls(data, meta=meta) + + + @classmethod + def reduced_psf_from_fits(cls, file_name, psf_ext, flux, cx, cy): + """ + Instantiates from a FITS file, then loads and bins a PSF at a given centre, then + softly and silently vanishes away. + """ + model = cls.create_from_fits(file_name) + return model.get_reduced_psf(psf_ext, flux, cx, cy) def from_fits(self, file_name): @@ -279,22 +291,6 @@ def to_fits(self, file_name): hdul.writeto(file_name, overwrite=True) - @staticmethod - def pix_limit(telescope, instrument, detector): - roman_det = defaultdict(lambda: {'low': 4, 'high': 4092}) - roman_ins = defaultdict(lambda: roman_det) - miri_det = defaultdict(lambda: {'low': 4, 'high': 1024}) - nircam_det = defaultdict(lambda: {'low': 0, 'high': 2040}) - detector_dict = { - 'roman': roman_ins, - 'jwst': { - 'nircam': nircam_det, - 'miri': miri_det - } - } - return detector_dict[telescope][instrument][detector] - - def make_gridded_model(self, data, **kwargs): """ Given data (which may be a GriddedPSFModel, an astropy NDData, or a numpy array), @@ -323,7 +319,74 @@ def make_gridded_model(self, data, **kwargs): for yp in psf_ygrid: meta['grid_xypos'].append((xp,yp)) return GriddedPSFModel(NDData(data, meta=meta), **kwargs) + + + def make_bin(self, x_size, y_size): + """ + Create X and Y lists that will result in creating lists of centre pixels from + which to create the lovely binned PSF model. + """ + oversample = self.meta['oversampling'] + xl = np.arange(x_size//2, oversample//2, -oversample) + xh = np.arange(x_size//2+oversample, x_size-oversample//2, oversample) + x_arr = list(np.append(np.sort(xl), np.sort(xh))) + yl = np.arange(y_size//2, oversample//2, -oversample) + yh = np.arange(y_size//2+oversample, y_size-oversample//2, oversample) + y_arr = list(np.append(np.sort(yl), np.sort(yh))) + return x_arr, y_arr + + + def get_reduced_psf(self, psf_ext, flux, cx, cy): + """ + Evaluate the appropriate (chosen based on psf_ext) GriddedPSFModel PSF, and then + bin it down to detector sampling, bearing in mind that the centre pixel always + needs to be the centre of a bin. + """ + psf_data = self.data[psf_ext] + ix, iy = int(np.floor(cx)), int(np.floor(cy)) + oversample = self.meta['oversampling'] + psf_pix = self.meta['psf_fov_pixels'] + + psf_y, psf_x = psf_data.data.shape[1], psf_data.data.shape[2] + ly, hy = iy - psf_y//2, iy + psf_y//2 + 1 + lx, hx = ix - psf_x//2, ix + psf_x//2 + 1 + y_arr,x_arr = np.mgrid[ly:hy, lx:hx] + psf_img = psf_data.evaluate(x_arr, y_arr, flux=flux, x_0=cx, y_0=cy) + + kernel = make_epsf_kernel(oversample) +# print("PSF Image Size: {}".format(psf_img.shape)) + x_list, y_list = self.make_bin(psf_x, psf_y) + binned_psf = np.zeros((len(y_list), len(x_list)), dtype=np.float64) + bx, by = 0, 0 + for y in y_list: + for x in x_list: + kyl, kyh = y - oversample//2, y + oversample//2 + 1 + kxl, kxh = x - oversample//2, x + oversample//2 + 1 +# print("(x,y) = ({},{})".format(x,y)) +# print("\tUsing img[{}:{},{}:{}]".format(kyl, kyh, kxl, kxh)) + binned_psf[by,bx] = np.sum(kernel[:,:]*psf_img[kyl:kyh,kxl:kxh]) + bx += 1 + by += 1 + bx = 0 + + return binned_psf + + @staticmethod + def pix_limit(telescope, instrument, detector): + roman_det = defaultdict(lambda: {'low': 4, 'high': 4092}) + roman_ins = defaultdict(lambda: roman_det) + miri_det = defaultdict(lambda: {'low': 4, 'high': 1024}) + nircam_det = defaultdict(lambda: {'low': 0, 'high': 2040}) + detector_dict = { + 'roman': roman_ins, + 'jwst': { + 'nircam': nircam_det, + 'miri': miri_det + } + } + return detector_dict[telescope][instrument][detector] + @staticmethod def get_instrument(telescope, instrument): diff --git a/stips/utilities/utilities.py b/stips/utilities/utilities.py index e51b973..7f61cb9 100644 --- a/stips/utilities/utilities.py +++ b/stips/utilities/utilities.py @@ -388,6 +388,6 @@ def getReducedPsf(psf_params, flux, cx, cy): psf_img = psf_data.evaluate(x_arr, y_arr, flux=flux, x_0=cx, y_0=cy) if np.sum(psf_img) == 0.: print("ERROR: ZERO FLUX PSF: {} ({},{})".format(psf_params, cx, cy)) - psf_img = psf_img.reshape(psf_y//os, os, psf_x//os, os).sum(axis=(1, 3))[1:-1,1:-1] + psf_img = psf_img.reshape(psf_y//os, os, psf_x//os, os).sum(axis=(1, 3)) return psf_img diff --git a/stips/version.py b/stips/version.py index e3cfb7d..086301b 100644 --- a/stips/version.py +++ b/stips/version.py @@ -1,2 +1,2 @@ -__version__ = '2.0.0dev7' +__version__ = '2.0.0dev8' __data__version__ = '2.0.0'