Photometric calibration

This page covers photometric calibration — matching detected objects to reference catalogues and fitting zero-point models. For object detection see Object detection and measurement, for flux measurement methods (aperture, optimal extraction, PSF photometry) see the measurement sections in Object detection and measurement and Point Spread Function (PSF) models, and for detection flags see Object flags.

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 stdpipe.photometry.match() function, and wrapped into a bit higher-level and easier to use 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:

\[{\rm Catalogue} = {\rm Instrumental} + C\cdot{\rm color} + ZP(x, y, {\rm spatial\_order}) + {\rm additive} + {\rm error}\]

where

  • \({\rm Instrumental}=-2.5\log_{10}({\rm ADU})\) is instrumental magnitude,

  • \({\rm Catalogue}\) is catalogue magnitude, e.g. \({\rm V}\), and

  • \({\rm color}\) is catalogue color, e.g. \({\rm B-V}\).

This is equivalent to the instrumenmtal photometric system defined as \({\rm V} - C\cdot({\rm B-V})\).

Zero point \(ZP\) is a spatially varying polynomial with the degree controlled by spatial_order parameter (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

\[{\rm additive} = -2.5/\log(10)/10^{-0.4\cdot{\rm Instrumental}} \cdot {\rm bg_{corr}}(x, y, {\rm bg\_order})\]

where \({\rm bg_{corr}}\) is a flux correction inside the aperture. This term is also spatially dependent, and is controlled by bg_order parameter. If bg_order=None, the fitting for this term is disabled.

The calibration routine performs an iterative weighted linear least square (or robust if 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.

# 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']
stdpipe.photometry.match(obj_ra, obj_dec, obj_mag, obj_magerr, obj_flags, cat_ra, cat_dec, cat_mag, cat_magerr=None, cat_color=None, sr=0.0008333333333333334, obj_x=None, obj_y=None, spatial_order=0, bg_order=None, nonlin=False, threshold=5.0, niter=10, accept_flags=0, cat_saturation=None, max_intrinsic_rms=0, sn=None, verbose=False, robust=True, scale_noise=False, use_color=True, force_color_term=None)[source]

Low-level photometric matching routine.

Builds a 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.

Parameters:
obj_raarray-like

Right Ascension values for the objects.

obj_decarray-like

Declination values for the objects.

obj_magarray-like

Instrumental magnitude values for the objects.

obj_magerrarray-like

Instrumental magnitude errors for the objects.

obj_flagsarray-like

Flags for the objects.

cat_raarray-like

Catalogue Right Ascension values.

cat_decarray-like

Catalogue Declination values.

cat_magarray-like

Catalogue magnitudes.

cat_magerrarray-like, optional

Catalogue magnitude errors.

cat_colorarray-like, optional

Catalogue color values.

srfloat

Matching radius in degrees.

obj_xarray-like, optional

X coordinates of objects on the image.

obj_yarray-like, optional

Y coordinates of objects on the image.

spatial_orderint

Order of zero point spatial polynomial (0 for constant).

bg_orderint or None

Order of additive flux term spatial polynomial. None to disable this term in the model.

nonlinbool

Whether to fit for simple non-linearity.

thresholdfloat

Rejection threshold (relative to magnitude errors) for object-catalogue pair to be rejected from the fit.

niterint

Number of iterations for the fitting.

accept_flagsint

Bitmask for acceptable object flags. Objects having any other bits set will be excluded from the model.

cat_saturationfloat or None

Saturation level for the catalogue. Stars brighter than this magnitude will be excluded from the fit.

max_intrinsic_rmsfloat

Maximal intrinsic RMS to use during the fitting. If set to 0, no intrinsic scatter is included in the noise model.

snfloat or None

Minimal acceptable signal to noise ratio (1/obj_magerr) for the objects to be included in the fit.

verbosebool or callable

Whether to show verbose messages during the run. May be either boolean, or a print-like function.

robustbool

Whether to use robust least squares fitting instead of weighted least squares.

scale_noisebool

Whether to re-scale the noise model (object and catalogue magnitude errors) to match actual scatter of the data points. Intrinsic scatter term is not scaled this way.

use_colorbool or int

Whether to use catalogue color for deriving the color term. If integer, it determines the color term order.

force_color_termfloat or None

Do not fit for the color term, but use this fixed value instead.

Returns:
dict or None

Dictionary with photometric results, or None if the fit failed. Contains the following keys:

  • oidx, cidx, dist – indices of positionally matched objects and catalogue stars, and their pairwise distances in degrees.

  • omag, omag_err, cmag, cmag_err – instrumental magnitudes of matched objects, corresponding catalogue magnitudes, and their errors. Array lengths equal 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 minus instrumental magnitudes) for every matched object, and 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 errors from the fit.

  • color_term – fitted color term. Instrumental photometric system is defined as 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, and optionally its error.

  • obj_zero – zero points for all input objects (not necessarily matched to the catalogue) computed via zero_fn, without color term.

  • params – internal parameters of the fitting 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 (not rejected during iterative thresholding, and passing initial quality cuts).

  • idx0 – same as idx but with only initial quality cuts applied.

Notes

The returned zero point function zero_fn has the signature:

zero_fn(xx, yy, mag=None, get_err=False, add_intrinsic_rms=False)

where xx and yy are coordinates on the image, mag is the object instrumental magnitude (needed to compute the additive flux term). If 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 the intrinsic scatter term.

The zero point returned by this function does not include the contribution of the color term. To derive the final calibrated magnitude for an object, manually add the color contribution:

mag_calibrated = mag_instrumental + color * color_term

where color is the true object color, and color_term is reported in the photometric results.

stdpipe.pipeline.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_flags='flags', 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)[source]

Higher-level photometric calibration routine.

Wraps stdpipe.photometry.match() with convenient defaults for typical tabular data.

Parameters:
objastropy.table.Table

Table of detected objects.

catastropy.table.Table

Reference photometric catalogue.

srfloat, optional

Matching radius in degrees. Defaults to half the median FWHM in sky units if pixscale is provided, otherwise 1 arcsec.

pixscalefloat, optional

Pixel scale in degrees/pixel; used to compute the default matching radius.

orderint, optional

Spatial polynomial order for the zero-point model (0 = constant).

bg_orderint or None, optional

Spatial polynomial order for an additive flux background term. None disables this term.

obj_col_magstr, optional

Column name for object instrumental magnitude.

obj_col_mag_errstr, optional

Column name for object magnitude error.

obj_col_rastr, optional

Column name for object Right Ascension.

obj_col_decstr, optional

Column name for object Declination.

obj_col_flagsstr, optional

Column name for object flags.

obj_col_xstr, optional

Column name for object x coordinate.

obj_col_ystr, optional

Column name for object y coordinate.

cat_col_magstr, optional

Column name for catalogue magnitude.

cat_col_mag_errstr, optional

Column name for catalogue magnitude error.

cat_col_mag1str, optional

First catalogue magnitude column for the color term (e.g. 'B').

cat_col_mag2str, optional

Second catalogue magnitude column for the color term (e.g. 'R').

cat_col_rastr, optional

Column name for catalogue Right Ascension.

cat_col_decstr, optional

Column name for catalogue Declination.

updatebool, optional

If True, add mag_calib and mag_calib_err columns to obj in-place with the calibrated magnitude and its error.

verbosebool or callable, optional

Whether to show verbose messages. May be boolean or a print-like callable.

**kwargs

Additional keyword arguments passed to stdpipe.photometry.match().

Returns:
dict or None

Dictionary with photometric results as returned by stdpipe.photometry.match(), or None on failure.

Plotting photometric match results

STDPipe also includes a dedicated plotting routine stdpipe.plots.plot_photometric_match() that may be used to visually inspect various aspects of the photometric match.

# 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()
stdpipe.plots.plot_photometric_match(m, ax=None, mode='mag', show_masked=True, show_final=True, **kwargs)[source]

Plot photometric match diagnostics.

Displays various representations of the results returned by stdpipe.photometry.match() or stdpipe.pipeline.calibrate_photometry().

Available modes:

  • 'mag' — photometric residuals vs catalogue magnitude

  • 'normed' — residuals divided by errors vs catalogue magnitude

  • 'color' — residuals vs catalogue color

  • 'zero' — spatial map of the empirical zero point

  • 'model' — spatial map of the zero-point model

  • 'residuals' — spatial map of zero-point fitting residuals

  • 'dist' — spatial map of angular separation (arcsec)

Parameters:
mdict

Photometric match results dictionary.

axmatplotlib.axes.Axes, optional

Axes to draw on; defaults to the current axes.

modestr, optional

Plotting mode (see above).

show_maskedbool, optional

If True, include masked objects in the plot.

show_finalbool, optional

If True, highlight objects used in the final fit.

**kwargs

Additional keyword arguments passed to binned_map() where applicable.

Detection limits and S/N modeling

STDPipe can estimate the detection limit magnitude by fitting a model of S/N versus magnitude and finding where it crosses a given threshold. The higher-level stdpipe.pipeline.get_detection_limit() works directly on calibrated object tables, while the lower-level functions give more control.

# Higher-level: estimate 5-sigma detection limit from calibrated objects
mag_lim = pipeline.get_detection_limit(obj, sn=5, verbose=True)
print('Detection limit: %.2f mag' % mag_lim)

# Lower-level: build S/N model and find limit manually
sn_model = photometry.make_sn_model(obj['mag_calib'], 1/obj['magerr'])
mag_lim, sn_model = photometry.get_detection_limit_sn(
    obj['mag_calib'], 1/obj['magerr'], sn=5, get_model=True)

# Plot S/N model
mags = np.linspace(np.min(obj['mag_calib']), mag_lim + 1, 100)
plt.scatter(obj['mag_calib'], 1/obj['magerr'], s=1, alpha=0.3)
plt.plot(mags, sn_model(mags), 'r-')
plt.axhline(5, ls='--', color='k', label='S/N = 5')
plt.axvline(mag_lim, ls='--', color='g', label='Detection limit')
plt.xlabel('Calibrated magnitude')
plt.ylabel('S/N')
plt.legend()
plt.show()
stdpipe.pipeline.get_detection_limit(obj, sn=5, method='sn', verbose=True)[source]

Estimate the detection limit magnitude.

The object table must contain calibrated magnitudes and errors in mag_calib and mag_calib_err columns.

Parameters:
objastropy.table.Table

Table with calibrated objects.

snfloat, optional

S/N threshold defining the detection limit.

methodstr, optional

Estimation method. Currently only 'sn' (S/N vs magnitude extrapolation) is implemented.

verbosebool or callable, optional

Whether to show verbose messages. May be boolean or a print-like callable.

Returns:
float or None

Magnitude corresponding to the detection limit at the given S/N level, or None if estimation fails.

stdpipe.photometry.make_sn_model(mag, sn)[source]

Build a model for signal to noise (S/N) ratio versus magnitude.

Assumes the noise comes from constant background noise plus Poissonian noise with constant gain.

Parameters:
magarray-like

Calibrated magnitudes.

snarray-like

S/N values corresponding to the magnitudes.

Returns:
callable

Function that accepts an array of magnitudes and returns the S/N model values for them.

Raises:
ValueError

If no finite positive S/N values are available to build the model.

stdpipe.photometry.get_detection_limit_sn(mag, mag_sn, sn=5, get_model=False, verbose=True)[source]

Estimate the detection limit using S/N vs magnitude method.

Parameters:
magarray-like

Calibrated magnitudes.

mag_snarray-like

S/N values corresponding to these magnitudes.

snfloat

S/N level for the detection limit.

get_modelbool

If True, also return the S/N model function.

verbosebool or callable

Whether to show verbose messages during the run. May be either boolean, or a print-like function.

Returns:
float or None

Magnitude corresponding to the detection limit at the given S/N level. None if the model cannot be built or the root is not found.

callable or None

S/N model function (only returned when get_model=True). None if the model cannot be built.