diff --git a/README.md b/README.md index 1006e25..32cc8df 100644 --- a/README.md +++ b/README.md @@ -74,11 +74,13 @@ Notes: ## Dev Guide - Adding notebooks to Read The Docs -- Copy notebook into `docs/notebooks` (N.B. the notebook must have at least one section header and be using the "Python 3 (ipykernel)" kernel, not some conda env kernel that may only be installed locally) +- Copy notebook into `docs/notebooks` (N.B. the notebook must have at least one section header* and be using the "Python 3 (ipykernel)" kernel, not some conda env kernel that may only be installed locally) - Update the toctree in the file `docs/notebooks.rst` - Ensure necessary requirements are declared in `pyproject.toml` and `docs/requirements.txt`. Also, make sure that the notebook being added to the docs is using the python3 (ipykernel) kernel, not some conda env kernel that may only be installed locally - To update the docs locally, from the `docs` dir run: `python -m sphinx -T -E -b html -d _build/doctrees -D language=en . ../_readthedocs/html` +* Multiple section headers from a notebook will also show up in the table of contents. + ## Dev Guide - Updating pyproject.toml If you are adding code that requires a new dependency, this needs to be included in pyproject.toml under the `[project]' section: diff --git a/docs/index.rst b/docs/index.rst index ec57dc4..f76c230 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -5,6 +5,11 @@ Welcome to adler's documentation! ======================================================================================== +Introduction +------------ + +Adler is an open source package which provides tools for detecting and characterising activity in observations of solar system objects. + Dev Guide - Getting Started --------------------------- diff --git a/docs/notebooks.rst b/docs/notebooks.rst index 61e3308..44fd5bf 100644 --- a/docs/notebooks.rst +++ b/docs/notebooks.rst @@ -5,3 +5,6 @@ Notebooks Introducing Jupyter Notebooks Adler phasecurve models + Adler plotting example + Adler outlier detection + Adler colour measurement \ No newline at end of file diff --git a/docs/notebooks/colour_functions_example.ipynb b/docs/notebooks/colour_functions_example.ipynb new file mode 100644 index 0000000..ec57c97 --- /dev/null +++ b/docs/notebooks/colour_functions_example.ipynb @@ -0,0 +1,391 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "6ebe3e22", + "metadata": {}, + "source": [ + "# Adler colour measurement\n", + "\n", + "Adler includes a function `col_obs_ref` which calculates the difference between the most recent brightness data point in the observation (obs) filter and a number of previous measurements in the reference (ref) filter." + ] + }, + { + "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" + ] + }, + { + "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 = \"../../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_data = AdlerData(ssoid, planetoid.filter_list)\n", + "\n", + "for filt in [\"g\", \"r\"]:\n", + " # get observations and define a phase curve model to fit\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_data.populate_phase_parameters(filt, **pc.ReturnModelDict())\n", + "\n", + " # add this phase curve to the figure using the Adler plotting function\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_data.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_data.get_phase_parameters_in_filter(\"g\", \"HG12_Pen16\").__dict__" + ] + }, + { + "cell_type": "markdown", + "id": "e07e86a8", + "metadata": {}, + "source": [ + "**Determine the apparitions (periods of observability) of the object.**\n", + "\n", + "Get the boundary times for each apparation of the object in the survey using the Adler helper function `apparition_gap_finder`.\n", + "In this example we will just look at changes in colour for a single apparition" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "73871e56", + "metadata": {}, + "outputs": [], + "source": [ + "# combine all measurements in r and g into one dataframe as apparitions are filter independent\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", + "# get the boundary times\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": "markdown", + "id": "08f25eb8", + "metadata": {}, + "source": [ + "Here we simulate observations coming night-by-night and the calculation of a g-r colour for the object\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bd0e0185", + "metadata": {}, + "outputs": [], + "source": [ + "# define colour function parameters\n", + "\n", + "# set number of reference observations to use for colour estimate, there are multiple options\n", + "# N_ref = 5\n", + "N_ref = 3\n", + "# N_ref = 1\n", + "# N_ref = None # selecting None uses all previous reference filter measurements\n", + "\n", + "# observation and filter field names\n", + "x_plot = \"midPointMjdTai\" # time column\n", + "y_plot = \"reduced_mag\" # magnitude column\n", + "yerr_plot = \"magErr\" # magnitude uncertainty column\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", + "\n", + " _df_obs_all = df_obs_all[\n", + " (df_obs_all[\"midPointMjdTai\"] >= time_min) & (df_obs_all[\"midPointMjdTai\"] < time_max)\n", + " ]\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_data.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_data.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 new filt_obs observation\n", + " x1 = time_min\n", + " for xi in range(len(df_obs)):\n", + " x2 = df_obs.iloc[xi][x_plot]\n", + "\n", + " # run the colour finding function here\n", + " col_dict = col_obs_ref(\n", + " planetoid,\n", + " adler_data,\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": [ + "# display the recorded colour parameters\n", + "df_col" + ] + }, + { + "cell_type": "markdown", + "id": "fae0b512", + "metadata": {}, + "source": [ + "Now we can plot how the colour changes as a function of time" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0cd867d9", + "metadata": {}, + "outputs": [], + "source": [ + "# 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", + "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", + "obs_ref_mean = np.mean(df_plot[y_plot])\n", + "obs_ref_std = np.std(df_plot[y_plot])\n", + "print(\"{}-{} mean = {}, std = {}\".format(filt_obs, filt_ref, obs_ref_mean, 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": "markdown", + "id": "2aad6ee0", + "metadata": {}, + "source": [ + "These colours can then be run through the previously written outlier detection functions.\n", + "We have recorded metadata which can help exclude erroneous colour measurements, such as the time difference between the obs and ref measurements." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "605a043a", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "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.9.19" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/docs/notebooks/outlier_detection_example.ipynb b/docs/notebooks/outlier_detection_example.ipynb new file mode 100644 index 0000000..984375d --- /dev/null +++ b/docs/notebooks/outlier_detection_example.ipynb @@ -0,0 +1,304 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "5b6823e7", + "metadata": {}, + "source": [ + "# Adler outlier detection\n", + "This notebook demonstrates some of the functions provide by Adler to assist with simple outlier detection" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f88d62a9", + "metadata": {}, + "outputs": [], + "source": [ + "from adler.dataclasses.AdlerPlanetoid import AdlerPlanetoid\n", + "from adler.science.PhaseCurve import PhaseCurve\n", + "from adler.dataclasses.AdlerData import AdlerData\n", + "import adler.utilities.science_utilities as utils\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" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "86e657fe", + "metadata": {}, + "outputs": [], + "source": [ + "# ssObjectId of object to analyse\n", + "ssoid = \"8268570668335894776\"\n", + "\n", + "# load object from local database\n", + "fname = \"../../notebooks/gen_test_data/adler_demo_testing_database.db\"\n", + "planetoid = AdlerPlanetoid.construct_from_SQL(ssoid, sql_filename=fname)\n", + "\n", + "# retrieve observations in the r filter\n", + "obs_r = planetoid.observations_in_filter(\"r\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f4131a10", + "metadata": {}, + "outputs": [], + "source": [ + "# define the phase curve model using the SSObject data\n", + "\n", + "sso_r = planetoid.SSObject_in_filter(\"r\")\n", + "r_H = sso_r.H\n", + "r_G12 = sso_r.G12\n", + "\n", + "pc = PhaseCurve(H=r_H * u.mag, phase_parameter_1=r_G12, model_name=\"HG12_Pen16\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0f625b5a", + "metadata": {}, + "outputs": [], + "source": [ + "# calculate data minus model residuals\n", + "res = obs_r.reduced_mag - pc.ReducedMag(obs_r.phaseAngle * u.degree).value" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "81e04bad", + "metadata": {}, + "outputs": [], + "source": [ + "res" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f573cce4", + "metadata": {}, + "outputs": [], + "source": [ + "# plot the observations with the SSObject phase curve\n", + "x_plot = \"phaseAngle\"\n", + "y_plot = \"reduced_mag\"\n", + "\n", + "x = getattr(obs_r, x_plot)\n", + "y = getattr(obs_r, y_plot)\n", + "yerr = obs_r.magErr\n", + "\n", + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "ax1.errorbar(x, y, yerr, fmt=\"o\")\n", + "\n", + "# plot the phase curve model\n", + "alpha = np.linspace(0, np.amax(obs_r.phaseAngle)) * u.deg\n", + "red_mag = pc.ReducedMag(alpha)\n", + "\n", + "# legend label for the phase curve model\n", + "pc_label = []\n", + "for x in pc.model_function.param_names:\n", + " pc_label.append(\"{}={:.2f}\".format(x, getattr(pc.model_function, x).value))\n", + "pc_label = \", \".join(pc_label)\n", + "\n", + "ax1.plot(alpha.value, red_mag.value, c=\"k\", label=pc_label)\n", + "\n", + "ax1.invert_yaxis()\n", + "ax1.set_xlabel(x_plot)\n", + "ax1.set_ylabel(y_plot)\n", + "ax1.legend()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "aee95371", + "metadata": {}, + "outputs": [], + "source": [ + "# plot the data - model residuals\n", + "x_plot = \"phaseAngle\"\n", + "\n", + "x = getattr(obs_r, x_plot)\n", + "yerr = obs_r.magErr\n", + "\n", + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "ax1.errorbar(x, res, yerr, fmt=\"o\")\n", + "\n", + "ax1.axhline(0, c=\"k\")\n", + "\n", + "# indicate the standard deviations of the residuals\n", + "res_std = np.std(res)\n", + "for i in range(1, 4):\n", + " ax1.axhline(res_std * i, ls=\":\", c=\"C{}\".format(i), label=\"{} sigma\".format(i))\n", + " ax1.axhline(-res_std * i, ls=\":\", c=\"C{}\".format(i))\n", + "\n", + "ax1.invert_yaxis()\n", + "ax1.set_xlabel(x_plot)\n", + "ax1.set_ylabel(\"data - model\")\n", + "ax1.legend()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "4318568a", + "metadata": {}, + "source": [ + "**Return a list of flags for outlying objects.**\n", + "\n", + "The Adler `utils.sigma_clip` function is a wrapper for `astropy.stats.sigma_clip`. We do this in order to return just the clip mask, and also to make it easier to call a \"zero\" central function." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "16c712ff", + "metadata": {}, + "outputs": [], + "source": [ + "# astropy sigma_clip normally uses the median as the central function\n", + "utils.sigma_clip(res, {\"cenfunc\": \"median\"})" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8983c9dd", + "metadata": {}, + "outputs": [], + "source": [ + "# assuming that the model is the ground truth, we use zero as the centroid for the residuals\n", + "utils.sigma_clip(res, {\"cenfunc\": utils.zero_func})" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "15bdb74e", + "metadata": {}, + "outputs": [], + "source": [ + "# use the standard deviation of the residuals to identify outliers\n", + "utils.outlier_std(res, res, std_cut=3.0)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "50276ce5", + "metadata": {}, + "outputs": [], + "source": [ + "# use a simple threshold value for residuals to find outliers\n", + "utils.outlier_diff(res, diff_cut=1.0)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "64e2f5bc", + "metadata": {}, + "outputs": [], + "source": [ + "# consider the residual compared to the uncertainty of the measurement\n", + "std_err = 5\n", + "utils.outlier_sigma_diff(res, yerr, std_err)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3f8cc586", + "metadata": {}, + "outputs": [], + "source": [ + "# plot the data - model residuals\n", + "x_plot = \"phaseAngle\"\n", + "\n", + "x = getattr(obs_r, x_plot)\n", + "yerr = obs_r.magErr\n", + "\n", + "fig = plt.figure()\n", + "gs = gridspec.GridSpec(1, 1)\n", + "ax1 = plt.subplot(gs[0, 0])\n", + "\n", + "ax1.errorbar(x, res, yerr, fmt=\"o\")\n", + "\n", + "ax1.axhline(0, c=\"k\")\n", + "\n", + "# indicate the standard deviations of the residuals\n", + "res_std = np.std(res)\n", + "for i in range(1, 4):\n", + " ax1.axhline(res_std * i, ls=\":\", c=\"C{}\".format(i), label=\"{} sigma\".format(i))\n", + " ax1.axhline(-res_std * i, ls=\":\", c=\"C{}\".format(i))\n", + "\n", + "mask = utils.outlier_sigma_diff(res, yerr, std_err)\n", + "ax1.scatter(x[mask], res[mask], edgecolor=\"r\", facecolor=\"none\", marker=\"o\", zorder=3)\n", + "\n", + "ax1.invert_yaxis()\n", + "ax1.set_xlabel(x_plot)\n", + "ax1.set_ylabel(\"data - model\")\n", + "ax1.legend()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "e66bfbee", + "metadata": {}, + "source": [ + "NB that for phase curve models, residuals can be much larger than the photometric uncertainty!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bcccc758", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "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.9.19" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/docs/notebooks/plotting_utilities_example.ipynb b/docs/notebooks/plotting_utilities_example.ipynb new file mode 100644 index 0000000..c18cccd --- /dev/null +++ b/docs/notebooks/plotting_utilities_example.ipynb @@ -0,0 +1,259 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "d1c2637e", + "metadata": {}, + "source": [ + "# Adler plotting utilities\n", + "For ease of use Adler provides helper functions which can be used to make common solar system object plots, such as light curves and phase curves." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9f5f4b7d", + "metadata": {}, + "outputs": [], + "source": [ + "from adler.dataclasses.AdlerPlanetoid import AdlerPlanetoid\n", + "from adler.science.PhaseCurve import PhaseCurve\n", + "from adler.utilities.plotting_utilities import plot_errorbar\n", + "\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import matplotlib.gridspec as gridspec\n", + "import astropy.units as u" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "28113970", + "metadata": {}, + "outputs": [], + "source": [ + "# retrieve the object data\n", + "ssoid = \"8268570668335894776\"\n", + "fname = \"../../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": "f27e1dec", + "metadata": {}, + "outputs": [], + "source": [ + "# use the adler plotting function to create a phase curve from the planetoid object\n", + "fig = plot_errorbar(planetoid, filt_list=[\"r\"])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fd41e5bf", + "metadata": {}, + "outputs": [], + "source": [ + "# with matplotlib we can access axes properties and update them after the fact\n", + "\n", + "# ax1.__dict__" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "72ec3def", + "metadata": {}, + "outputs": [], + "source": [ + "# access the axes object to update attributes\n", + "ax1 = fig.axes[0]\n", + "ax1.set_xlabel(\"phaseAngle (degrees)\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3b751ba2", + "metadata": {}, + "outputs": [], + "source": [ + "# replot the figure\n", + "fig" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2ee0ae9e", + "metadata": {}, + "outputs": [], + "source": [ + "# fit a phase curve model to the data\n", + "\n", + "filt = \"r\"\n", + "sso = planetoid.SSObject_in_filter(filt)\n", + "obs = planetoid.observations_in_filter(filt)\n", + "\n", + "H = sso.H\n", + "G12 = sso.G12\n", + "\n", + "pc = PhaseCurve(H=H * u.mag, phase_parameter_1=G12, model_name=\"HG12_Pen16\")\n", + "alpha = np.linspace(0, np.amax(obs.phaseAngle)) * u.deg\n", + "red_mag = pc.ReducedMag(alpha)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c8c4608d", + "metadata": {}, + "outputs": [], + "source": [ + "# add this phase curve to the figure\n", + "ax1.plot(alpha.value, pc.ReducedMag(alpha).value, label=\"{} {}\".format(filt, pc.model_name))\n", + "ax1.legend() # udpate the figure legend" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ab0df0b1", + "metadata": {}, + "outputs": [], + "source": [ + "fig" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6cae1b45", + "metadata": {}, + "outputs": [], + "source": [ + "# we can also pass the fig object to the plotting function again to add more data\n", + "fig2 = plot_errorbar(planetoid, fig=fig, filt_list=[\"g\", \"i\"], label_list=[\"g\", \"i\"])\n", + "\n", + "# update the legend\n", + "ax1 = fig2.axes[0]\n", + "ax1.legend()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3d92c46d", + "metadata": {}, + "outputs": [], + "source": [ + "fig2" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f5a183bb", + "metadata": {}, + "outputs": [], + "source": [ + "# inspect the different items that have been plotted" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f356cbf7", + "metadata": {}, + "outputs": [], + "source": [ + "ax1._children" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4bf0f661", + "metadata": {}, + "outputs": [], + "source": [ + "ax1.containers" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c847b938", + "metadata": {}, + "outputs": [], + "source": [ + "# add the legend label to the r filter data\n", + "ax1.containers[0]._label = \"r\"\n", + "ax1.legend()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "07982dd4", + "metadata": {}, + "outputs": [], + "source": [ + "fig2" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a56d5711", + "metadata": {}, + "outputs": [], + "source": [ + "# we use `plot_errorbar` to save the figure, without adding anything extra to the figure\n", + "fig3 = plot_errorbar(planetoid, fig=fig2, filename=\"phase_curve_{}.png\".format(ssoid))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1c89f107", + "metadata": {}, + "outputs": [], + "source": [ + "fig3" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b76a0b12", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "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.9.19" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}