142 colour functions (#154)

* add x,y,z pos of observation

* include topocentric coords

* add tperi to mpcorb query

* add ecliptic coordinates to query

* include MPCORB fullDesignation

* update fullDesignation unit test

* fix fit covariance matrix

* add more tests

* work with quantity objects

* no units test

* update PhaseCurve parameter names

* fix InitModelDict

* initial apparition finder

* get obs as df func

* initial commit for colour code

* initial colour function

* update colour func

* add colour demo nb

* colour test data

* add unit colour unit tests

* [pre-commit.ci lite] apply automatic fixes

* demo nb

* add description to doc string

---------

Co-authored-by: pre-commit-ci-lite[bot] <117423508+pre-commit-ci-lite[bot]@users.noreply.github.com>

notebooks/colour_functions_demo.ipynb

@@ -0,0 +1,391 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d1ef57c2",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib notebook"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "61f50e0d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib inline"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "82537a61",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from adler.dataclasses.AdlerPlanetoid import AdlerPlanetoid\n",
+    "from adler.dataclasses.AdlerData import AdlerData\n",
+    "from adler.science.PhaseCurve import PhaseCurve\n",
+    "from adler.science.Colour import col_obs_ref\n",
+    "from adler.utilities.plotting_utilities import plot_errorbar\n",
+    "from adler.utilities.science_utilities import apparition_gap_finder, get_df_obs_filt\n",
+    "\n",
+    "import numpy as np\n",
+    "import pandas as pd\n",
+    "import matplotlib.pyplot as plt\n",
+    "import matplotlib.gridspec as gridspec\n",
+    "import astropy.units as u\n",
+    "from scipy.interpolate import splrep, BSpline, CubicSpline\n",
+    "from astropy.modeling.fitting import SLSQPLSQFitter\n",
+    "from astropy.time import Time\n",
+    "from astroquery.jplhorizons import Horizons"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b9f6dcda",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# ssObjectId of object to analyse\n",
+    "ssoid = \"6098332225018\"  # good MBA test object"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c6f705c1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# retrieve the object data\n",
+    "# fname = \"../tests/data/testing_database.db\"\n",
+    "fname = \"/Users/jrobinson/lsst-adler/notebooks/gen_test_data/adler_demo_testing_database.db\"\n",
+    "planetoid = AdlerPlanetoid.construct_from_SQL(ssoid, sql_filename=fname)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "21e0b93e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# check orbit parameters\n",
+    "e = planetoid.MPCORB.e\n",
+    "incl = planetoid.MPCORB.incl\n",
+    "q = planetoid.MPCORB.q\n",
+    "a = q / (1.0 - e)\n",
+    "Q = a * (1.0 + e)\n",
+    "print(a, e, incl)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ee4064e7",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# plot and fit the phase curves in g and r filters\n",
+    "\n",
+    "fig = plt.figure()\n",
+    "gs = gridspec.GridSpec(1, 1)\n",
+    "ax1 = plt.subplot(gs[0, 0])\n",
+    "\n",
+    "adler_cols = AdlerData(ssoid, planetoid.filter_list)\n",
+    "\n",
+    "for filt in [\"g\", \"r\"]:\n",
+    "    # get observations and define a phase curve model to fit\n",
+    "\n",
+    "    sso = planetoid.SSObject_in_filter(filt)\n",
+    "    obs = planetoid.observations_in_filter(filt)\n",
+    "\n",
+    "    H = sso.H\n",
+    "    G12 = sso.G12\n",
+    "    pc = PhaseCurve(H=H * u.mag, phase_parameter_1=G12, model_name=\"HG12_Pen16\")\n",
+    "\n",
+    "    alpha = np.linspace(0, np.amax(obs.phaseAngle)) * u.deg\n",
+    "\n",
+    "    pc_fit = pc.FitModel(\n",
+    "        np.array(getattr(obs, \"phaseAngle\")) * u.deg,\n",
+    "        np.array(getattr(obs, \"reduced_mag\")) * u.mag,\n",
+    "    )\n",
+    "    pc = pc.InitModelSbpy(pc_fit)\n",
+    "    red_mag = pc.ReducedMag(alpha)\n",
+    "\n",
+    "    adler_cols.populate_phase_parameters(filt, **pc.ReturnModelDict())\n",
+    "\n",
+    "    # add this phase curve to the figure\n",
+    "    fig = plot_errorbar(planetoid, filt_list=[filt], fig=fig)\n",
+    "    ax1 = fig.axes[0]\n",
+    "    ax1.plot(\n",
+    "        alpha.value,\n",
+    "        pc.ReducedMag(alpha).value,\n",
+    "        label=\"{}: H={:.2f}, G12={:.2f}\".format(filt, pc.H, pc.phase_parameter_1),\n",
+    "    )\n",
+    "    ax1.legend()\n",
+    "\n",
+    "ax1.invert_yaxis()\n",
+    "\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d0df3c51",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# inspect the r filter phase curve model\n",
+    "adler_cols.get_phase_parameters_in_filter(\"r\", \"HG12_Pen16\").__dict__"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c46fb9eb",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# inspect the g filter phase curve model\n",
+    "adler_cols.get_phase_parameters_in_filter(\"g\", \"HG12_Pen16\").__dict__"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "73871e56",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# get the boundary times for each apparation of the object in the survey\n",
+    "df_obs_all = pd.DataFrame()\n",
+    "for filt in [\"r\", \"g\"]:\n",
+    "    obs = planetoid.observations_in_filter(filt)\n",
+    "    _df_obs = pd.DataFrame(obs.__dict__)\n",
+    "    df_obs_all = pd.concat([df_obs_all, _df_obs])\n",
+    "df_obs_all = df_obs_all.sort_values(\"midPointMjdTai\")\n",
+    "\n",
+    "t_app = apparition_gap_finder(np.array(df_obs_all[\"midPointMjdTai\"]))\n",
+    "print(t_app)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "63980010",
+   "metadata": {},
+   "source": [
+    "Now we can inpsect how the colour of the object varies (or not) as a function of time. The adler function `col_obs_ref` will compare the latest observation in a given filter with observations in another filter. By setting parameter `N_ref` one can set how many past obsevrations to use when calculating the latest colour."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "bd0e0185",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Here we simulate observations coming night-by-night and the calculation of a g-r colour for the object\n",
+    "\n",
+    "# define colour function parameters\n",
+    "\n",
+    "# set number of reference observations to use for colour estimate\n",
+    "# N_ref = 5\n",
+    "N_ref = 3\n",
+    "# N_ref = 1\n",
+    "# N_ref = None\n",
+    "\n",
+    "# observation and filter field names\n",
+    "x_plot = \"midPointMjdTai\"\n",
+    "y_plot = \"AbsMag\"\n",
+    "y_plot = \"reduced_mag\"\n",
+    "yerr_plot = \"magErr\"\n",
+    "filt_obs = \"g\"  # observation filter\n",
+    "filt_ref = \"r\"  # reference filter (we are calculating a filt_obs - filt_ref colour)\n",
+    "\n",
+    "# define colour field names\n",
+    "colour = \"{}-{}\".format(filt_obs, filt_ref)\n",
+    "colErr = \"{}-{}Err\".format(filt_obs, filt_ref)\n",
+    "delta_t_col = \"delta_t_{}\".format(colour)\n",
+    "y_ref_col = \"{}_{}\".format(y_plot, filt_ref)\n",
+    "x1_ref_col = \"{}1_{}\".format(x_plot, filt_ref)\n",
+    "x2_ref_col = \"{}2_{}\".format(x_plot, filt_ref)\n",
+    "\n",
+    "fig = plt.figure()\n",
+    "gs = gridspec.GridSpec(1, 1)\n",
+    "ax1 = plt.subplot(gs[0, 0])\n",
+    "\n",
+    "col_dict_list = []\n",
+    "for app_i in range(len(t_app) - 1):\n",
+    "    # consider only one apparition\n",
+    "    if app_i != 3:\n",
+    "        continue\n",
+    "\n",
+    "    time_min = t_app[app_i]\n",
+    "    time_max = t_app[app_i + 1]\n",
+    "    print(time_min, time_max)\n",
+    "\n",
+    "    _df_obs_all = df_obs_all[\n",
+    "        (df_obs_all[\"midPointMjdTai\"] >= time_min) & (df_obs_all[\"midPointMjdTai\"] < time_max)\n",
+    "    ]\n",
+    "    print(app_i, len(_df_obs_all))\n",
+    "    _time_max = np.amax(_df_obs_all[\"midPointMjdTai\"])\n",
+    "\n",
+    "    # get the phase curve model and observations for each filter\n",
+    "\n",
+    "    # get the stored AdlerData parameters for the observation filter\n",
+    "    ad_g = adler_cols.get_phase_parameters_in_filter(filt_obs, \"HG12_Pen16\")\n",
+    "    pc_g = PhaseCurve().InitModelDict(ad_g.__dict__)  # make the PhaseCurve object from AdlerData\n",
+    "    # get the phase curve model for the reference filter\n",
+    "    ad_r = adler_cols.get_phase_parameters_in_filter(filt_ref, \"HG12_Pen16\")\n",
+    "    pc_r = PhaseCurve().InitModelDict(ad_r.__dict__)\n",
+    "    # get the observations in both filters\n",
+    "    df_obs = get_df_obs_filt(\n",
+    "        planetoid, filt_obs, x1=time_min, x2=_time_max, col_list=[y_plot, yerr_plot], pc_model=pc_g\n",
+    "    )\n",
+    "    df_obs_ref = get_df_obs_filt(\n",
+    "        planetoid, filt_ref, x1=time_min, x2=_time_max, col_list=[y_plot, yerr_plot], pc_model=pc_r\n",
+    "    )\n",
+    "\n",
+    "    ax1.errorbar(df_obs[x_plot], df_obs[y_plot], df_obs[yerr_plot], fmt=\"o\", label=filt_obs)\n",
+    "    ax1.errorbar(df_obs_ref[x_plot], df_obs_ref[y_plot], df_obs_ref[yerr_plot], fmt=\"o\", label=filt_ref)\n",
+    "\n",
+    "    # simulate stepping through each filt_obs observation\n",
+    "    x1 = time_min\n",
+    "    for xi in range(len(df_obs)):\n",
+    "        x2 = df_obs.iloc[xi][x_plot]\n",
+    "        print(xi, x1, x2)\n",
+    "\n",
+    "        # do the colour finding function here\n",
+    "        col_dict = col_obs_ref(\n",
+    "            planetoid,\n",
+    "            adler_cols,\n",
+    "            filt_obs=filt_obs,\n",
+    "            filt_ref=filt_ref,\n",
+    "            N_ref=N_ref,\n",
+    "            x_col=x_plot,\n",
+    "            y_col=y_plot,\n",
+    "            yerr_col=yerr_plot,\n",
+    "            x1=x1,\n",
+    "            x2=x2,\n",
+    "        )\n",
+    "        col_dict_list.append(col_dict)\n",
+    "\n",
+    "        # plot some lines to show the colour and mean reference\n",
+    "        ax1.vlines(df_obs.iloc[xi][x_plot], df_obs.iloc[xi][y_plot], col_dict[y_ref_col], color=\"k\", ls=\":\")\n",
+    "        ax1.hlines(col_dict[y_ref_col], col_dict[x1_ref_col], col_dict[x2_ref_col], color=\"k\", ls=\"--\")\n",
+    "\n",
+    "# store running colour parameters as a dataframe\n",
+    "df_col = pd.DataFrame(col_dict_list)\n",
+    "df_col = pd.concat([df_obs, df_col], axis=1)\n",
+    "\n",
+    "ax1.set_xlabel(x_plot)\n",
+    "ax1.set_ylabel(y_plot)\n",
+    "ax1.legend()\n",
+    "ax1.invert_yaxis()\n",
+    "\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a7734122",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# recorded colour parameters\n",
+    "df_col"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0cd867d9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Plot how the colour changes as a function of time\n",
+    "# find filt_obs - filt_ref of newest filt_obs observation to the mean of the previous N_ref filt_ref observations\n",
+    "# colour code by time diff between obs and most recent obs_ref\n",
+    "\n",
+    "x_plot = \"midPointMjdTai\"\n",
+    "y_plot = colour\n",
+    "y_plot_err = colErr\n",
+    "c_plot = delta_t_col\n",
+    "df_plot = df_col\n",
+    "\n",
+    "fig = plt.figure()\n",
+    "gs = gridspec.GridSpec(1, 1)\n",
+    "ax1 = plt.subplot(gs[0, 0])\n",
+    "\n",
+    "\n",
+    "s1 = ax1.scatter(df_plot[x_plot], df_plot[y_plot], c=df_plot[c_plot], zorder=3)\n",
+    "cbar1 = plt.colorbar(s1)\n",
+    "ax1.errorbar(df_plot[x_plot], df_plot[y_plot], df_plot[yerr_plot], fmt=\".\", zorder=1)\n",
+    "\n",
+    "# TODO: leave the most recent obs out of this mean?\n",
+    "# loop through new obs and calculate current mean and std?\n",
+    "obs_ref_mean = np.mean(df_plot[y_plot])\n",
+    "obs_ref_median = np.median(df_plot[y_plot])\n",
+    "obs_ref_std = np.std(df_plot[y_plot])\n",
+    "print(obs_ref_median, obs_ref_std)\n",
+    "\n",
+    "ax1.axhline(obs_ref_mean, c=\"k\")\n",
+    "ax1.axhspan(obs_ref_mean - obs_ref_std, obs_ref_mean + obs_ref_std, zorder=0, color=\"k\", alpha=0.2)\n",
+    "\n",
+    "ax1.set_xlabel(x_plot)\n",
+    "ax1.set_ylabel(y_plot)\n",
+    "cbar1.set_label(c_plot)\n",
+    "\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "05abfee3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# These colours can then be run through the previously written outlier detection functions"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "605a043a",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "adler-dev",
+   "language": "python",
+   "name": "adler-dev"
+  },
+  "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.10.14"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}