diff --git a/notebooks/simulated_asteroid_ephem_rsp.ipynb b/notebooks/simulated_asteroid_ephem_rsp.ipynb new file mode 100644 index 0000000..7c4b0df --- /dev/null +++ b/notebooks/simulated_asteroid_ephem_rsp.ipynb @@ -0,0 +1,657 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Testing DP0.3 ephemerides and coordinates\n", + "By Jamie Robinson" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This notebook retrieves the observations and orbital parameters of an object in DP0.3. We compare the ephemerides and coordinates from DP0.3 to the output from propagating the orbit." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:15.026373Z", + "iopub.status.busy": "2024-06-18T19:01:15.026079Z", + "iopub.status.idle": "2024-06-18T19:01:17.359920Z", + "shell.execute_reply": "2024-06-18T19:01:17.359236Z", + "shell.execute_reply.started": "2024-06-18T19:01:15.026354Z" + } + }, + "outputs": [], + "source": [ + "import pandas as pd\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import matplotlib.gridspec as gridspec\n", + "from astropy import units as u\n", + "from astropy.coordinates import SkyCoord, GCRS\n", + "from sbpy.data import Orbit\n", + "from astropy.time import Time\n", + "from astropy.coordinates import solar_system_ephemeris\n", + "from sbpy.photometry import HG\n", + "from astropy.table import QTable\n", + "\n", + "from lsst.rsp import get_tap_service" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:17.361390Z", + "iopub.status.busy": "2024-06-18T19:01:17.360891Z", + "iopub.status.idle": "2024-06-18T19:01:17.438804Z", + "shell.execute_reply": "2024-06-18T19:01:17.437930Z", + "shell.execute_reply.started": "2024-06-18T19:01:17.361368Z" + } + }, + "outputs": [], + "source": [ + "service = get_tap_service(\"ssotap\")\n", + "assert service is not None" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:17.439847Z", + "iopub.status.busy": "2024-06-18T19:01:17.439631Z", + "iopub.status.idle": "2024-06-18T19:01:17.442630Z", + "shell.execute_reply": "2024-06-18T19:01:17.442049Z", + "shell.execute_reply.started": "2024-06-18T19:01:17.439830Z" + } + }, + "outputs": [], + "source": [ + "# ssObjectId of object to analyse\n", + "ssoid = \"6098332225018\" # good test object" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:17.443664Z", + "iopub.status.busy": "2024-06-18T19:01:17.443428Z", + "iopub.status.idle": "2024-06-18T19:01:18.508689Z", + "shell.execute_reply": "2024-06-18T19:01:18.508031Z", + "shell.execute_reply.started": "2024-06-18T19:01:17.443646Z" + } + }, + "outputs": [], + "source": [ + "query = \"\"\"\n", + "SELECT\n", + " *\n", + "FROM\n", + " dp03_catalogs_10yr.DiaSource as dia\n", + "INNER JOIN\n", + " dp03_catalogs_10yr.SSSource as sss\n", + "ON\n", + " dia.diaSourceId = sss.diaSourceId\n", + "WHERE\n", + " dia.ssObjectId={}\n", + "ORDER by dia.ssObjectId\n", + "\"\"\".format(\n", + " ssoid\n", + ")\n", + "\n", + "df = service.search(query).to_table().to_pandas()\n", + "df" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:18.509807Z", + "iopub.status.busy": "2024-06-18T19:01:18.509560Z", + "iopub.status.idle": "2024-06-18T19:01:18.661531Z", + "shell.execute_reply": "2024-06-18T19:01:18.660953Z", + "shell.execute_reply.started": "2024-06-18T19:01:18.509790Z" + } + }, + "outputs": [], + "source": [ + "results = service.search(\"SELECT * FROM dp03_catalogs_10yr.MPCORB \" \"WHERE ssObjectId={}\".format(ssoid))\n", + "df_orb = results.to_table().to_pandas()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:18.663655Z", + "iopub.status.busy": "2024-06-18T19:01:18.663419Z", + "iopub.status.idle": "2024-06-18T19:01:18.683552Z", + "shell.execute_reply": "2024-06-18T19:01:18.682809Z", + "shell.execute_reply.started": "2024-06-18T19:01:18.663638Z" + } + }, + "outputs": [], + "source": [ + "df_orb" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:18.685016Z", + "iopub.status.busy": "2024-06-18T19:01:18.684704Z", + "iopub.status.idle": "2024-06-18T19:01:18.699123Z", + "shell.execute_reply": "2024-06-18T19:01:18.698473Z", + "shell.execute_reply.started": "2024-06-18T19:01:18.684988Z" + } + }, + "outputs": [], + "source": [ + "# put obs in time order\n", + "df = df.sort_values(\"midPointMjdTai\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Resample the sparse observations by propagating orbit with sbpy & oorb" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:18.700389Z", + "iopub.status.busy": "2024-06-18T19:01:18.700015Z", + "iopub.status.idle": "2024-06-18T19:01:18.715818Z", + "shell.execute_reply": "2024-06-18T19:01:18.715179Z", + "shell.execute_reply.started": "2024-06-18T19:01:18.700367Z" + } + }, + "outputs": [], + "source": [ + "# check orbit\n", + "e = df_orb[\"e\"]\n", + "incl = df_orb[\"incl\"]\n", + "q = df_orb[\"q\"]\n", + "a = q / (1.0 - e)\n", + "Q = a * (1.0 + e)\n", + "print(a, e, incl)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:18.716909Z", + "iopub.status.busy": "2024-06-18T19:01:18.716682Z", + "iopub.status.idle": "2024-06-18T19:01:18.743743Z", + "shell.execute_reply": "2024-06-18T19:01:18.743021Z", + "shell.execute_reply.started": "2024-06-18T19:01:18.716892Z" + } + }, + "outputs": [], + "source": [ + "df_orb[\"a\"] = a\n", + "df_orb[\"Q\"] = Q\n", + "df_orb" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:18.745120Z", + "iopub.status.busy": "2024-06-18T19:01:18.744684Z", + "iopub.status.idle": "2024-06-18T19:01:18.761040Z", + "shell.execute_reply": "2024-06-18T19:01:18.760412Z", + "shell.execute_reply.started": "2024-06-18T19:01:18.745100Z" + } + }, + "outputs": [], + "source": [ + "# Calculate some extra orbital elements\n", + "df_orb[\"P\"] = df_orb[\"a\"] ** (3.0 / 2.0) # orbital period in years\n", + "df_orb[\"n\"] = 360.0 / (df_orb[\"P\"] * 365.25) # mean motion in deg/day\n", + "df_orb[\"M\"] = (\n", + " df_orb[\"n\"] * (df_orb[\"epoch\"] - df_orb[\"tperi\"])\n", + ") % 360 # angles must be in correct range otherwise sbpy/pyoorb freak out\n", + "df_orb[[\"a\", \"P\", \"n\", \"M\"]]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:18.762022Z", + "iopub.status.busy": "2024-06-18T19:01:18.761789Z", + "iopub.status.idle": "2024-06-18T19:01:18.775973Z", + "shell.execute_reply": "2024-06-18T19:01:18.775117Z", + "shell.execute_reply.started": "2024-06-18T19:01:18.762000Z" + } + }, + "outputs": [], + "source": [ + "# rename columns for consistency\n", + "df_orb = df_orb.rename(columns={\"node\": \"Omega\", \"peri\": \"w\"})" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:18.777363Z", + "iopub.status.busy": "2024-06-18T19:01:18.776889Z", + "iopub.status.idle": "2024-06-18T19:01:18.796868Z", + "shell.execute_reply": "2024-06-18T19:01:18.796165Z", + "shell.execute_reply.started": "2024-06-18T19:01:18.777340Z" + } + }, + "outputs": [], + "source": [ + "df_orb[[\"a\", \"e\", \"incl\", \"Omega\", \"w\", \"M\"]]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:18.798419Z", + "iopub.status.busy": "2024-06-18T19:01:18.797782Z", + "iopub.status.idle": "2024-06-18T19:01:18.816304Z", + "shell.execute_reply": "2024-06-18T19:01:18.815461Z", + "shell.execute_reply.started": "2024-06-18T19:01:18.798397Z" + } + }, + "outputs": [], + "source": [ + "# create an sbpy oorb object from dataframe via QTable\n", + "tab = QTable.from_pandas(\n", + " df_orb[[\"a\", \"e\", \"incl\", \"Omega\", \"w\", \"M\"]],\n", + " units={\"a\": u.au, \"incl\": u.deg, \"Omega\": u.deg, \"w\": u.deg, \"M\": u.deg},\n", + ")\n", + "orbit = Orbit.from_table(tab)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:18.817772Z", + "iopub.status.busy": "2024-06-18T19:01:18.817316Z", + "iopub.status.idle": "2024-06-18T19:01:18.833560Z", + "shell.execute_reply": "2024-06-18T19:01:18.832844Z", + "shell.execute_reply.started": "2024-06-18T19:01:18.817749Z" + } + }, + "outputs": [], + "source": [ + "# oorb requires certain extra fields\n", + "orbit[\"epoch\"] = Time(Time(df_orb[\"epoch\"], format=\"mjd\").jd, format=\"jd\")\n", + "orbit[\"targetname\"] = np.array(df_orb[\"ssObjectId\"]).astype(str)\n", + "orbit[\"H\"] = df_orb[\"mpcH\"] * u.mag\n", + "orbit[\"G\"] = df_orb[\"mpcG\"] * u.dimensionless_unscaled" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:18.834791Z", + "iopub.status.busy": "2024-06-18T19:01:18.834470Z", + "iopub.status.idle": "2024-06-18T19:01:18.848589Z", + "shell.execute_reply": "2024-06-18T19:01:18.847950Z", + "shell.execute_reply.started": "2024-06-18T19:01:18.834772Z" + } + }, + "outputs": [], + "source": [ + "orbit" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:18.849539Z", + "iopub.status.busy": "2024-06-18T19:01:18.849320Z", + "iopub.status.idle": "2024-06-18T19:01:18.865983Z", + "shell.execute_reply": "2024-06-18T19:01:18.865340Z", + "shell.execute_reply.started": "2024-06-18T19:01:18.849513Z" + } + }, + "outputs": [], + "source": [ + "# define a set of JD times to propagate the orbital elements to\n", + "N = 1000\n", + "times = Time(\n", + " Time(np.linspace(np.amin(df[\"midPointMjdTai\"]), np.amax(df[\"midPointMjdTai\"]), N), format=\"mjd\").jd,\n", + " format=\"jd\",\n", + ")\n", + "times[0]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:18.867204Z", + "iopub.status.busy": "2024-06-18T19:01:18.866796Z", + "iopub.status.idle": "2024-06-18T19:01:18.880152Z", + "shell.execute_reply": "2024-06-18T19:01:18.879521Z", + "shell.execute_reply.started": "2024-06-18T19:01:18.867185Z" + } + }, + "outputs": [], + "source": [ + "# create an empty dataframe to hold resampled observations\n", + "df_dense = pd.DataFrame()\n", + "df_dense[\"midPointMjdTai\"] = times.mjd" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:18.881123Z", + "iopub.status.busy": "2024-06-18T19:01:18.880917Z", + "iopub.status.idle": "2024-06-18T19:01:51.982325Z", + "shell.execute_reply": "2024-06-18T19:01:51.981730Z", + "shell.execute_reply.started": "2024-06-18T19:01:18.881108Z" + } + }, + "outputs": [], + "source": [ + "# propagate the orbit forward in time.\n", + "# probably a better way to do this but I can't get oo_propagate to work with a Time list right now\n", + "# see: https://github.com/NASA-Planetary-Science/sbpy/issues/341\n", + "\n", + "df_pos = pd.DataFrame() # empty dataframe to hold cartesian coordinates\n", + "\n", + "for i in range(len(times)):\n", + " print(i)\n", + " prop_elem = orbit.oo_propagate(times[i]) # propagate the orbit to the selected time step\n", + " del prop_elem.table[\n", + " \"orbtype\"\n", + " ] # orbtype is added as int, sbpy freaks out so delete the orbtype and then _to_oo works it out\n", + " print(\"propagate\")\n", + " statevec = prop_elem.oo_transform(\"CART\") # transform from orbital elements to cartesian\n", + " print(\"transform\")\n", + "\n", + " # append new cartesian coordinates to the dataframe\n", + " _df_statevec = statevec.table.to_pandas()\n", + " df_pos = pd.concat((df_pos, _df_statevec))\n", + "\n", + "df_pos.reset_index(drop=True, inplace=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:51.983360Z", + "iopub.status.busy": "2024-06-18T19:01:51.983151Z", + "iopub.status.idle": "2024-06-18T19:01:52.000157Z", + "shell.execute_reply": "2024-06-18T19:01:51.999473Z", + "shell.execute_reply.started": "2024-06-18T19:01:51.983345Z" + } + }, + "outputs": [], + "source": [ + "df_pos" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Use astropy coordinates to transform between all the coordinate systems" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:52.001741Z", + "iopub.status.busy": "2024-06-18T19:01:52.001163Z", + "iopub.status.idle": "2024-06-18T19:01:52.010357Z", + "shell.execute_reply": "2024-06-18T19:01:52.009551Z", + "shell.execute_reply.started": "2024-06-18T19:01:52.001719Z" + } + }, + "outputs": [], + "source": [ + "# define heliocentric cartesian coordinates\n", + "c_xyz_hel = SkyCoord(\n", + " x=np.array(df_pos[\"x\"]),\n", + " y=np.array(df_pos[\"y\"]),\n", + " z=np.array(df_pos[\"z\"]),\n", + " unit=\"AU\",\n", + " representation_type=\"cartesian\",\n", + " frame=\"heliocentrictrueecliptic\",\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:52.013233Z", + "iopub.status.busy": "2024-06-18T19:01:52.012889Z", + "iopub.status.idle": "2024-06-18T19:01:52.025365Z", + "shell.execute_reply": "2024-06-18T19:01:52.024660Z", + "shell.execute_reply.started": "2024-06-18T19:01:52.013210Z" + } + }, + "outputs": [], + "source": [ + "# transform to heliocentric ecliptic coords\n", + "c_ecl_hel = c_xyz_hel.copy()\n", + "c_ecl_hel.representation_type = \"spherical\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:52.026475Z", + "iopub.status.busy": "2024-06-18T19:01:52.026241Z", + "iopub.status.idle": "2024-06-18T19:01:52.095939Z", + "shell.execute_reply": "2024-06-18T19:01:52.095204Z", + "shell.execute_reply.started": "2024-06-18T19:01:52.026457Z" + } + }, + "outputs": [], + "source": [ + "# transform to geocentric equatorial coords (times required to calculate Earth position)\n", + "with solar_system_ephemeris.set(\"jpl\"):\n", + " c_eq_geo = c_xyz_hel.transform_to(GCRS(obstime=times))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:52.096956Z", + "iopub.status.busy": "2024-06-18T19:01:52.096748Z", + "iopub.status.idle": "2024-06-18T19:01:52.100730Z", + "shell.execute_reply": "2024-06-18T19:01:52.099911Z", + "shell.execute_reply.started": "2024-06-18T19:01:52.096940Z" + } + }, + "outputs": [], + "source": [ + "# transform to geocentric cartesian coords\n", + "c_xyz_geo = c_eq_geo.copy()\n", + "c_xyz_geo.representation_type = \"cartesian\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:52.101667Z", + "iopub.status.busy": "2024-06-18T19:01:52.101450Z", + "iopub.status.idle": "2024-06-18T19:01:52.118900Z", + "shell.execute_reply": "2024-06-18T19:01:52.118257Z", + "shell.execute_reply.started": "2024-06-18T19:01:52.101652Z" + } + }, + "outputs": [], + "source": [ + "# transform from geo equatorial (ra, dec) to geo ecliptic (lon, lat)\n", + "c_ecl_geo = c_eq_geo.transform_to(\"geocentrictrueecliptic\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:52.120101Z", + "iopub.status.busy": "2024-06-18T19:01:52.119878Z", + "iopub.status.idle": "2024-06-18T19:01:52.364967Z", + "shell.execute_reply": "2024-06-18T19:01:52.364185Z", + "shell.execute_reply.started": "2024-06-18T19:01:52.120086Z" + } + }, + "outputs": [], + "source": [ + "# plot the propagated cartesian positions against the database values\n", + "\n", + "x_plot = \"midPointMjdTai\"\n", + "y_plot1 = \"heliocentricX\"\n", + "y_plot2 = \"heliocentricY\"\n", + "y_plot3 = \"heliocentricZ\"\n", + "df_plot = df\n", + "\n", + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "x = ax1.scatter(df_plot[x_plot], df_plot[y_plot1], label=y_plot1)\n", + "x = ax1.scatter(df_plot[x_plot], df_plot[y_plot2], label=y_plot2)\n", + "x = ax1.scatter(df_plot[x_plot], df_plot[y_plot3], label=y_plot3)\n", + "\n", + "ax1.plot(times.mjd, c_xyz_hel.x)\n", + "ax1.plot(times.mjd, c_xyz_hel.y)\n", + "ax1.plot(times.mjd, c_xyz_hel.z)\n", + "\n", + "ax1.set_xlabel(x_plot)\n", + "ax1.set_ylabel(\"distance\")\n", + "ax1.legend()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "There are some deviations of x and y positions, probably due to slightly different reference frames and methods of propagating orbits.\n", + "\n", + "There is something wrong with database z positions, will be fixed soon!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "execution": { + "iopub.execute_input": "2024-06-18T19:01:52.366083Z", + "iopub.status.busy": "2024-06-18T19:01:52.365860Z", + "iopub.status.idle": "2024-06-18T19:01:52.558412Z", + "shell.execute_reply": "2024-06-18T19:01:52.557808Z", + "shell.execute_reply.started": "2024-06-18T19:01:52.366066Z" + } + }, + "outputs": [], + "source": [ + "# the ecliptic coordinates look good!\n", + "\n", + "x_plot = \"midPointMjdTai\"\n", + "y_plot1 = \"eclipticLambda\"\n", + "y_plot2 = \"eclipticBeta\"\n", + "df_plot = df\n", + "\n", + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "x = ax1.scatter(df_plot[x_plot], df_plot[y_plot1], label=y_plot1)\n", + "x = ax1.scatter(df_plot[x_plot], df_plot[y_plot2], label=y_plot2)\n", + "\n", + "ax1.plot(times.mjd, c_ecl_geo.lon.degree)\n", + "ax1.plot(times.mjd, c_ecl_geo.lat.degree)\n", + "\n", + "ax1.set_xlabel(x_plot)\n", + "ax1.set_ylabel(\"angle\")\n", + "ax1.legend()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "lsst-sbpy", + "language": "python", + "name": "lsst-sbpy" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.7" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}