Skip to content

Commit c5d97ea

Browse files
MelissaGrahamMelissaGraham
authored
Merge pull request #116 from lsst/tickets/SP-2880
SP-2880: incorporate review comments into NBs 207 1, 2
2 parents 9efccb4 + d2e5f0c commit c5d97ea

2 files changed

Lines changed: 81 additions & 49 deletions

File tree

‎DP1/200_Data_Products/207_Timeseries_features/207_1_Timeseries_values.ipynb‎

Lines changed: 71 additions & 43 deletions
Original file line numberDiff line numberDiff line change
@@ -22,7 +22,7 @@
2222
"Data Release: [Data Preview 1](https://dp1.lsst.io)\\\n",
2323
"Container Size: Large\\\n",
2424
"LSST Science Pipelines version: r29.2.0\\\n",
25-
"Last verified to run: 2025-12-30\\\n",
25+
"Last verified to run: 2026-02-03\\\n",
2626
"Repository: [github.com/lsst/tutorial-notebooks](https://github.com/lsst/tutorial-notebooks)\\\n",
2727
"DOI: [10.11578/rubin/dc.20250909.20](https://doi.org/10.11578/rubin/dc.20250909.20)"
2828
]
@@ -85,7 +85,7 @@
8585
"* `<f>_psfFluxNdata` : The number of associated `diaSources`.\n",
8686
"* `<f>_psfFluxErrMean` : Mean of the `diaSource` PSF flux errors.\n",
8787
"* `<f>_psfFluxMin`, `Max` : Minimum and maximum `diaSource` PSF flux.\n",
88-
"* `<f>_psfFluxMean`, `MeanErr` : Weighted mean of `diaSource` PSF flux, and its standard error.\n",
88+
"* `<f>_psfFluxMean`, `MeanErr` : Inverse variance weighted mean of `diaSource` PSF flux, and its standard error.\n",
8989
"* `<f>_psfFluxSigma` : Standard deviation of the distribution of `<f>_psfFlux`.\n",
9090
"* `<f>_psfFluxMAD` : Median absolute deviation (MAD) of `diaSource` PSF flux. Does not include scale factor for comparison to sigma.\n",
9191
"* `<f>_psfFluxChi2` : $\\chi^2$ statistic for the scatter of `<f>_psfFlux` around `<f>_psfFluxMean`.\n",
@@ -339,7 +339,7 @@
339339
"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.\n",
340340
"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.\n",
341341
"\n",
342-
"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."
342+
"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."
343343
]
344344
},
345345
{
@@ -592,7 +592,7 @@
592592
"Note that this feature uses the *unweighted* average flux ($\\bar{f}$) in its calculation:\n",
593593
"$\\sigma_f = \\sqrt{\\frac{\\sum{(f - \\bar{f})^2}}{N-1}}$.\n",
594594
"\n",
595-
"The `<f>_psfFluxMAD`, the median absolute deviation (MAD), is the median value of an array composed of the absolute values of the differences between the `psfFlux` and the median value of the `psfFlux`, for a given filter.\n",
595+
"The `<f>_psfFluxMAD`, the median absolute deviation (MAD), is the median value of an array composed of the absolute values of the differences between the `psfFlux` and the median value of the `psfFlux`, for a given filter, where the median flux is `<f>_psfFluxPercentile50` (see Section 3.8).\n",
596596
"\n",
597597
"Recalculate the `<f>_psfFluxSigma` and `MAD` columns using the `DiaSource` table, and confirm they match the `DiaObject` table."
598598
]
@@ -625,42 +625,59 @@
625625
"id": "8ba9a97d-2df3-48c4-92af-ebf95b8444ea",
626626
"metadata": {},
627627
"source": [
628-
"Show the minimum, maximum, mean, sigma, and MAD values as lines on the lightcurve and the flux distribution.\n",
629-
"\n",
630-
"Define the labels, linewidths, and linestyles of the lines to represent the statistical features, then create the plot."
628+
"Show the minimum and maximum, mean and sigma, and median and MAD values as lines on the lightcurve and the flux distribution."
631629
]
632630
},
633631
{
634632
"cell_type": "code",
635633
"execution_count": null,
636-
"id": "e7dc1ab3-1bf5-4acb-a962-2b18735734f8",
634+
"id": "3484a8a8-bb89-4b0d-b19d-dd97997f12c4",
637635
"metadata": {},
638636
"outputs": [],
639637
"source": [
640-
"s_lb = ['min/max', None, 'mean', 'sigma', None, 'MAD', None]\n",
641-
"s_lw = [1, 1, 1, 2, 2, 1, 1]\n",
642-
"s_ls = ['dashed', 'dashed', 'solid', 'dotted', 'dotted', '-.', '-.']\n",
643638
"fig, ax = plt.subplots(2, 2, figsize=(8, 5), sharex='col')\n",
644-
"for f, filt in enumerate(use_filters):\n",
645-
" tx = np.where(star_diaSources['band'] == filt)[0]\n",
646-
" ax[f, 0].errorbar(star_diaSources['midpointMjdTai'][tx], star_diaSources['psfFlux'][tx],\n",
647-
" yerr=star_diaSources['psfFluxErr'][tx], fmt=f_sym[filt],\n",
648-
" ms=7, alpha=0.5, mew=0, color=f_col[filt], label=filt + ', star')\n",
649-
" ax[f, 1].hist(star_diaSources['psfFlux'][tx], bins=30,\n",
650-
" alpha=0.5, color=f_col[filt], label=filt)\n",
651-
" s_vl = [star_diaObject[filt + '_psfFluxMin'], star_diaObject[filt + '_psfFluxMax'],\n",
652-
" star_diaObject[filt + '_psfFluxMean'],\n",
653-
" star_diaObject[filt + '_psfFluxMean'] - star_diaObject[filt + '_psfFluxSigma'],\n",
654-
" star_diaObject[filt + '_psfFluxMean'] + star_diaObject[filt + '_psfFluxSigma'],\n",
655-
" star_diaObject[filt + '_psfFluxMean'] - star_diaObject[filt + '_psfFluxMAD'],\n",
656-
" star_diaObject[filt + '_psfFluxMean'] + star_diaObject[filt + '_psfFluxMAD']]\n",
657-
" for s in range(len(s_vl)):\n",
658-
" ax[f, 0].axhline(s_vl[s], lw=s_lw[s], ls=s_ls[s], color='grey', label=s_lb[s])\n",
659-
" ax[f, 1].axvline(s_vl[s], lw=s_lw[s], ls=s_ls[s], color='grey', label=s_lb[s])\n",
660-
" ax[f, 0].set_ylabel('psfFlux [nJy]')\n",
661-
" ax[f, 1].set_ylabel('N(diaSources)')\n",
662-
" ax[f, 1].legend(bbox_to_anchor=(1.7, 1), loc='upper right')\n",
663-
" del tx, s_vl\n",
639+
"\n",
640+
"f = 0\n",
641+
"filt = 'g'\n",
642+
"s_lw = [1, 1, 1, 1, 1]\n",
643+
"s_ls = ['solid', 'solid', 'dashed', 'dotted', 'dotted']\n",
644+
"\n",
645+
"s_lb = ['min/max', None, 'mean', 'sigma', None]\n",
646+
"tx = np.where(star_diaSources['band'] == filt)[0]\n",
647+
"ax[0, 0].errorbar(star_diaSources['midpointMjdTai'][tx], star_diaSources['psfFlux'][tx],\n",
648+
" yerr=star_diaSources['psfFluxErr'][tx], fmt=f_sym[filt],\n",
649+
" ms=7, alpha=0.5, mew=0, color=f_col[filt], label=filt + ', star')\n",
650+
"ax[0, 1].hist(star_diaSources['psfFlux'][tx], bins=30,\n",
651+
" alpha=0.5, color=f_col[filt], label=filt)\n",
652+
"s_vl = [star_diaObject[filt + '_psfFluxMin'], star_diaObject[filt + '_psfFluxMax'],\n",
653+
" star_diaObject[filt + '_psfFluxMean'],\n",
654+
" star_diaObject[filt + '_psfFluxMean'] - star_diaObject[filt + '_psfFluxSigma'],\n",
655+
" star_diaObject[filt + '_psfFluxMean'] + star_diaObject[filt + '_psfFluxSigma']]\n",
656+
"for s in range(len(s_vl)):\n",
657+
" ax[0, 0].axhline(s_vl[s], lw=s_lw[s], ls=s_ls[s], color='grey', label=s_lb[s])\n",
658+
" ax[0, 1].axvline(s_vl[s], lw=s_lw[s], ls=s_ls[s], color='grey', label=s_lb[s])\n",
659+
"ax[0, 1].legend(bbox_to_anchor=(1.7, 1), loc='upper right')\n",
660+
"ax[0, 0].set_ylabel('psfFlux [nJy]')\n",
661+
"del tx, s_vl, s_lb\n",
662+
"\n",
663+
"s_lb = ['min/max', None, 'median', 'MAD', None]\n",
664+
"tx = np.where(star_diaSources['band'] == filt)[0]\n",
665+
"ax[1, 0].errorbar(star_diaSources['midpointMjdTai'][tx], star_diaSources['psfFlux'][tx],\n",
666+
" yerr=star_diaSources['psfFluxErr'][tx], fmt=f_sym[filt],\n",
667+
" ms=7, alpha=0.5, mew=0, color=f_col[filt], label=filt + ', star')\n",
668+
"ax[1, 1].hist(star_diaSources['psfFlux'][tx], bins=30,\n",
669+
" alpha=0.5, color=f_col[filt], label=filt)\n",
670+
"s_vl = [star_diaObject[filt + '_psfFluxMin'], star_diaObject[filt + '_psfFluxMax'],\n",
671+
" star_diaObject[filt + '_psfFluxPercentile50'],\n",
672+
" star_diaObject[filt + '_psfFluxPercentile50'] - star_diaObject[filt + '_psfFluxMAD'],\n",
673+
" star_diaObject[filt + '_psfFluxPercentile50'] + star_diaObject[filt + '_psfFluxMAD']]\n",
674+
"for s in range(len(s_vl)):\n",
675+
" ax[1, 0].axhline(s_vl[s], lw=s_lw[s], ls=s_ls[s], color='grey', label=s_lb[s])\n",
676+
" ax[1, 1].axvline(s_vl[s], lw=s_lw[s], ls=s_ls[s], color='grey', label=s_lb[s])\n",
677+
"ax[1, 0].set_ylabel('psfFlux [nJy]')\n",
678+
"ax[1, 1].set_ylabel('N(diaSources)')\n",
679+
"ax[1, 1].legend(bbox_to_anchor=(1.7, 1), loc='upper right')\n",
680+
"del tx, s_vl\n",
664681
"ax[1, 0].set_xlabel('MJD [d]')\n",
665682
"ax[1, 1].set_xlabel('psfFlux [nJy]')\n",
666683
"plt.suptitle('Statistical features on the lightcurve and flux distribution of a variable star')\n",
@@ -674,7 +691,7 @@
674691
"id": "42e723a1-be50-4aa2-978b-1b76a4985486",
675692
"metadata": {},
676693
"source": [
677-
"> **Figure 2:** The $g$- and $r$-band lightcurves (left) and flux distributions (right) for the variable star, with the minimum and maximum (dashed), mean (solid), and standard deviation (sigma; dotted) of the flux values overplotted as grey lines."
694+
"> **Figure 2:** The $g$-band lightcurves (left) and flux distributions (right) for the variable star, with the minimum and maximum (solid), mean and standard deviation (top row; dashed and dotted), and median and MAD (bottom row; dashed and dotted) of the flux values overplotted as grey lines. "
678695
]
679696
},
680697
{
@@ -687,7 +704,7 @@
687704
"The `<f>_psfFluxChi2` for a given filter is $\\chi^2 = \\sum\\left(\\frac{f - \\bar{f}}{e_f}\\right)^2$,\n",
688705
"where $f$ is the `psfFlux`, $\\bar{f}$ is the weighted mean of the `psfFlux`, and $e_f$ is the `psfFluxErr`.\n",
689706
"\n",
690-
"Notice that the equation for $\\chi^2$ is quite similar to the one for $\\sigma_f$, except in $\\chi^2$ the flux differences from the mean are weighted by the flux error.\n",
707+
"Notice that the equation for $\\chi^2$ is quite similar to the one for $\\sigma_f$, except in $\\chi^2$ the flux differences from the mean are weighted by the flux error, whereas $\\sigma_f$ uses a denominator of $N-1$.\n",
691708
"\n",
692709
"In a sense, the $\\chi^2$ is a measure of whether the flux variability is more than that expected from random variations due to measurement uncertainty.\n",
693710
"\n",
@@ -966,28 +983,31 @@
966983
"\n",
967984
"The `<f>_psfFluxLinearSlope` and `<f>_psfFluxLinearIntercept` are the result of fitting a linear regression using the `scipy.optimize.lsq_linear` function.\n",
968985
"\n",
969-
"The `psfFluxMaxSlope` is the maximum \"instantaneous\" slope, or the maximum ration of the series of time-ordered values of $\\Delta f / \\Delta t$.\n",
986+
"The `psfFluxMaxSlope` is the maximum \"instantaneous\" slope, or the maximum ratio of the series of time-ordered values of $\\Delta f / \\Delta t$.\n",
970987
"\n",
971988
"Demonstrate how to derive the linear fit parameters."
972989
]
973990
},
974991
{
975992
"cell_type": "code",
976993
"execution_count": null,
977-
"id": "9b8e89d5-a363-4095-a70b-9d670621946f",
994+
"id": "b3ea4d9a-5816-4dc9-bd64-31494f790e7b",
978995
"metadata": {},
979996
"outputs": [],
980997
"source": [
981998
"for f, filt in enumerate(use_filters):\n",
982999
" diao_b = float(star_diaObject[filt + '_psfFluxLinearIntercept'])\n",
9831000
" diao_m = float(star_diaObject[filt + '_psfFluxLinearSlope'])\n",
9841001
" diao_mm = float(star_diaObject[filt + '_psfFluxMaxSlope'])\n",
985-
" tx = np.where(star_diaSources['band'] == filt)[0]\n",
986-
" A = np.array([star_diaSources['midpointMjdTai'][tx]\n",
987-
" / star_diaSources['psfFluxErr'][tx],\n",
988-
" 1.0 / star_diaSources['psfFluxErr'][tx]]).transpose()\n",
989-
" dias_m, dias_b = lsq_linear(A, star_diaSources['psfFlux'][tx]\n",
990-
" / star_diaSources['psfFluxErr'][tx]).x\n",
1002+
" tx = np.where((star_diaSources['band'] == filt)\n",
1003+
" & (~np.isnan(star_diaSources['psfFlux']))\n",
1004+
" & (~np.isnan(star_diaSources['psfFluxErr']))\n",
1005+
" & (~np.isnan(star_diaSources['midpointMjdTai'])))[0]\n",
1006+
" fluxes = star_diaSources['psfFlux'][tx]\n",
1007+
" errors = star_diaSources['psfFluxErr'][tx]\n",
1008+
" times = star_diaSources['midpointMjdTai'][tx]\n",
1009+
" A = np.array([times / errors, 1 / errors]).transpose()\n",
1010+
" dias_m, dias_b = lsq_linear(A, fluxes / errors).x\n",
9911011
" sx = np.argsort(star_diaSources['midpointMjdTai'][tx])\n",
9921012
" dias_mm = (np.diff(star_diaSources['psfFlux'][tx[sx]])\n",
9931013
" / np.diff(star_diaSources['midpointMjdTai'][tx[sx]])).max()\n",
@@ -997,15 +1017,15 @@
9971017
" print('Max Slope %10.0f %10.0f ' % (diao_mm, dias_mm))\n",
9981018
" if f == 0:\n",
9991019
" print(' ')\n",
1000-
" del diao_b, diao_m, diao_mm, tx, A, dias_b, dias_m, sx, dias_mm"
1020+
" del diao_b, diao_m, diao_mm, tx, fluxes, errors, times, A, dias_b, dias_m, sx, dias_mm"
10011021
]
10021022
},
10031023
{
10041024
"cell_type": "markdown",
10051025
"id": "515aaa7f-6566-4a71-b85c-a73509522900",
10061026
"metadata": {},
10071027
"source": [
1008-
"**> Notice:** For the linear fits, the re-derived parameters from the `DiaSource` table are not an exact match to the parameters stored in the `DiaObject` table.\n",
1028+
">**Notice:** For the linear fits, the re-derived parameters from the `DiaSource` table are not an exact match to the parameters stored in the `DiaObject` table. It is unclear why as the code above uses the same calculation as in the [diaCalculationPlugin](https://github.com/lsst/meas_base/blob/main/python/lsst/meas/base/diaCalculationPlugins.py) for class `LinearFitDiaPsfFlux`.\n",
10091029
"\n",
10101030
"The known variable star is not expected to exhibit a linear slope.\n",
10111031
"Plot the lightcurve and show the best fit line."
@@ -1099,6 +1119,14 @@
10991119
"metadata": {},
11001120
"outputs": [],
11011121
"source": []
1122+
},
1123+
{
1124+
"cell_type": "code",
1125+
"execution_count": null,
1126+
"id": "22f49979-af1e-4884-9436-e4b44d143e8a",
1127+
"metadata": {},
1128+
"outputs": [],
1129+
"source": []
11021130
}
11031131
],
11041132
"metadata": {

‎DP1/200_Data_Products/207_Timeseries_features/207_2_Timeseries_distributions.ipynb‎

Lines changed: 10 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -22,7 +22,7 @@
2222
"Data Release: [Data Preview 1](https://dp1.lsst.io)\\\n",
2323
"Container Size: Large\\\n",
2424
"LSST Science Pipelines version: r29.2.0\\\n",
25-
"Last verified to run: 2025-12-30\\\n",
25+
"Last verified to run: 2026-02-03\\\n",
2626
"Repository: [github.com/lsst/tutorial-notebooks](https://github.com/lsst/tutorial-notebooks)\\\n",
2727
"DOI: [10.11578/rubin/dc.20250909.20](https://doi.org/10.11578/rubin/dc.20250909.20)"
2828
]
@@ -326,7 +326,11 @@
326326
"id": "16eeb034-6728-406e-bfc8-7b875c3a9416",
327327
"metadata": {},
328328
"source": [
329-
"Retrieve the $r$-band difference-image (`psfFlux`) and direct-image (`scienceFlux`) measurements from the `DiaSource` table, their errors, and the midpoint of the exposure time as a Modified Julian Date (MJD)."
329+
"Retrieve the $r$-band difference-image (`psfFlux`) and direct-image (`scienceFlux`) measurements from the `DiaSource` table, their errors, and the midpoint of the exposure time as a Modified Julian Date (MJD).\n",
330+
"\n",
331+
">**Notice:** In general, for scientific analyses of lightcurves it is recommended to use the *forced* flux measurements in the `ForcedSourceOnDiaObject` table.\n",
332+
"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.\n",
333+
"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."
330334
]
331335
},
332336
{
@@ -527,7 +531,7 @@
527531
"source": [
528532
"> **Figure 2**: The `psfFlux` (difference image flux; top) and `scienceFlux` (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 $r$-band. The errors bars on the flux are, for most data points, relatively too small to be seen.\n",
529533
"\n",
530-
"Notice that visits for which this object's brightness was similar to that of its brightness in the template image are missing from the lightcurve above (i.e., visits around MJD $\\sim 606404$ days, when `psfFlux` was near 0 nJy).\n",
534+
"Notice that visits for which this object's brightness was similar to that of its brightness in the template image are missing from the lightcurve above (i.e., visits around MJD $\\sim 60640$ days, when `psfFlux` was near 0 nJy).\n",
531535
"\n",
532536
"Notice also that because this object is an astrophysical transient embedded in its host galaxy, the `scienceFlux` measurements are contaminated by host galaxy light.\n",
533537
"This contamination is not simply a \"pedestal\" (i.e., are not a constant amount added to the `psfFlux`) because the forced photometry uses the PSF of the image.\n",
@@ -542,7 +546,7 @@
542546
"source": [
543547
"### 2.3. Samples of diaObjects\n",
544548
"\n",
545-
"Retrieve two samples of `diaObjects`: one from the same DP1 field as the star (the \"Low-$b$\" field), and one from the same DP1 field as the transient (the ECDFS field).\n",
549+
"Retrieve two samples of `diaObjects`: one from the same DP1 field as the star (the low Galactic latitude or \"Low-$b$\" field), and one from the same DP1 field as the transient (the ECDFS field).\n",
546550
"\n",
547551
"To define the cone search, use the coordinates of the star and transient, and a 2.0 deg radius, which will encompass the relevant DP1 field.\n",
548552
"Print the coordinates of the star and transient."
@@ -589,7 +593,7 @@
589593
"metadata": {},
590594
"source": [
591595
"Apply the constraints in the ADQL query statement, and add a requirement that only `diaObjects` which were detected at least 10 times in the $r$-band be retrieved.\n",
592-
"This will constraint will mean that the timeseries features used in this tutorial are all measured from at least 10 detections.\n",
596+
"This constraint will mean that the timeseries features used in this tutorial are all measured from at least 10 detections.\n",
593597
"\n",
594598
"Create the query for the star's field (the \"Low-$b$\" field) and submit the query job."
595599
]
@@ -716,7 +720,7 @@
716720
"id": "8c66e3e1-2aec-464a-9e7d-36a8782cc36e",
717721
"metadata": {},
718722
"source": [
719-
"> **Figure 3:** The number of `diaObjects` with a given number of `diaSources` (number of difference image detections) in the $r$-band. While the star is in the tail of its population's distribution, the transient is not."
723+
"> **Figure 3:** The number of `diaObjects` with a given number of `diaSources` (number of difference image detections) in the $r$-band. The vertical lines mark the values for the star and transient respectively (as in the legend). While the star is in the tail of its population's distribution, the transient is not."
720724
]
721725
},
722726
{

0 commit comments

Comments
 (0)