Commit 5dfe9db0 authored by karpov-sv's avatar karpov-sv
Browse files

Photometric calibration documented

parent 424599d4
Photometric calibration
=======================
Photometric calibration is performed by positionally matching detected objects with catalogue stars, and then building the photometric model for their instrumental magnitudes.
The model includes:
- catalogue magnitude
- (optionally spatially varying) zero point
- (optionally) catalogue color
- (very optionally) additive flux term, corresponding e.g. to biased background estimation
- (very very optionally) intrinsic scatter on top of measurement errors
and its fitting is implemented in :func:`stdpipe.photometry.match` function, and wrapped into a bit higher-level and easier to use :func:`stdpipe.pipeline.calibrate_photometry`. The routines perform an iterative fitting with rejection of pairs deviating too much (more than `threshold` sigmas) from the model that is defined like that:
.. math::
{\rm Catalogue} = {\rm Instrumental} + C\cdot{\rm color} + ZP(x, y, {\rm spatial_order}) + {\rm additive} + {\rm error}
where
* :math:`{\rm Instrumental}=-2.5\log_{10}({\rm ADU})` is instrumental magnitude,
* :math:`{\rm Catalogue}` is catalogue magnitude, e.g. :math:`{\rm V}`, and
* :math:`{\rm color}` is catalogue color, e.g. :math:`{\rm B-V}`.
This is equivalent to the instrumenmtal photometric system defined as :math:`{\rm V} - C\cdot({\rm B-V})`.
Zero point :math:`ZP` is a spatially varying polynomial with the degree controlled by `spatial_order` parameter (:code:`spatial_order=0` for a constant zero point)
Additive flux term is defined by linearizing the additional flux in every photometric aperture (e.g. due to incorrect background level determination) and has a form
.. math::
{\rm additive} = -2.5/\log(10)/10^{-0.4\cdot{\rm Instrumental}} \cdot {\rm bg_corr}(x, y, {\rm bg_order})
where :math:`{\rm bg_corr}` is a flux correction inside the aperture. This term is also spatially dependent, and is controlled by `bg_order` parameter. If :code:`bg_order=None`, the fitting for this term is disabled.
The calibration routine performs an iterative weighted linear least square (or robust if :code:`robust=True`) fitting with rejection of pairs deviating too much (more than `threshold` sigmas) from the model. Optional intrinsic scatter (specified through `max_intrinsic_rms` parameter) may also be fitted for, and may help accounting for the effects of e.g. multiplicative noise (flatfielding, subpixel sensitivity variations, etc).
.. attention::
The calibration (or zero point, or systematic) error is currently estimated as a fitted model error, and is typically way too small. Most probably the method is wrong! So I welcome any input or discussions on this topic.
.. code-block:: python
# Photometric calibration using 2 arcsec matching radius, r magnitude, g-r color and second order spatial variations
m = pipeline.calibrate_photometry(obj, cat, sr=2/3600, cat_col_mag='rmag', cat_col_mag1='gmag', cat_col_mag2='rmag', max_intrinsic_rms=0.02, order=2, verbose=True)
# The code above automatically augments the object list with calibrated magnitudes, but we may also do it manually
obj['mag_calib'] = obj['mag'] + m['zero_fn'](obj['x'], obj['y'])
obj['mag_calib_err'] = np.hypot(obj['magerr'], m['zero_fn'](obj['x'], obj['y'], get_err=True))
# Now, if we happen to know object colors, we may also compute proper color-corrected magnitudes:
obj['mag_calib_color'] = obj['mag_calib'] + obj['color']*m['color_term']
.. autofunction:: stdpipe.photometry.match
:noindex:
.. autofunction:: stdpipe.pipeline.calibrate_photometry
:noindex:
Plotting photometric match results
----------------------------------
*STDPipe* also includes a dedicated plotting routine :func:`stdpipe.plots.plot_photometric_match` that may be used to visually inspect various aspects of the photometric match.
.. code-block:: python
# Photometric residuals as a function of catalogue magnitude
plt.subplot(211)
plots.plot_photometric_match(m)
plt.ylim(-0.5, 0.5)
# Photometric residuals as a function of catalogue color
plt.subplot(212)
plots.plot_photometric_match(m, mode='color')
plt.ylim(-0.5, 0.5)
plt.xlim(0.0, 1.5)
plt.show()
# Zero point (difference between catalogue and instrumental magnitudes for every star) map
plots.plot_photometric_match(m, mode='zero', bins=6, show_dots=True, aspect='equal')
plt.show()
# Fitted zero point model map
plots.plot_photometric_match(m, mode='model', bins=6, show_dots=True, aspect='equal')
plt.show()
# Residuals between empirical zero point and fitted model
plots.plot_photometric_match(m, mode='residuals', bins=6, show_dots=True, aspect='equal')
plt.show()
# Astrometric displacement between objects and matched catalogue stars
plots.plot_photometric_match(m, mode='dist', bins=6, show_dots=True, aspect='equal')
plt.show()
.. autofunction:: stdpipe.plots.plot_photometric_match
:noindex:
......@@ -481,6 +481,62 @@ def get_intrinsic_scatter(y, yerr, min=0, max=None):
return C.x[2]
def match(obj_ra, obj_dec, obj_mag, obj_magerr, obj_flags, cat_ra, cat_dec, cat_mag, cat_magerr=None, cat_color=None, sr=3/3600, obj_x=None, obj_y=None, spatial_order=0, bg_order=None, threshold=5.0, niter=10, accept_flags=0, cat_saturation=None, max_intrinsic_rms=0, sn=None, verbose=False, robust=True, scale_noise=False):
"""Low-level photometric matching routine.
It tries to build the photometric model for objects detected on the image that includes catalogue magnitude, positionally-dependent zero point, linear color term, optional additive flux term, and also takes into account possible intrinsic magnitude scatter on top of measurement errors.
:param obj_ra: Array of Right Ascension values for the objects
:param obj_dec: Array of Declination values for the objects
:param obj_mag: Array of instrumental magnitude values for the objects
:param obj_magerr: Array of instrumental magnitude errors for the objects
:param obj_flags: Array of flags for the objects
:param cat_ra: Array of catalogue Right Ascension values
:param cat_dec: Array of catalogue Declination values
:param cat_mag: Array of catalogue magnitudes
:param cat_magerr: Array of catalogue magnitude errors
:param cat_color: Array of catalogue color values, optional
:param sr: Matching radius, degrees
:param obj_x: Array of `x` coordinates of objects on the image, optional
:param obj_y: Array of `y` coordinates of objects on the image, optional
:param spatial_order: Order of zero point spatial polynomial (0 for constant).
:param bg_order: Order of additive flux term spatial polynomial (None to disable this term in the model)
:param threshold: Rejection threshold (relative to magnitude errors) for object-catalogue pair to be rejected from the fit
:param niter: Number of iterations for the fitting
:param accept_flags: Bitmask for acceptable object flags. Objects having any other
:param cat_saturation: Saturation level for the catalogue - stars brighter than this magnitude will be excluded from the fit
:param max_intrinsic_rms: Maximal intrinsic RMS to use during the fitting. If set to 0, no intrinsic scatter is included in the noise model.
:param sn: Minimal acceptable signal to noise ratio (1/obj_magerr) for the objects to be included in the fit
:param verbose: Whether to show verbose messages during the run of the function or not. May be either boolean, or a `print`-like function.
:param robust: Whether to use robust least squares fitting routine instead of weighted least squares
:param scale_noise: Whether to re-scale the noise model (object and catalogue magnitude errors) to match actual scatter of the data points or not. Intrinsic scatter term is not being scaled this way.
:returns: The dictionary with photometric results, as described below.
The results of photometric matching are returned in a dictionary with the following fields:
- `oidx`, `cidx`, `dist` - indices of positionally matched objects and catalogue stars, as well as their pairwise distances in degrees
- `omag`, `omagerr`, `cmag`, `cmagerr` - arrays of object instrumental magnitudes of matched objects, corresponding catalogue magnitudes, and their errors. Array lengths are equal to the number of positional matches.
- `color` - catalogue colors corresponding to the matches, or zeros if no color term fitting is requested
- `ox`, `oy`, `oflags` - coordinates of matched objects on the image, and their flags
- `zero`, `zero_err` - empirical zero points (catalogue - instrumental magnitudes) for every matched object, as well as its errors, derived as a hypotenuse of their corresponding errors.
- `zero_model`, `zero_model_err` - modeled "full" zero points (including color terms) for matched objects, and their corresponding errors from the fit
- `color_term` - fitted color term. Instrumental photometric system is defined as :code:`obj_mag = cat_mag - color*color_term`
- `zero_fn` - function to compute the zero point (without color term) at a given position and for a given instrumental magnitude of object, and optionally its error.
- `obj_zero` - zero points for all input objects (not necessarily matched to the catalogue) computed through aforementioned function, i.e. without color term
- `params` - Internal parameters of the fittting polynomial
- `intrinsic_rms`, `error_scale` - fitted values of intrinsic scatter and noise scaling
- `idx` - boolean index of matched objects/catalogue stars used in the final fit (i.e. not rejected during iterative thresholding, and passing initial quality cuts
- `idx0` - the same but with just initial quality cuts taken into account
Returned zero point computation function has the following signature:
:obj:`zero_fn(xx, yy, mag=None, get_err=False, add_intrinsic_rms=False)`
where `xx` and `yy` are coordinates on the image, `mag` is object instrumental magnitude (needed to compute additive flux term). If :code:`get_err=True`, the function returns estimated zero point error instead of zero point, and `add_intrinsic_rms` controls whether this error estimation should also include intrinsic scatter term or not.
The zero point returned by this function does not include the contribution of color term. Therefore, in order to derive the final calibrated magnitude for the object, you will need to manually add the color contribution: :code:`mag_calibrated = mag_instrumental + color*color_term`, where `color` is a true object color, and `color_term` is reported in the photometric results.
"""
# Simple wrapper around print for logging in verbose mode only
log = (verbose if callable(verbose) else print) if verbose else lambda *args,**kwargs: None
......@@ -625,7 +681,7 @@ def match(obj_ra, obj_dec, obj_mag, obj_magerr, obj_flags, cat_ra, cat_dec, cat_
'params': C.params,
'error_scale': np.sqrt(C.scale),
'intrinsic_rms': intrinsic_rms,
'obj_zero': zero_fn(obj_x, obj_y),
'obj_zero': zero_fn(obj_x, obj_y, mag=obj_mag),
'ox': ox, 'oy': oy, 'oflags': oflags,
'idx': idx, 'idx0': idx0}
......
......@@ -241,14 +241,41 @@ def filter_transient_candidates(obj, sr=None, pixscale=None, time=None,
def calibrate_photometry(obj, cat, sr=None, pixscale=None, order=0, bg_order=None,
obj_col_mag='mag', obj_col_mag_err='magerr',
obj_col_ra='ra', obj_col_dec='dec',
obj_col_x='x', obj_col_y='y',
cat_col_mag='R', cat_col_mag_err=None,
cat_col_mag1=None, cat_col_mag2=None,
cat_col_ra='RAJ2000', cat_col_dec='DEJ2000',
update=True, verbose=False, **kwargs):
"""
Higher-level photometric calibration routine
"""
"""Higher-level photometric calibration routine.
It wraps :func:`stdpipe.photometry.match` routine with some convenient defaults so that it is easier to use with typical tabular data.
:param obj: Table of detected objects
:param cat: Reference photometric catalogue
:param sr: Matching radius in degrees, optional
:param pixscale: Pixel scale, degrees per pixel. If specified, and `sr` is not set, then median value of half of FWHM, multiplied by pixel scale, is used as a matching radius.
:param order: Order of zero point spatial polynomial (0 for constant).
:param bg_order: Order of additive flux term spatial polynomial (None to disable this term in the model)
:param obj_col_mag: Column name for object instrumental magnitude
:param obj_col_mag_err: Column name for object magnitude error
:param obj_col_ra: Column name for object Right Ascension
:param obj_col_dec: Column name for object Declination
:param obj_col_x: Column name for object x coordinate
:param obj_col_y: Column name for object y coordinate
:param cat_col_mag: Column name for catalogue magnitude
:param cat_col_mag_err: Column name for catalogue magnitude error
:param cat_col_mag1: Column name for the first catalogue magnitude defining the stellar color
:param cat_col_mag2: Column name for the second catalogue magnitude defining the stellar color
:param cat_col_ra: Column name for catalogue Right Ascension
:param cat_col_dec: Column name for catalogue Declination
:param update: If True, `mag_calib` and `mag_calib_err` columns with calibrated magnitude (without color term) and its error will be added to the object table
:param verbose: Whether to show verbose messages during the run of the function or not. May be either boolean, or a `print`-like function
:param kwargs: The rest of keyword arguments will be directly passed to :func:`stdpipe.photometry.match`.
:returns: The dictionary with photometric results, as returned by :func:`stdpipe.photometry.match`.
"""
# Simple wrapper around print for logging in verbose mode only
log = (verbose if callable(verbose) else print) if verbose else lambda *args,**kwargs: None
......
......@@ -30,7 +30,17 @@ def colorbar(obj=None, ax=None, size="5%", pad=0.1):
ax.get_figure().sca(ax)
def imshow(image, qq=None, show_colorbar=True, show_axis=True, stretch='linear', ax=None, **kwargs):
"""Simple wrapper around pyplot.imshow with histogram-based intensity scaling"""
"""Simple wrapper around pyplot.imshow with histogram-based intensity scaling, optional colorbar, etc.
:param image: Numpy 2d array to display
:param qq: two-element tuple (or list) with quantiles that define lower and upper limits for image intensity normalization. Default is `[0.5, 99.5]`. Will be superseded by manually provided `vmin` and `vmax` arguments.
:param show_colorbar: Whether to show a colorbar alongside the image
:param show_axis: Whether to show the axes around the image
:param stretch: Image intensity stretching mode - e.g. `linear`, `log`, `asinh`, or anything else supported by Astropy visualization layer
:param ax: Matplotlib Axes object to be used for plotting, optional
:param kwargs: The rest of parameters will be directly passed to :func:`matplotlib.pyplot.imshow`
"""
if ax is None:
ax = plt.gca()
......@@ -73,6 +83,22 @@ def imshow(image, qq=None, show_colorbar=True, show_axis=True, stretch='linear',
return img
def binned_map(x, y, value, bins=16, statistic='mean', qq=[0.5, 97.5], show_colorbar=True, show_axis=True, show_dots=False, ax=None, **kwargs):
"""Plots various statistical estimators binned onto regular grid from the set of irregular data points (`x`, `y`, `value`).
:param x: Abscissae of the data points
:param y: Ordinates of the data points
:param value: Values of the data points
:param bins: Number of bins per axis
:param statistic: Statistical estimator to plot, may be `mean`, `median`, or a function
:param qq: two-element tuple (or list) with quantiles that define lower and upper limits for image intensity normalization. Default is `[0.5, 97.5]`. Will be superseded by manually provided `vmin` and `vmax` arguments.
:param show_colorbar: Whether to show a colorbar alongside the image
:param show_axis: Whether to show the axes around the image
:param show_dots: Whether to overlay the positions of data points onto the plot
:param ax: Matplotlib Axes object to be used for plotting, optional
:param kwargs: The rest of parameters will be directly passed to :func:`matplotlib.pyplot.imshow`
:returns: None
"""
gmag0, xe, ye, binnumbers = binned_statistic_2d(x, y, value, bins=bins, statistic=statistic)
vmin1,vmax1 = np.percentile(gmag0[np.isfinite(gmag0)], qq)
......@@ -153,6 +179,28 @@ def plot_cutout(cutout, planes=['image', 'template', 'diff', 'mask'], fig=None,
fig.suptitle(title)
def plot_photometric_match(m, ax=None, mode='mag', show_masked=True, show_final=True, **kwargs):
"""Convenience plotting routine for photometric match results.
It plots various representations of the photometric match results returned by :func:`stdpipe.photometry.match` or :func:`stdpipe.pipeline.calibrate_photometry`, depending on the `mode` parameter:
- `mag` - displays photometric residuals as a function of catalogue magnitude
- `color` - displays photometric residuals as a function of catalogue color
- `zero` - displays the map of empirical zero point, i.e. difference of catalogue and instrumental magnitudes for all matched objects
- `model` - displays the map of zero point model
- `residuals` - displays fitting residuals between zero point and its model
- `dist` - displays the map of angular separation between matched objects and stars, in arcseconds
The parameter `show_dots` controls whether to overlay the positions of the matched objects onto the maps, when applicable.
:param m: Dictionary with photometric match results
:param ax: Matplotlib Axes object to be used for plotting, optional
:param mode: plotting mode - one of `mag`, `color`, `zero`, `model`, `residuals`, or `dist`
:param show_masked: Whether to show masked objects
:param show_final: Whether to additionally highlight the objects used for the final fit, i.e. not rejected during iterative thresholding
:param kwargs: the rest of parameters will be directly passed to :func:`stdpipe.plots.binned_map` when applicable.
:returns: None
"""
if ax is None:
ax = plt.gca()
......@@ -224,10 +272,24 @@ from contextlib import contextmanager
@contextmanager
def figure_saver(filename=None, show=False, tight_layout=True, **kwargs):
'''
Simple matplotlib Figure() wrapper, implemented as a context manager.
"""Simple matplotlib Figure() wrapper, implemented as a context manager.
It stores the figure to specified file, and optionally displays it interactively if run inside Jupyter.
'''
Intended to be used as:
.. code-block:: python
with figure_saver('/tmp/figure.png', show=True, figsize=(10, 6)) as fig:
ax = fig.add_subplot(111)
ax.plot(x, y, '.-')
:param filename: Name of a file where to store the image. May be in any format supported by Matplotlib
:param show: Whether to also display the figure inside Jupuyter notebook
:param tight_layout: Whether to call :code:`fig.tight_layout()` on the figure before saving/displaying it
:param kwargs: The rest of parameters will be directly passed to :func:`matplotlib.pyplot.Figure`
"""
fig = plt.Figure(**kwargs)
try:
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment