Source code for photutils.datasets.images

# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""
This module provides tools for making simulated images for documentation
examples and tests.
"""

import astropy.units as u
import numpy as np
from astropy.convolution import discretize_model
from astropy.modeling import Model
from astropy.nddata.utils import NoOverlapError
from astropy.table import Table

from photutils.utils._parameters import as_pair
from photutils.utils._progress_bars import add_progress_bar
from photutils.utils.cutouts import _overlap_slices as overlap_slices

__all__ = ['make_model_image']


[docs] def make_model_image(shape, model, params_table, *, model_shape=None, bbox_factor=None, x_name='x_0', y_name='y_0', params_map=None, discretize_method='center', discretize_oversample=10, progress_bar=False): """ Make a 2D image containing sources generated from a user-specified astropy 2D model. The model parameters for each source are taken from the input ``params_table`` table. By default, the table is searched for column names that match model parameter names and the values specified by ``x_name`` and ``y_name``. However, the user can specify a different mapping between model parameter names and column names using the ``params_map`` keyword. Parameters ---------- shape : 2-tuple of int The shape of the output image. model : 2D `astropy.modeling.Model` The 2D model to be used to render the sources. The model must be two-dimensional where it accepts 2 inputs (i.e., (x, y)) and has 1 output. The model must have parameters for the x and y positions of the sources. Typically, these parameters are named 'x_0' and 'y_0', but the parameter names can be specified using the ``x_name`` and ``y_name`` keywords. params_table : `~astropy.table.Table` A table containing the model parameters for each source. Each row of the table corresponds to a source whose model parameters are defined by the column names, which must match the model parameter names. The table must contain columns for the x and y positions of the sources. The column names for the x and y positions can be specified using the ``x_name`` and ``y_name`` keywords. Model parameters not defined in the table or ``params_maps`` will be set to the ``model`` default value. To attach units to model parameters, ``params_table`` must be input as a `~astropy.table.QTable`. If the table contains a column named 'model_shape', then the values in that column will be used to override the ``model_shape`` keyword and model ``bounding_box`` for each source. This can be used to render each source with a different shape. If the table contains a column named 'local_bkg', then the per-pixel local background values in that column will be used to added to each model source over the region defined by its ``model_shape``. The 'local_bkg' column must have the same flux units as the output image (e.g., if the input ``model`` has 'amplitude' or 'flux' parameters with units). Including 'local_bkg' should be used with care, especially in crowded fields where the ``model_shape`` of sources overlap (see Notes below). Except for ``model_shape`` and ``local_bkg`` column names, column names that do not match model parameters will be ignored unless ``params_map`` is input. model_shape : 2-tuple of int, int, or `None`, optional The shape around the (x, y) center of each source that will used to evaluate the ``model``. If ``model_shape`` is a scalar integer, then a square shape of size ``model_shape`` will be used. If `None`, then the bounding box of the model will be used (which can optionally be scaled using the ``bbox_factor`` keyword if the model supports it). This keyword must be specified if the model does not have a ``bounding_box`` attribute. If specified, this keyword overrides the model ``bounding_box`` attribute. To use a different shape for each source, include a column named ``'model_shape'`` in the ``params_table``. For that case, this keyword is ignored. bbox_factor : `None` or float, optional The multiplicative factor to pass to the model ``bounding_box`` method to determine the model shape. If the model ``bounding_box`` method does not accept a ``factor`` keyword, then this keyword is ignored. If `None`, the default model bounding box will be used. This keyword is ignored if ``model_shape`` is specified or if the ``params_table`` contains a ``'model_shape'`` column. Note that some Photutils PSF models have a ``bbox_factor`` keyword that is be used to define the model bounding box. In that case, this keyword is ignored. x_name : str, optional The name of the ``model`` parameter that corresponds to the x position of the sources. If ``param_map`` is not input, then this value must also be a column name in ``params_table``. y_name : str, optional The name of the ``model`` parameter that corresponds to the y position of the sources. If ``param_map`` is not input, then this value must also be a column name in ``params_table``. params_map : dict or None, optional A dictionary mapping the model parameter names to the column names in the input ``params_table``. The dictionary keys are the model parameter names and the values are the column names in the input ``params_table``. This can be used to map column names to model parameter names that are different. For example, if the input column name is 'flux_f200w' and the model parameter name is 'flux', then use ``column_map={'flux': 'flux_f200w'}``. This table may also be used if you want to map the model x and y parameters to different columns than ``x_name`` and ``y_name``, but the ``x_name`` and ``y_name`` keys must be included in the dictionary. discretize_method : {'center', 'interp', 'oversample', 'integrate'}, \ optional One of the following methods for discretizing the model on the pixel grid: * ``'center'`` (default) Discretize model by taking the value at the center of the pixel bins. This method should be used for ePSF/PRF single or gridded models. * ``'interp'`` Discretize model by bilinearly interpolating between the values at the corners of the pixel bins. * ``'oversample'`` Discretize model by taking the average of model values in the pixel bins on an oversampled grid. Use the ``discretize_oversample`` keyword to set the integer oversampling factor. * ``'integrate'`` Discretize model by integrating the model over the pixel bins using `scipy.integrate.quad`. This mode conserves the model integral on a subpixel scale, but it is *extremely* slow. discretize_oversample : int, optional The integer oversampling factor used when ``descretize_method='oversample'``. This keyword is ignored otherwise. progress_bar : bool, optional Whether to display a progress bar while adding the sources to the image. The progress bar requires that the `tqdm <https://tqdm.github.io/>`_ optional dependency be installed. Note that the progress bar does not currently work in the Jupyter console due to limitations in ``tqdm``. Returns ------- array : 2D `~numpy.ndarray` The rendered image containing the model sources. Examples -------- .. plot:: :include-source: import matplotlib.pyplot as plt from astropy.modeling.models import Moffat2D from photutils.datasets import (make_model_image, make_random_models_table) model = Moffat2D() n_sources = 25 shape = (100, 100) param_ranges = {'amplitude': [100, 200], 'x_0': [0, shape[1]], 'y_0': [0, shape[0]], 'gamma': [1, 2], 'alpha': [1, 2]} params = make_random_models_table(n_sources, param_ranges, seed=0) model_shape = (15, 15) data = make_model_image(shape, model, params, model_shape=model_shape) plt.imshow(data, origin='lower') plt.tight_layout() .. plot:: :include-source: import matplotlib.pyplot as plt import numpy as np from astropy.modeling.models import Gaussian2D from photutils.datasets import make_model_image, make_model_params model = Gaussian2D() shape = (500, 500) n_sources = 100 params = make_model_params(shape, n_sources, x_name='x_mean', y_name='y_mean', min_separation=25, amplitude=(100, 500), x_stddev=(1, 3), y_stddev=(1, 3), theta=(0, np.pi)) model_shape = (25, 25) data = make_model_image(shape, model, params, model_shape=model_shape, x_name='x_mean', y_name='y_mean') plt.imshow(data, origin='lower') plt.tight_layout() Notes ----- The local background value around each source is optionally included using the ``local_bkg`` column in the input ``params_table``. This local background added to each source over its ``model_shape`` region. In regions where the ``model_shape`` of source overlap, the local background will be added multiple times. This is not an issue if the sources are well-separated, but for crowded fields, this option should be used with care. """ if not isinstance(shape, tuple) or len(shape) != 2: raise ValueError('shape must be a 2-tuple') if not isinstance(model, Model): raise TypeError('model must be a Model instance') if model.n_inputs != 2 or model.n_outputs != 1: raise ValueError('model must be a 2D model') if not isinstance(params_table, Table): raise TypeError('params_table must be an astropy Table') xypos_map = {x_name: x_name, y_name: y_name} # by default, use the model parameter names as the column names # if they are in the table params_to_set = set(params_table.colnames) & set(model.param_names) xypos_map.update({param: param for param in params_to_set}) if params_map is not None: # params_map takes precedence over x_name and y_name and # any matching column names in params_table xypos_map.update(params_map) params_map = xypos_map for key, value in params_map.items(): if key not in model.param_names: raise ValueError(f'key "{key}" not in model parameter names') if value not in params_table.colnames: raise ValueError(f'value "{value}" not in params_table column ' 'names') if model_shape is not None: model_shape = as_pair('model_shape', model_shape, lower_bound=(0, 1)) variable_shape = False if 'model_shape' in params_table.colnames: model_shape = np.array(params_table['model_shape']) if model_shape.ndim == 1: model_shape = np.array([as_pair('model_shape', shape) for shape in model_shape]) variable_shape = True if model_shape is None: try: _ = model.bounding_box except NotImplementedError as exc: raise ValueError('model_shape must be specified if the model ' 'does not have a bounding_box attribute') from exc if 'local_bkg' in params_table.colnames: local_bkg = params_table['local_bkg'] else: local_bkg = np.zeros(len(params_table)) # copy the input model to leave it unchanged model = model.copy() if progress_bar: # pragma: no cover desc = 'Add model sources' params_table = add_progress_bar(params_table, desc=desc) image = np.zeros(shape, dtype=float) for i, source in enumerate(params_table): for key, param in params_map.items(): setattr(model, key, source[param]) # This assumes that if the user also uses params_table to # override the (x/y)_name mapping that the x_name and y_name # values are correct (i.e., the mapping keys include x_name and # y_name). There is no good way to check/enforce this. x0 = getattr(model, x_name).value y0 = getattr(model, y_name).value if variable_shape: mod_shape = model_shape[i] elif model_shape is None: # the bounding box size generally depends on model parameters, # so needs to be calculated for each source mod_shape = _model_shape_from_bbox(model, bbox_factor=bbox_factor) else: mod_shape = model_shape try: slc_lg, _ = overlap_slices(shape, mod_shape, (y0, x0), mode='trim') if discretize_method == 'center': yy, xx = np.mgrid[slc_lg] subimg = model(xx, yy) else: if discretize_method == 'interp': discretize_method = 'linear_interp' x_range = (slc_lg[1].start, slc_lg[1].stop) y_range = (slc_lg[0].start, slc_lg[0].stop) subimg = discretize_model(model, x_range=x_range, y_range=y_range, mode=discretize_method, factor=discretize_oversample) if i == 0 and isinstance(subimg, u.Quantity): image <<= subimg.unit try: image[slc_lg] += subimg + local_bkg[i] except u.UnitConversionError as exc: raise ValueError('The local_bkg column must have the same ' 'flux units as the output image') from exc except NoOverlapError: continue return image
def _model_shape_from_bbox(model, bbox_factor=None): """ Calculate the model shape from the model bounding box. Parameters ---------- model : 2D `astropy.modeling.Model` The 2D model to be used to render the sources. bbox_factor : `None` or float, optional The multiplicative factor to pass to the model ``bounding_box`` method to determine the model shape. If the model ``bounding_box`` method does not accept a ``factor`` keyword, then this keyword is ignored. If `None`, the default model bounding box will be used. Returns ------- model_shape : 2-tuple of int The shape around the (x, y) center of the model that will used to evaluate the model. Raises ------ ValueError If the model does not have a bounding_box attribute. """ try: hasattr(model, 'bounding_box') except NotImplementedError as exc: msg = 'model does not have a bounding_box attribute' raise ValueError(msg) from exc if bbox_factor is not None: try: bbox = model.bounding_box(factor=bbox_factor) except NotImplementedError: bbox = model.bounding_box.bounding_box() else: bbox = model.bounding_box.bounding_box() return (int(np.ceil(bbox[0][1] - bbox[0][0])), int(np.ceil(bbox[1][1] - bbox[1][0])))