207.1. Timeseries values#
207.1. Timeseries values¶
For the Rubin Science Platform at data.lsst.cloud.
Data Release: Data Preview 2
Container Size: Large
LSST Science Pipelines version: r30.0.10
Last verified to run: 2026-07-24
Repository: github.com/lsst/tutorial-notebooks
DOI: 10.11578/rubin/dc.20250909.20
Learning objective: Understand how the values of the timeseries features in the DiaObject table are derived.
LSST data products: DiaObject, DiaSource
Packages: scipy
Credit: Originally developed by the Rubin Community Science team. Please consider acknowledging them if this notebook is used for the preparation of journal articles, software releases, or other notebooks.
Get Support: Everyone is encouraged to ask questions or raise issues in the Support Category of the Rubin Community Forum. Rubin staff will respond to all questions posted there.
1. Introduction¶
A set of timeseries features, also called variability characterization parameters, are computed for all rows of the DiaObject table using flux measurements from the DiaSource table.
This tutorial serves as a reference for interpreting the timeseries features of a given DiaObject.
For each feature, this tutorial provides descriptions and equations, code to recalculate the parameter using the DiaSource flux measurements, and an illustration of what the feature indicates about the lightcurve.
A known pulsating variable star is used as the example DiaObject.
It is important to note that although all timeseries features are calculated for all DiaObjects,
they are not all appropriate for the interpretation of all DiaObjects.
For example, a linear slope is calculated for all DiaObjects, but many types of transients and variables are not expected to exhibit a linear slope.
Why are there fewer features for DP2 than DP1? For DP2, the Data Release pipeline's calculation of timeseries features was adjusted to match the Prompt pipeline's codes for timeseries features. The intent is now that both pipelines will match, going forward, and build up the number of timeseries features.
For future data releases, additional timeseries features are likely to be included, as discussed in Data Management Tech Note 118 Review of Timeseries Features (DMTN-118).
Related tutorials: The 200-level tutorials on the DiaObject and DiaSource tables, and the other tutorial in the 207 series on timeseries features.
1.1. Features summary¶
The timeseries features are calculated based on the fluxes in the DiaSource table, and this table only includes visits (observations) for which there was a SNR $>5$ detection in the difference image.
For example, for periodic variables, visits for which the star had a brightness similar to its brightness in the template image - and thus very little flux in the difference image - do not contribute to the timeseries features' values.
The following parameters are calculated per filter (<f>), using the difference image photometry (psfFlux) column from the diaSources that are associated with a given DiaObject.
<f>_psfFluxNdata: The number of associateddiaSources.<f>_psfFluxErrMean: Mean of thediaSourcePSF flux errors.<f>_psfFluxMin,Max: Minimum and maximumdiaSourcePSF flux.<f>_psfFluxMean,MeanErr: Inverse variance weighted mean ofdiaSourcePSF flux, and its standard error.<f>_psfFluxSigma: Standard deviation of the distribution of<f>_psfFlux.<f>_psfFluxMaxSlope: Maximum ratio of time ordered $\Delta f / \Delta t$.
The following parameters are calculated per filter with forced photometry on the science image (also called the direct or visit image) at the diaSource positions (the scienceFlux column).
<f>_scienceFluxMean: Weighted mean of the PSF flux (forced photometered on the visit image).<f>_scienceFluxMeanErr: Standard error on<f>_scienceFluxMean.
Although the scienceFlux features are forced photometry, because they are from the DiaSource table they only include images in which the astrophysical object was detected with SNR$>5$ in the difference image.
These timeseries features do not use the forced photometry that is performed at the location of all diaObjects in all visit and difference images, which is stored in the ForcedSourceOnDiaObject table.
1.2. Import packages¶
Import numpy, a fundamental package for scientific computing with arrays in Python (numpy.org), and matplotlib, a comprehensive library for data visualization (matplotlib.org; matplotlib gallery).
Also import scipy.stats, a package containing a large number of statistical functions (scipy.stats).
From the lsst package, import modules for accessing the Table Access Protocol (TAP) service and making colorblind-friedly plots.
import numpy as np
import matplotlib.pyplot as plt
from lsst.rsp import RSPDiscovery
from lsst.utils.plotting import (get_multiband_plot_colors,
get_multiband_plot_symbols,
get_multiband_plot_linestyles)
1.3. Define parameters¶
Define colors, symbols, and linestyles to represent the six LSST filters, $ugrizy$.
f_col = get_multiband_plot_colors()
f_name = f_col.keys()
f_sym = get_multiband_plot_symbols()
f_lin = get_multiband_plot_linestyles()
Get an instance of the TAP service.
discovery = RSPDiscovery("dp2")
tap_service = discovery.get_tap_client()
For this tutorial, use only the $g$- and $r$-band to illustrate the timeseries features.
use_filters = ['g', 'r']
2. Retrieve data for a known variable star¶
The first step to obtain light curve data is to define the coordinates (ra and dec) of a known RR Lyrae variable star captured in DP2: Gaia DR3 6005347144427786624.
This star is classified in Gaia DR3's vari_rrlyrae table as an RRab-type variable with a period of 0.5180638410804721 days.
star_diaObjectId = 743582465276248313
2.1. The diaObject record¶
Retrieve all of the $g$- and $r$-band timeseries features from the DiaObject table.
Create a comma-separated list of columns to retrieve from the DiaObject table.
Include the $g$- and $r$-band timeseries features, plus the nDiaSources column.
query = "SELECT column_name " \
"FROM tap_schema.columns " \
"WHERE table_name = 'dp2.DiaObject'"
results = tap_service.search(query).to_table()
columns_list = ''
for filt in use_filters:
for name in results['column_name']:
if name.find(filt + '_') == 0:
columns_list += name + ', '
columns_list += 'nDiaSources'
del query, results
Define a query to retrieve all columns in the list from the DiaObject table for the row corresponding to the variable star, and submit the query to the TAP service.
query = """SELECT {} FROM dp2.DiaObject WHERE diaObjectId = {}
""".format(columns_list, star_diaObjectId)
job = tap_service.submit_job(query)
job.run()
job.wait(phases=['COMPLETED', 'ERROR'])
print('Job phase is', job.phase)
if job.phase == 'ERROR':
job.raise_if_error()
del query
Job phase is COMPLETED
Retrieve the results from the TAP service as an astropy table named star_diaObject.
assert job.phase == 'COMPLETED'
star_diaObject = job.fetch_result().to_table()
assert len(star_diaObject) == 1
Option to view star_diaObject.
# star_diaObject
job.delete()
2.2. The diaSource records¶
In general, for scientific analyses of lightcurves it is recommended to use the forced flux measurements in the ForcedSourceOnDiaObject table.
However, because the timeseries features in the DiaObject table are calculated using the detection (unforced) flux measurements from the DiaSource table, this tutorial also uses them.
The key difference to be aware of is that the DiaSource table only contains flux measurements for visits in which the diaObject was detected with SNR$\geq5$ in the difference image, whereas the ForcedSourceOnDiaObject table contains forced flux measurements for all visits.
Define a query to retrieve the $g$- and $r$-band difference-image (psfFlux) and direct-image (scienceFlux) measurements, the flux errors, and the exposure time midpoint Modified Julian Date (midpointMjdTai) from the DiaSource table.
query = """SELECT psfFlux, psfFluxErr, scienceFlux, scienceFluxErr, band, midpointMjdTai
FROM dp2.DiaSource
WHERE (band = '{}' OR band = '{}')
AND diaObjectId = {}
""".format(str(use_filters[0]), str(use_filters[1]), str(star_diaObjectId))
job = tap_service.submit_job(query)
job.run()
job.wait(phases=['COMPLETED', 'ERROR'])
print('Job phase is', job.phase)
if job.phase == 'ERROR':
job.raise_if_error()
del query
Job phase is COMPLETED
Retrieve the results from the TAP service as an astropy table named star_diaSources.
assert job.phase == 'COMPLETED'
star_diaSources = job.fetch_result().to_table()
print(len(star_diaSources))
352
job.delete()
Visualize the measured fluxes from the DiaSource table as a lightcurve.
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(6, 4), sharex=True)
for filt in use_filters:
tx = np.where(star_diaSources['band'] == filt)[0]
ax1.errorbar(star_diaSources[tx]['midpointMjdTai'], star_diaSources[tx]['psfFlux'],
yerr=star_diaSources[tx]['psfFluxErr'], fmt=f_sym[filt], ms=7,
alpha=0.5, mew=0, color=f_col[filt], label=filt)
ax2.errorbar(star_diaSources[tx]['midpointMjdTai'], star_diaSources[tx]['scienceFlux'],
yerr=star_diaSources[tx]['scienceFluxErr'], fmt=f_sym[filt], ms=7,
alpha=0.5, mew=0, color=f_col[filt], label=filt)
del tx
ax1.set_ylabel('psfFlux [nJy]')
ax2.set_ylabel('scienceFlux [nJy]')
ax2.set_xlabel('MJD [d]')
ax1.legend(loc='upper left')
plt.tight_layout()
plt.show()
Figure 1: The
psfFlux(difference image flux; top) andscienceFlux(direct image flux; bottom) for all visits in which the star was detected with SNR>5 in the difference image, versus the MJD of the visit, for the $g$- and $r$-bands only. The errors bars on the flux are, for most data points, relatively too small to be seen.
3. Timeseries features¶
The following subsections demonstrate how each of the timeseries features in the DiaObject table are derived by recalculating them from the observations in the DiaSource table.
For scientific analyses it is not necessary to recalculate these features; this tutorial does it only as a demonstration.
All features are calculated from the PSF flux measured on the difference images, psfFlux, except the three features in the last sub-section which are calculated from the PSF flux measured on the science images, scienceFlux.
3.1. Ndata¶
The <f>_psfFluxNdata column is the number of rows in the DiaSource table that are associated with a given diaObject, for a given filter.
In other words, this is the number of difference-image detections of a given diaObject in a given band.
Recalculate the <f>_psfFluxNdata column using the DiaSource table, and confirm it matches the DiaObject table.
for f, filt in enumerate(use_filters):
diao_ndata = int(star_diaObject[filt+'_psfFluxNdata'])
tx = np.where(star_diaSources['band'] == filt)[0]
dias_ndata = len(tx)
print(filt + '-band diaObject diaSource')
print('Ndata %10i %10i' % (diao_ndata, dias_ndata))
if f == 0:
print(' ')
del diao_ndata, tx, dias_ndata
g-band diaObject diaSource Ndata 204 204 r-band diaObject diaSource Ndata 148 148
3.2. ErrMean¶
The <f>_psfFluxErrMean column is the unweighted average of the psfFluxErr values: $\frac{1}{N}\sum{e_{f}}$, where $e_f$ is the flux error (psfFluxErr) and $N$ is the number of flux measurements.
Recalculate the <f>_psfFluxErrMean column using the DiaSource table, and confirm it matches the DiaObject table.
for f, filt in enumerate(use_filters):
diao_emean = float(star_diaObject[filt + '_psfFluxErrMean'])
tx = np.where(star_diaSources['band'] == filt)[0]
dias_emean = np.sum(star_diaSources['psfFluxErr'][tx])/len(tx)
print(filt + '-band diaObject diaSource')
print('ErrMean %10.0f %10.0f ' % (diao_emean, dias_emean))
if f == 0:
print(' ')
del diao_emean, tx, dias_emean
g-band diaObject diaSource ErrMean 364 364 r-band diaObject diaSource ErrMean 527 527
3.3. Min, Max¶
The <f>_psfFluxMin and Max values are the lowest and highest values of the psfFlux column (difference-image fluxes) for the diaObject.
Recalculate the <f>_psfFluxMin and Max columns using the DiaSource table, and confirm they match the DiaObject table.
for f, filt in enumerate(use_filters):
diao_min = float(star_diaObject[filt + '_psfFluxMin'])
diao_max = float(star_diaObject[filt + '_psfFluxMax'])
tx = np.where(star_diaSources['band'] == filt)[0]
dias_min = np.min(star_diaSources['psfFlux'][tx])
dias_max = np.max(star_diaSources['psfFlux'][tx])
print(filt + '-band diaObject diaSource')
print('Min %10.0f %10.0f ' % (diao_min, dias_min))
print('Max %10.0f %10.0f ' % (diao_max, dias_max))
if f == 0:
print(' ')
del diao_min, diao_max, tx, dias_min, dias_max
g-band diaObject diaSource Min -8579 -8579 Max 178302 178302 r-band diaObject diaSource Min -7607 -7607 Max 79400 79400
3.4. Mean, MeanErr¶
The <f>_psfFluxMean is the weighted mean of the measured difference-image fluxes: $\bar{f_w} = \frac{\sum{f\times w}}{\sum{w}}$.
The weights ($w$) are the inverse of the square of the psfFluxErr ($e_f$) values: $w = \frac{1}{e_f^2}$.
The <f>_psfFluxMeanErr is the error on the weighted mean flux, and it is the inverse of the root of the sum of the weights: $\epsilon = \frac{1}{\sqrt{\sum{w}}}$.
Recalculate the <f>_psfFluxMean and MeanErr columns using the DiaSource table, and confirm they match the DiaObject table.
for f, filt in enumerate(use_filters):
diao_mean = float(star_diaObject[filt + '_psfFluxMean'])
diao_meane = float(star_diaObject[filt + '_psfFluxMeanErr'])
tx = np.where(star_diaSources['band'] == filt)[0]
dias_weights = 1.0/(star_diaSources['psfFluxErr'][tx]**2)
dias_wmean = np.sum(star_diaSources['psfFlux'][tx] * dias_weights)/np.sum(dias_weights)
dias_meane = 1.0 / np.sqrt(np.sum(dias_weights))
print(filt + '-band diaObject diaSource')
print('Mean %10.0f %10.0f ' % (diao_mean, dias_wmean))
print('MeanErr %10.2f %10.2f ' % (diao_meane, dias_meane))
if f == 0:
print(' ')
del diao_mean, diao_meane, tx, dias_weights, dias_wmean, dias_meane
g-band diaObject diaSource Mean 23140 23140 MeanErr 24.08 24.08 r-band diaObject diaSource Mean 20796 20796 MeanErr 36.40 36.40
3.5. Sigma¶
The <f>_psfFluxSigma is the standard deviation in the measured difference-image fluxes ($\sigma_f$).
Note that this feature uses the unweighted average flux ($\bar{f}$) in its calculation:
$\sigma_f = \sqrt{\frac{\sum{(f - \bar{f})^2}}{N-1}}$.
Recalculate the <f>_psfFluxSigma column using the DiaSource table, and confirm it matches the DiaObject table.
for f, filt in enumerate(use_filters):
diao_sigma = float(star_diaObject[filt + '_psfFluxSigma'])
tx = np.where(star_diaSources['band'] == filt)[0]
dias_mean = np.sum(star_diaSources['psfFlux'][tx])/len(tx)
dias_sigma = np.sqrt(np.sum((star_diaSources['psfFlux'][tx] - dias_mean)**2)/(len(tx)-1))
print(filt + '-band diaObject diaSource')
print('Sigma %10.0f %10.0f ' % (diao_sigma, dias_sigma))
if f == 0:
print(' ')
del diao_sigma, tx, dias_mean, dias_sigma
g-band diaObject diaSource Sigma 49076 49076 r-band diaObject diaSource Sigma 25754 25754
Show the minimum and maximum, mean and sigma, and median and MAD values as lines on the lightcurve and the flux distribution.
fig, ax = plt.subplots(2, 2, figsize=(8, 5), sharex='col')
s_lw = [1, 1, 1, 1, 1]
s_ls = ['solid', 'solid', 'dashed', 'dotted', 'dotted']
s_lb = ['min/max', None, 'mean', 'sigma', None]
for f, filt in enumerate(['g', 'r']):
tx = np.where(star_diaSources['band'] == filt)[0]
ax[f, 0].errorbar(star_diaSources['midpointMjdTai'][tx], star_diaSources['psfFlux'][tx],
yerr=star_diaSources['psfFluxErr'][tx], fmt=f_sym[filt],
ms=7, alpha=0.5, mew=0, color=f_col[filt], label=filt + ', star')
ax[f, 1].hist(star_diaSources['psfFlux'][tx], bins=30,
alpha=0.5, color=f_col[filt], label=filt)
s_vl = [star_diaObject[filt + '_psfFluxMin'], star_diaObject[filt + '_psfFluxMax'],
star_diaObject[filt + '_psfFluxMean'],
star_diaObject[filt + '_psfFluxMean'] - star_diaObject[filt + '_psfFluxSigma'],
star_diaObject[filt + '_psfFluxMean'] + star_diaObject[filt + '_psfFluxSigma']]
for s in range(len(s_vl)):
ax[f, 0].axhline(s_vl[s], lw=s_lw[s], ls=s_ls[s], color='grey', label=s_lb[s])
ax[f, 1].axvline(s_vl[s], lw=s_lw[s], ls=s_ls[s], color='grey', label=s_lb[s])
ax[f, 1].legend(bbox_to_anchor=(1.7, 1), loc='upper right')
ax[f, 0].set_ylabel('psfFlux [nJy]')
del tx, s_vl, s_lb
plt.suptitle('Statistical features on the lightcurve and flux distribution of a variable star')
plt.tight_layout()
plt.show()
Figure 2: The $g$- (top) and $r$-band (bottom) light curves (left) and flux distributions (right) for the variable star, with the minimum and maximum (solid), and mean and standard deviation (dashed and dotted) statistics overplotted as grey lines.
3.6. Max slope¶
The psfFluxMaxSlope is the maximum "instantaneous" slope, or the maximum ratio of the series of time-ordered values of $\Delta f / \Delta t$.
for f, filt in enumerate(use_filters):
diao_mm = float(star_diaObject[filt + '_psfFluxMaxSlope'])
tx = np.where((star_diaSources['band'] == filt)
& (~np.isnan(star_diaSources['psfFlux']))
& (~np.isnan(star_diaSources['psfFluxErr']))
& (~np.isnan(star_diaSources['midpointMjdTai'])))[0]
sx = np.argsort(star_diaSources['midpointMjdTai'][tx])
dias_mm = (np.diff(star_diaSources['psfFlux'][tx[sx]])
/ np.diff(star_diaSources['midpointMjdTai'][tx[sx]])).max()
print(filt + 'band diaObject diaSource')
print('Max Slope %10.0f %10.0f ' % (diao_mm, dias_mm))
if f == 0:
print(' ')
del diao_mm, tx, sx, dias_mm
gband diaObject diaSource Max Slope 11886967 11887349 rband diaObject diaSource Max Slope 6746364 6746335
Notice: For the linear fits, the re-derived parameters from the
DiaSourcetable are not an exact match to the parameters stored in theDiaObjecttable. It is unclear why as the code above uses the same calculation as in the diaCalculationPlugin for classLinearFitDiaPsfFlux.
3.7. Science flux features¶
The <f>_scienceFluxMean is the weighted mean of the forced PSF fluxes measured on the science (direct) images.
The <f>_scienceFluxMeanErr is the error on the weighted mean flux.
Both features use the same formulae as quoted in Section 3.4.
Demonstrate how to derive these statistics for the scienceFlux.
for f, filt in enumerate(use_filters):
diao_sci_mean = float(star_diaObject[filt + '_scienceFluxMean'])
diao_sci_meane = float(star_diaObject[filt + '_scienceFluxMeanErr'])
tx = np.where(star_diaSources['band'] == filt)[0]
dias_sci_w = 1.0/(star_diaSources['scienceFluxErr'][tx]**2)
dias_sci_wmean = np.sum(star_diaSources['scienceFlux'][tx] * dias_sci_w) / np.sum(dias_sci_w)
dias_sci_meane = 1.0 / np.sqrt(np.sum(dias_sci_w))
print(filt + '-band diaObject diaSource')
print('Mean %10.0f %10.0f ' % (diao_sci_mean, dias_sci_wmean))
print('MeanErr %10.2f %10.2f ' % (diao_sci_meane, dias_sci_meane))
if f == 0:
print(' ')
del diao_sci_mean, diao_sci_meane
del tx, dias_sci_w, dias_sci_wmean, dias_sci_meane
g-band diaObject diaSource Mean 120780 120780 MeanErr 24.19 24.19 r-band diaObject diaSource Mean 153782 153782 MeanErr 36.12 36.12