Merge pull request #97 from lsst/tickets/SP-2540

MelissaGraham · web-flow · commit 189e886e9812 · 2025-11-25T18:34:53.000-08:00
tickets/SP-2540: add photometry evaluation
diff --git a/MPC/RFL_Solar_System_MPC_Retrieval.ipynb b/MPC/RFL_Solar_System_MPC_Retrieval.ipynb
@@ -18,12 +18,13 @@
     "\n",
     "</div>\n",
     "\n",
-    "For the Rubin Science Platform at data.lsst.cloud. <br>\n",
-    "For Solar System object discoveries released as part of Rubin First Look. <br>\n",
-    "Container Size: large <br>\n",
-    "LSST Science Pipelines version: r29.2.0 <br>\n",
-    "Last verified to run: 2025-09-02 <br>\n",
-    "Repository: <a href=\"https://github.com/lsst/tutorial-notebooks\">github.com/lsst/tutorial-notebooks</a> <br>\n"
+    "For the Rubin Science Platform at data.lsst.cloud.\\\n",
+    "For Solar System object discoveries released as part of Rubin First Look.\\\n",
+    "Container Size: large\\\n",
+    "LSST Science Pipelines version: r29.2.0\\\n",
+    "Last verified to run: 2025-11-25\\\n",
+    "Repository: [github.com/lsst/tutorial-notebooks](https://github.com/lsst/tutorial-notebooks)\\\n",
+    "DOI: [10.11578/rubin/dc.20250909.20](https://doi.org/10.11578/rubin/dc.20250909.20)"
    ]
   },
   {
@@ -35,14 +36,14 @@
     "\n",
     "**LSST data products:** Small body observations made by Rubin, queried from the Minor Planet Center (MPC), and orbits computed by the MPC from those Rubin observations are used in the notebook. No direct LSST data products are used in this notebook.\n",
     "\n",
-    "**Packages:** `get_orb` and `get_obs` <a href=\"https://www.minorplanetcenter.net/mpcops/documentation/\">APIs from the Minor Planet Center</a>.\n",
+    "**Packages:** `get_orb` and `get_obs` [APIs from the Minor Planet Center](https://www.minorplanetcenter.net/mpcops/documentation/).\n",
     "\n",
     "**Credit:**\n",
     "Originally developed by Sarah Greenstreet and Mario Juric from the Rubin Community Science team, Solar System Processing (SSP) team and Solar System Commissioning science unit. Please consider acknowledging them if this notebook is used for the preparation of journal articles, software releases, or other notebooks.\n",
     "\n",
     "**Get Support:**\n",
-    "Everyone is encouraged to ask questions or raise issues in the <a href=\"https://community.lsst.org/c/support\">Support Category</a> \n",
-    "of the Rubin Community Forum. Rubin staff will respond to all questions posted there."
+    "Everyone is encouraged to ask questions or raise issues in the [Support Category](https://community.lsst.org/c/support) of the Rubin Community Forum.\n",
+    "Rubin staff will respond to all questions posted there."
    ]
   },
   {
@@ -66,11 +67,8 @@
    "source": [
     "### 1.1. Import packages\n",
     "\n",
-    "Import `numpy`, a fundamental package for scientific computing with arrays in Python\n",
-    "(<a href=\"https://numpy.org\">numpy.org</a>), and\n",
-    "`matplotlib`, a comprehensive library for data visualization\n",
-    "(<a href=\"https://matplotlib.org/\">matplotlib.org</a>; \n",
-    "<a href=\"https://matplotlib.org/stable/gallery/index.html\">matplotlib gallery</a>)."
+    "Import `numpy`, a fundamental package for scientific computing with arrays in Python ([numpy.org](https://numpy.org)), and\n",
+    "`matplotlib`, a comprehensive library for data visualization ([matplotlib.org](https://matplotlib.org/); [matplotlib gallery](https://matplotlib.org/stable/gallery/index.html))."
    ]
   },
   {
@@ -84,7 +82,9 @@
     "import numpy as np\n",
     "import pandas as pd\n",
     "import matplotlib.pyplot as plt\n",
-    "from astropy.time import Time"
+    "from astropy.time import Time\n",
+    "from lsst.utils.plotting import (get_multiband_plot_colors,\n",
+    "                                 get_multiband_plot_symbols)"
    ]
   },
   {
@@ -143,6 +143,26 @@
     "colors = prop_cycle.by_key()['color']"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "fd91df61-fb09-48e4-88a1-bd9fa6786aec",
+   "metadata": {},
+   "source": [
+    "Define the colors and symbols to use in plots to represent the six LSST filters, $ugrizy$."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f9d6274f-1d2f-428e-b84b-2d3a8fb79fb9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "filter_names = ['u', 'g', 'r', 'i', 'z', 'y']\n",
+    "filter_colors = get_multiband_plot_colors()\n",
+    "filter_symbols = get_multiband_plot_symbols()"
+   ]
+  },
   {
    "cell_type": "markdown",
    "id": "6d5dae6e-e177-497a-b61f-b078809ae987",
@@ -160,7 +180,7 @@
    "id": "dd380a4c-858b-42cf-bd2c-43fb9311c96a",
    "metadata": {},
    "source": [
-    "### 2.1 Query the Minor Planet Center\n",
+    "### 2.1. Query the Minor Planet Center\n",
     "\n",
     "Use the Minor Planet Center's `get_orb` API to retrieve the orbit information for all 2,098 objects in the Rubin First Look Solar System discoveries and store the list of retrieved dictionaries in a DataFrame called `discovery_orbits`. Return fields for the object's provisional designation, absolute magnitude, semimajor axis, eccentricity, inclination, argument of perihelion, longitude of ascending node, time of perihelion passage, the total number of observations used in the orbit fit, and the orbit classification information. The request may take roughly 10 minutes or more to complete.\n",
     "\n",
@@ -238,15 +258,16 @@
    "id": "b57066e8-f27d-4f1b-932f-529f4bb75cb7",
    "metadata": {},
    "source": [
-    "### 2.2 Plot orbit classifications"
+    "### 2.2. Plot orbit classifications"
    ]
   },
   {
    "cell_type": "markdown",
    "id": "a928a607-9977-4a09-b689-9efeb49b8437",
    "metadata": {},
    "source": [
-    "Plot the semimajor axes, eccentricities, and inclinations of the 2,098 Rubin First Look Solar System object discoveries. Use the orbit types defined by the MPC (https://minorplanetcenter.net/mpcops/documentation/orbit-types/) to distinguish each orbital classification in the plot.\n"
+    "Plot the semimajor axes, eccentricities, and inclinations of the 2,098 Rubin First Look Solar System object discoveries.\n",
+    "Use the [orbit types defined by the MPC](https://minorplanetcenter.net/mpcops/documentation/orbit-types/) to distinguish each orbital classification in the plot."
    ]
   },
   {
@@ -256,7 +277,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "a, e, inc = discovery_orbits[\"semimajor_axis\"], discovery_orbits[\"eccentricity\"], discovery_orbits[\"inclination\"]\n",
+    "a = discovery_orbits[\"semimajor_axis\"]\n",
+    "e = discovery_orbits[\"eccentricity\"]\n",
+    "inc = discovery_orbits[\"inclination\"]\n",
     "inner_orbits = discovery_orbits[a < 7]\n",
     "outer_orbits = discovery_orbits[a >= 7]\n",
     "\n",
@@ -279,7 +302,7 @@
     "    99: \"Other\"\n",
     "}\n",
     "\n",
-    "fig, axes = plt.subplots(1, 2, figsize=(8,4))\n",
+    "fig, axes = plt.subplots(1, 2, figsize=(8, 4))\n",
     "plt.subplots_adjust(wspace=0.25)\n",
     "\n",
     "for ax, oo in zip(axes, (inner_orbits, outer_orbits)):\n",
@@ -291,11 +314,12 @@
     "            s = 10\n",
     "            marker, markers = markers[0], markers[1:]\n",
     "        label = orbit_type_map[otype]\n",
-    "        ax.scatter(orb[\"semimajor_axis\"], orb[\"eccentricity\"], s=s, marker=marker, label=f\"{label} ({len(orb)})\")\n",
+    "        ax.scatter(orb[\"semimajor_axis\"], orb[\"eccentricity\"], s=s,\n",
+    "                   marker=marker, label=f\"{label} ({len(orb)})\")\n",
     "\n",
     "    for handle in ax.legend().legend_handles:\n",
     "        handle.set_sizes([10])\n",
-    "    ax.legend(loc = 'upper right', bbox_to_anchor=(0.92, 1), fontsize=6)\n",
+    "    ax.legend(loc='upper right', bbox_to_anchor=(0.92, 1), fontsize=6)\n",
     "    ax.set_xlabel(\"semimajor axis (au)\")\n",
     "    ax.set_ylabel(\"eccentricity\")\n",
     "\n",
@@ -313,11 +337,12 @@
     "            s = 10\n",
     "            marker, markers = markers[0], markers[1:]\n",
     "        label = orbit_type_map[otype]\n",
-    "        ax.scatter(orb[\"semimajor_axis\"], orb[\"inclination\"], s=s, marker=marker, label=f\"{label} ({len(orb)})\")\n",
+    "        ax.scatter(orb[\"semimajor_axis\"], orb[\"inclination\"], s=s,\n",
+    "                   marker=marker, label=f\"{label} ({len(orb)})\")\n",
     "\n",
     "    for handle in ax.legend().legend_handles:\n",
     "        handle.set_sizes([10])\n",
-    "    ax.legend(loc = 'upper right', fontsize=6)\n",
+    "    ax.legend(loc='upper right', fontsize=6)\n",
     "    ax.set_xlabel(\"semimajor axis (au)\")\n",
     "    ax.set_ylabel(\"inclination (deg)\")"
    ]
@@ -351,7 +376,7 @@
    "id": "2edbdaba-551d-4b72-a346-16227a9c45fb",
    "metadata": {},
    "source": [
-    "### 3.1 Query the Minor Planet Center\n",
+    "### 3.1. Query the Minor Planet Center\n",
     "\n",
     "Use the Minor Planet Center's `get_obs` API to retrieve all observations for the 2,098 objects in the Rubin First Look Solar System discoveries and store them in a pandas DataFrame called `discovery_observations`, reducing the number of columns to those needed and converting selected strings to floats using numpy's `to_numeric` module. The request may take up to 40 minutes or more to complete."
    ]
@@ -383,7 +408,7 @@
     "\n",
     "for i, obj in enumerate(targets[1:]):\n",
     "    response = requests.get(\"https://data.minorplanetcenter.net/api/get-obs\",\n",
-    "                            json={\"desigs\": [obj], \"output_format\":[\"ADES_DF\"]})\n",
+    "                            json={\"desigs\": [obj], \"output_format\": [\"ADES_DF\"]})\n",
     "    if response.ok:\n",
     "        ades_df = pd.DataFrame(response.json()[0]['ADES_DF'])\n",
     "    else:\n",
@@ -399,11 +424,12 @@
     "                                        'poscov13', 'poscov22', 'poscov23',\n",
     "                                        'poscov33', 'precdec', 'precra',\n",
     "                                        'prectime', 'prog', 'rastar', 'rcv',\n",
-    "                                        'ref', 'remarks', 'rmscorr','rmsdelay',\n",
+    "                                        'ref', 'remarks', 'rmscorr', 'rmsdelay',\n",
     "                                        'rmsdist', 'rmsdoppler', 'rmsfit',\n",
     "                                        'rmspa', 'rmstime', 'seeing', 'subfmt',\n",
     "                                        'subfrm', 'sys', 'trkid', 'trx', 'unctime'], axis=1)\n",
-    "    discovery_observations = pd.concat([discovery_observations, ingest_observations], ignore_index=True)\n",
+    "    discovery_observations = pd.concat([discovery_observations, ingest_observations],\n",
+    "                                       ignore_index=True)\n",
     "    discovery_observations['dec'] = pd.to_numeric(discovery_observations['dec'])\n",
     "    discovery_observations['ra'] = pd.to_numeric(discovery_observations['ra'])\n",
     "    discovery_observations['rmsdec'] = pd.to_numeric(discovery_observations['rmsdec'])\n",
@@ -482,7 +508,7 @@
    "id": "2a885ec6-9861-435b-9025-998005fa6735",
    "metadata": {},
    "source": [
-    "### 3.2 Examine the magnitude variation for an example discovery\n",
+    "### 3.2. Examine the magnitude variation for an example discovery\n",
     "\n",
     "Sort `discovery_orbits` by number of observations in descending order to find the objects with the largest number of observations. Print the 10 objects with the highest number of observations."
    ]
@@ -494,7 +520,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "discovery_orbits_observations_sorted_desc = discovery_orbits.sort_values(by='nobs_total', ascending=False)\n",
+    "discovery_orbits_observations_sorted_desc = discovery_orbits.sort_values(by='nobs_total',\n",
+    "                                                                         ascending=False)\n",
     "discovery_orbits_observations_sorted_desc.head(10)"
    ]
   },
@@ -513,7 +540,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "select = np.isin(Rubin_discovery_observations['provid'], '2025 MX21', invert=False)\n",
+    "object_name = '2025 MX21'\n",
+    "select = np.isin(Rubin_discovery_observations['provid'], object_name, invert=False)\n",
     "selected_Rubin_discovery_observations = Rubin_discovery_observations[select].copy()\n",
     "len(selected_Rubin_discovery_observations)"
    ]
@@ -534,26 +562,37 @@
    "outputs": [],
    "source": [
     "MJDmin, MJDmax = 60797, 60797.2\n",
-    "MJD, mag = selected_Rubin_discovery_observations['epoch'], selected_Rubin_discovery_observations['mag']\n",
+    "MJD = selected_Rubin_discovery_observations['epoch']\n",
+    "mag = selected_Rubin_discovery_observations['mag']\n",
     "full_obs = selected_Rubin_discovery_observations\n",
     "single_night_obs = selected_Rubin_discovery_observations[(MJDmin < MJD) & (MJD < MJDmax)]\n",
     "\n",
     "fig, axes = plt.subplots(1, 2, figsize=(14, 4))\n",
     "\n",
+    "axes[0].axvline(MJDmin + 0.1, lw=10, alpha=0.2, color='grey')\n",
+    "a = 0\n",
     "for ax, mags in zip(axes, (full_obs, single_night_obs)):\n",
-    "    markers = \"s^*+xDh\"*2\n",
-    "    for label, mag in mags.groupby('band'):\n",
-    "        s = 10\n",
-    "        marker, markers = markers[0], markers[1:]\n",
-    "        ax.scatter(mag['epoch'], mag['mag'], s=s, marker=marker, label=f\"{label} ({len(mag)})\")\n",
+    "    for filt in ['g', 'r', 'i']:\n",
+    "        select_filt = np.isin(mags['band'], filt, invert=False)\n",
+    "        tx = np.where(select_filt == 1)[0]\n",
+    "        ax.scatter(mags['epoch'][select_filt], mags['mag'][select_filt],\n",
+    "                   s=10, marker=filter_symbols[filt],\n",
+    "                   color=filter_colors[filt],\n",
+    "                   label=f\"{filt} ({len(tx)})\")\n",
+    "        if a == 1:\n",
+    "            ax.errorbar(mags['epoch'][select_filt], mags['mag'][select_filt],\n",
+    "                        yerr=mags['rmsmag'][select_filt], fmt='None',\n",
+    "                        ecolor=filter_colors[filt])\n",
+    "        del select_filt, tx\n",
+    "    a += 1\n",
     "\n",
     "    for handle in ax.legend().legend_handles:\n",
     "        handle.set_sizes([10])\n",
     "    ax.invert_yaxis()\n",
     "    ax.set_xlabel('MJD')\n",
     "    ax.set_ylabel('Magnitude')\n",
     "    ax.minorticks_on()\n",
-    "fig.suptitle(f\"Rubin First Look Solar System Object Discovery {selected_Rubin_discovery_observations['provid'].values[0]} Magnitude Variation\")\n",
+    "fig.suptitle(f\"Rubin First Look Solar System Object Discovery {object_name} Magnitude Variation\")\n",
     "\n",
     "plt.show()"
    ]
@@ -566,10 +605,106 @@
     "> **Figure 2:** Magnitude over time for the selected Rubin First Look Solar System object discovery both for the full observation period (left) and on a single night (right) to see the magnitude variation. Observations were taken in $g$-, $r$-, and $i$-band."
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "8af622a5-3d13-4a27-9c99-393360927989",
+   "metadata": {},
+   "source": [
+    "### 3.3. Evaluate photometric quality\n",
+    "\n",
+    "In a full Rubin data release, the photometry for moving objects in the `DiaSource` catalog comes with flag parameters to indicate potential issues with the photometric measurement or calibration (e.g., contamination from cosmic ray flux, sources at the edge of a detector), and the difference- and direct-images are available for visual checks on the source being measured.\n",
+    "Observational metadata such as seeing, sky background, and limiting metadata are also available.\n",
+    "These aspects of a full Rubin data release enable users to identify and reject potential poor-quality measurements from an analysis.\n",
+    "\n",
+    "However, for the RFL photometry released to the MPC, the only indicator of photometric quality available is the `rmsmag` column."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "590e0d1f-d15b-4944-8a36-76d1091fdd73",
+   "metadata": {},
+   "source": [
+    "Show the distribution of all `rmsmag` values in order to approximate the boundary between the distribution of \"acceptable\" or \"normal\" measurements, and the tail of \"outlier\" measurements."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "738421d1-3180-4a67-9bf5-4537e083294b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig = plt.figure(figsize=(8, 3))\n",
+    "plt.axvline(0.25, ls='dashed', alpha=0.7, color='black')\n",
+    "for filt in ['g', 'r', 'i']:\n",
+    "    select_filt = np.isin(Rubin_discovery_observations['band'],\n",
+    "                          filt, invert=False)\n",
+    "    plt.hist(Rubin_discovery_observations['rmsmag'][select_filt], 80,\n",
+    "             histtype='step', log=True, color=filter_colors[filt], label=filt)\n",
+    "plt.xlabel('rms mag')\n",
+    "plt.ylabel('log(N)')\n",
+    "plt.legend(loc='upper right')\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "4a447a05-0d3a-435a-a28a-2b1f45865444",
+   "metadata": {},
+   "source": [
+    "> **Figure 3:** The log of the number of photometric measurements with a given `rmsmag` value. The dashed grey line is drawn at $0.25$ mag, the approximate (by-eye) division between the distribution of `rmsmag` values and the tail of \"outliers\"."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "dd3811a9-a6cb-4649-aba1-e4cc060de302",
+   "metadata": {},
+   "source": [
+    "Show scatter plots of the `rmsmag` values as a function of `mag` values, per filter."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "94881aa2-32d8-4647-92c7-e6bbde01209b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, axes = plt.subplots(1, 3, figsize=(10, 5))\n",
+    "for f, filt in enumerate(['g', 'r', 'i']):\n",
+    "    select = np.isin(Rubin_discovery_observations['band'], filt, invert=False)\n",
+    "    axes[f].axhline(0.25, alpha=0.5, color='grey')\n",
+    "    axes[f].plot(Rubin_discovery_observations['mag'][select],\n",
+    "                 Rubin_discovery_observations['rmsmag'][select],\n",
+    "                 'o', ms=2, mew=0, alpha=0.7, color=filter_colors[filt])\n",
+    "    axes[f].set_xlabel('mag')\n",
+    "    axes[f].set_title(filt)\n",
+    "axes[0].set_ylabel('rms mag')\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "1d196cd2-5a3a-4e66-8308-a0cd1f23106c",
+   "metadata": {},
+   "source": [
+    "> **Figure 4:** The `rmsmag` as a function of `mag` for the $g$-, $r$-, and $i$-band photometric measurements. The main locus of points in each plot represents \"normal\" photometric measurements, and the outliers are clearly visible. The horizontal lines drawn at `rmsmag` $=0.25$ mag represent about the upper bound on the \"inliers\"."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "07571984-17db-4c2d-b297-96d79f088909",
+   "metadata": {},
+   "source": [
+    "The above demonstrates a simple, by-eye evaluation that `rmsmag` $> 0.25$ mag indicates an \"outlier\" photometric meaurement for the RFL data sent to the MPC.\n",
+    "\n",
+    "Inclusion of these data points in any analysis should be done with caution."
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "9eb55b51-a84b-4404-bea1-ff593692aa69",
+   "id": "ab073769-fba2-4821-82c9-141dca24a75d",
    "metadata": {},
    "outputs": [],
    "source": []
