Merge pull request #39 from lsst/u/jcarlin/fsodo_300_fix

Update queries in lightcurve NBs

Author: jeffcarlin

DP1/300_Science_Demos/305_Galactic_variables_and_transients/305_1_Variable_star_light_curves.ipynb

@@ -18,7 +18,7 @@
     "For the Rubin Science Platform at data.lsst.cloud. <br>\n",
     "Data Release: <a href=\"https://dp1.lsst.io/\">Data Preview 1</a> <br>\n",
     "Container Size: large <br>\n",
-    "LSST Science Pipelines version: Release r29.1.1 <br>\n",
+    "LSST Science Pipelines version: r29.1.1 <br>\n",
     "Last verified to run: 2025-06-20 <br>\n",
     "Repository: <a href=\"https://github.com/lsst/tutorial-notebooks\">github.com/lsst/tutorial-notebooks</a> <br>"
    ]
@@ -30,7 +30,7 @@
    "source": [
     "**Learning objective:** Learn how to extract and plot a light curve for a variable star from DP1.\n",
     "\n",
-    "**LSST data products:** `ForcedSourceOnDiaObject`, `Visit`\n",
+    "**LSST data products:** `ForcedSourceOnDiaObject`, `DiaObject`, `Visit`\n",
     "\n",
     "**Packages:** `matplotlib`, `numpy`, `lsst.rsp`, `lsst.utils.plotting`\n",
     "\n",
@@ -156,13 +156,15 @@
     "\n",
     "The units used to plot the light curve are in flux (nJy) and not magnitudes since the negative flux measurements would be omitted when converting to magnitudes.\n",
     "\n",
+    "### 2.1 Look up diaObjectId from coordinates\n",
+    "\n",
     "The first step to obtain light curve data is to define the coordinates (ra and dec) of a variable star captured in DP1: Gaia DR3 2912281258855051520."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "35077a54-59bd-44fa-9d31-114e3ddbd578",
+   "id": "2b4fc159-bca1-40e8-b9b1-58b39a151bf8",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -172,7 +174,7 @@
   },
   {
    "cell_type": "markdown",
-   "id": "1f486e43-cbbe-4cab-b131-27eadc69c3b5",
+   "id": "2d09a784-726c-4723-9df6-2f6ad99f6396",
    "metadata": {},
    "source": [
     "Define a 0.5 arcsecond search radius, and convert it to degrees."
@@ -181,7 +183,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "fe86ad87-b4ee-45cc-ad6c-71d8537c249e",
+   "id": "29a7fea7-0329-419d-af1f-fbd9b5e9eb4f",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -190,31 +192,105 @@
   },
   {
    "cell_type": "markdown",
-   "id": "3154ec3a-cd86-45ad-b215-0946624902b7",
+   "id": "ba3dc2e2-aa27-49bb-a4de-446538973955",
+   "metadata": {},
+   "source": [
+    "Search at this position in the `DiaObject` table to retrieve the `diaObjectId` of the corresponding object."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "11e3579b-6812-479a-916b-fdad65e01baa",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "query = \"SELECT ra, dec, diaObjectId \"\\\n",
+    "        \"FROM dp1.DiaObject \"\\\n",
+    "        \"WHERE CONTAINS(POINT('ICRS', ra, dec), \"\\\n",
+    "        \"CIRCLE('ICRS', \" + str(ra_var) + \", \"\\\n",
+    "        + str(dec_var) + \", \" + str(search_rad) + \")) = 1\"\n",
+    "print(query)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "08a4d843-d206-4f35-b635-e4f9253603d7",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "job = service.submit_job(query)\n",
+    "job.run()\n",
+    "job.wait(phases=['COMPLETED', 'ERROR'])\n",
+    "print('Job phase is', job.phase)\n",
+    "if job.phase == 'ERROR':\n",
+    "    job.raise_if_error()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "f14304ef-f1b2-4d2e-9b2f-c8099a592cc8",
+   "metadata": {
+    "execution": {
+     "iopub.execute_input": "2025-06-30T22:25:25.363571Z",
+     "iopub.status.busy": "2025-06-30T22:25:25.362713Z",
+     "iopub.status.idle": "2025-06-30T22:25:25.368671Z",
+     "shell.execute_reply": "2025-06-30T22:25:25.367868Z",
+     "shell.execute_reply.started": "2025-06-30T22:25:25.363544Z"
+    }
+   },
+   "source": [
+    "Fetch the results of the TAP search in table form and print the unique entries for `diaObjectId` in the output table."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9b7285dd-4b63-4482-baf0-8f63650f8834",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "assert job.phase == 'COMPLETED'\n",
+    "DiaObj = job.fetch_result().to_table()\n",
+    "print(np.unique(DiaObj['diaObjectId']))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "4f3fcea5-957d-4984-8f93-e49e1970e48b",
    "metadata": {},
    "source": [
-    "Perform a TAP search on the coordinates of the transient in the `ForcedSourceOnDiaObject` table.\n",
+    "### 2.2. Search based on diaObjectId\n",
     "\n",
-    "Use a table JOIN with the `Visit` table to extract visit timing info (`expMidptMJD`) for each entry in the `ForcedSourceOnDiaObject`."
+    "Provide this `diaObjectId` as a constraint to a TAP query to obtain the time series from `ForcedSourceOnDiaObject`. The use of an ID in this manner is always the recommended way to access data from a ForcedSource table."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "aee30c7f-38cf-4637-a7e9-dd476a69ddbc",
+   "id": "0ca36902-76c6-4acd-a02b-c7e609f2cb82",
    "metadata": {},
    "outputs": [],
    "source": [
-    "query = \"SELECT fsodo.diaObjectId, fsodo.coord_ra, fsodo.coord_dec, \"\\\n",
-    "        \"fsodo.visit, fsodo.detector, fsodo.band, \"\\\n",
+    "DiaObjID = DiaObj['diaObjectId'][0]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3fb7afdc-9709-4d8e-a3e6-fa46075fe0ef",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "query = \"SELECT fsodo.coord_ra, fsodo.coord_dec, \"\\\n",
+    "        \"fsodo.diaObjectId, fsodo.visit, fsodo.band, \"\\\n",
     "        \"fsodo.psfDiffFlux, fsodo.psfDiffFluxErr, \"\\\n",
     "        \"fsodo.psfFlux as psfFlux, fsodo.psfFluxErr, \"\\\n",
     "        \"vis.expMidptMJD \"\\\n",
     "        \"FROM dp1.ForcedSourceOnDiaObject as fsodo \"\\\n",
     "        \"JOIN dp1.Visit as vis ON vis.visit = fsodo.visit \"\\\n",
-    "        \"WHERE CONTAINS (POINT('ICRS', coord_ra, coord_dec), \"\\\n",
-    "        \"CIRCLE('ICRS', \" + str(ra_var) + \", \"\\\n",
-    "        + str(dec_var) + \", \" + str(search_rad) + \")) = 1 \"\n",
+    "        \"WHERE fsodo.diaObjectId = \"+str(DiaObjID)\n",
     "print(query)"
    ]
   },
@@ -223,13 +299,13 @@
    "id": "600390ba-63e6-47b6-a91b-1c19a4f5ff87",
    "metadata": {},
    "source": [
-    "Run the TAP query, then extract the results of the TAP search to an Astropy table."
+    "Run the TAP query."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "8685e9ba-5206-4ab9-9317-4f581ace2144",
+   "id": "fca45813-9e8f-4571-ac51-ca10dac2046a",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -238,18 +314,27 @@
     "job.wait(phases=['COMPLETED', 'ERROR'])\n",
     "print('Job phase is', job.phase)\n",
     "if job.phase == 'ERROR':\n",
-    "    job.raise_if_error()\n",
-    "assert job.phase == 'COMPLETED'"
+    "    job.raise_if_error()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "24628f61-7888-4df3-a72e-d2b18c6963c8",
+   "metadata": {},
+   "source": [
+    "Fetch the results of the TAP search in table form and print the number of rows."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "b60aec9d-54ff-43fb-bacf-4d498b453887",
+   "id": "8685e9ba-5206-4ab9-9317-4f581ace2144",
    "metadata": {},
    "outputs": [],
    "source": [
-    "forced_sources = job.fetch_result().to_table()"
+    "assert job.phase == 'COMPLETED'\n",
+    "forced_source = job.fetch_result().to_table()\n",
+    "print(len(forced_source))"
    ]
   },
   {
@@ -267,7 +352,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# forced_sources"
+    "# forced_source"
    ]
   },
   {
@@ -286,7 +371,7 @@
    "outputs": [],
    "source": [
     "for band in filter_names:\n",
-    "    print(f\"{band}: {np.sum((forced_sources['band'] == band))} visits\")"
+    "    print(f\"{band}: {np.sum((forced_source['band'] == band))} visits\")"
    ]
   },
   {
@@ -312,11 +397,11 @@
     "fig.subplots_adjust(wspace=0, hspace=0)\n",
     "\n",
     "for band in filter_names:\n",
-    "    pickband = np.where(forced_sources['band'] == band)\n",
-    "    ax1.plot(forced_sources[pickband]['expMidptMJD'], forced_sources[pickband]['psfDiffFlux'],\n",
+    "    pickband = np.where(forced_source['band'] == band)\n",
+    "    ax1.plot(forced_source[pickband]['expMidptMJD'], forced_source[pickband]['psfDiffFlux'],\n",
     "             marker=filter_symbols[band], color=filter_colors[band],\n",
     "             label=f\"{band}\", markersize=7, linestyle='none', fillstyle='none')\n",
-    "    ax2.plot(forced_sources[pickband]['expMidptMJD'], forced_sources[pickband]['psfFlux'],\n",
+    "    ax2.plot(forced_source[pickband]['expMidptMJD'], forced_source[pickband]['psfFlux'],\n",
     "             marker=filter_symbols[band], color=filter_colors[band], linestyle='none',\n",
     "             markersize=7, fillstyle='none')\n",
     "\n",
@@ -374,10 +459,10 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "pickr = np.where(forced_sources['band'] == 'r')\n",
-    "ind_max = np.argmax(forced_sources['psfFlux'][pickr])\n",
-    "t0 = forced_sources['expMidptMJD'][pickr][ind_max]\n",
-    "mjd_norm = (forced_sources['expMidptMJD'] - t0) / period\n",
+    "pickr = np.where(forced_source['band'] == 'r')\n",
+    "ind_max = np.argmax(forced_source['psfFlux'][pickr])\n",
+    "t0 = forced_source['expMidptMJD'][pickr][ind_max]\n",
+    "mjd_norm = (forced_source['expMidptMJD'] - t0) / period\n",
     "phase = np.mod(mjd_norm, 1.0)"
    ]
   },
@@ -400,12 +485,12 @@
     "fig.subplots_adjust(wspace=0, hspace=0)\n",
     "\n",
     "for band in filter_names:\n",
-    "    pickband = np.where((forced_sources['band'] == band)\n",
-    "                        & (forced_sources['psfDiffFlux'] > -50000))\n",
-    "    ax1.plot(phase[pickband], forced_sources[pickband]['psfDiffFlux'],\n",
+    "    pickband = np.where((forced_source['band'] == band)\n",
+    "                        & (forced_source['psfDiffFlux'] > -50000))\n",
+    "    ax1.plot(phase[pickband], forced_source[pickband]['psfDiffFlux'],\n",
     "             marker=filter_symbols[band], color=filter_colors[band],\n",
     "             label=f\"{band}\", markersize=7, linestyle='none', fillstyle='none')\n",
-    "    ax2.plot(phase[pickband], forced_sources[pickband]['psfFlux'],\n",
+    "    ax2.plot(phase[pickband], forced_source[pickband]['psfFlux'],\n",
     "             marker=filter_symbols[band], color=filter_colors[band], linestyle='none',\n",
     "             markersize=7, fillstyle='none')\n",
     "\n",
@@ -443,7 +528,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "mag = forced_sources['psfFlux'].to(u.ABmag).value"
+    "mag = forced_source['psfFlux'].to(u.ABmag).value"
    ]
   },
   {
@@ -456,7 +541,7 @@
     "fig, ax = plt.subplots(1, 1, figsize=(7, 5))\n",
     "\n",
     "for band in filter_names:\n",
-    "    pickband = np.where(forced_sources['band'] == band)\n",
+    "    pickband = np.where(forced_source['band'] == band)\n",
     "    ax.plot(phase[pickband], mag[pickband],\n",
     "            marker=filter_symbols[band], color=filter_colors[band],\n",
     "            label=f\"{band}\", markersize=7, linestyle='none', fillstyle='none')\n",