diff --git a/notebooks/simulated_asteroid_ephem_rsp.ipynb b/notebooks/simulated_asteroid_ephem_rsp.ipynb deleted file mode 100644 index 7c4b0df..0000000 --- a/notebooks/simulated_asteroid_ephem_rsp.ipynb +++ /dev/null @@ -1,657 +0,0 @@ -{ - "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 -}