updates

Author: MelissaGraham

MPC/RFL_Solar_System_MPC_Retrieval.ipynb

@@ -84,7 +84,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 +145,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 +182,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,7 +260,7 @@
    "id": "b57066e8-f27d-4f1b-932f-529f4bb75cb7",
    "metadata": {},
    "source": [
-    "### 2.2 Plot orbit classifications"
+    "### 2.2. Plot orbit classifications"
    ]
   },
   {
@@ -256,7 +278,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 +303,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 +315,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 +338,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 +377,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 +409,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 +425,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 +509,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 +521,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 +541,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 +563,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 +606,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 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": []