{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "0789d321-e15e-44b8-b96e-34b4e993cf16"
    }
   },
   "source": [
    "# Looking at the \"Object Catalogs\": merged tract-patch catalogs in DC2 Run 1.1p\n",
    "<br>Owner: **Michael Wood-Vasey** ([@wmwv](https://github.com/LSSTDESC/DC2-analysis/issues/new?body=@wmwv))\n",
    "<br>Last Verified to Run: **2018-07-24**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Purpose:\n",
    "Demonstrate reading a static-sky object catalog from photometry on the coadd images based on the static coadd catalog.  Shows the example of plotting a color-color stellar locus against the data to demonstrate some basic use of the magnitudes and extendedness columns.\n",
    "\n",
    "### Learning Objectives:\n",
    "After working through and studying this Notebook you should be able to\n",
    "1. Read in the DC2 Run 1.1p object catalog directly from the flat files.\n",
    "2. Construct a color-color diagram for stars\n",
    "3. Overplot a stellar locus model\n",
    "\n",
    "### Logistics\n",
    "This is intended to be runnable at NERSC through the https://jupyter-dev.nersc.gov interface from a local git clone of https://github.com/LSSTDESC/DC2-analysis in your NERSC directory.  But you can also run it wherever, with appropriate adjustment of the 'repo' location to point to a place where you have a Butler repo will all of the images. Instructions for setting up the proper python kernel can be found here: https://confluence.slac.stanford.edu/x/1_ubDQ\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "69694b18-dcc8-463d-8e86-e4fcd572f680"
    }
   },
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "%matplotlib inline \n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# cmap = 'Oranges'\n",
    "cmap = 'viridis_r'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "00c38da4-1a09-46e2-a3dc-b63886072989"
    }
   },
   "outputs": [],
   "source": [
    "merged_tract_data_dir = '/global/projecta/projectdirs/lsst/global/in2p3/Run1.1/object_catalog/'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the LSST Data Management processing, the sky is divided into tracts.  Each tract is divided into patches (8x8 in this processing).  The object catalog files are saved by tract.  These files are HDF5 files which are divided into group keys by patch.\n",
    "\n",
    "For more information see: [How to the Run 1.1p skymap](Plotting_the_Run1.1p_skymap.ipynb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "1d68d81a-790f-460d-900c-3f447c8309a9"
    }
   },
   "outputs": [],
   "source": [
    "tract, patch = 4849, 16\n",
    "basename = 'merged_tract_%d.hdf5' % tract\n",
    "key = 'coadd_%d_%2d' % (tract, patch)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "ee575d25-f444-4410-9b18-9d263b4d43fd"
    }
   },
   "outputs": [],
   "source": [
    "merged_tract_file = os.path.join(merged_tract_data_dir, basename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "a37e0a04-e349-4229-ba8a-bb96f857d68e"
    }
   },
   "outputs": [],
   "source": [
    "df = pd.read_hdf(merged_tract_file, key=key)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "f27a27c2-5123-420c-8cce-a0242b0b9966"
    }
   },
   "outputs": [],
   "source": [
    "print(len(df), \"objects\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "e1b3faf8-518f-4162-853e-7da255b2d7bd"
    }
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8,8))\n",
    "plt.scatter(df['g_mag'] - df['r_mag'], df['r_mag'] - df['i_mag'], marker='.')\n",
    "plt.xlim(-1, +2)\n",
    "plt.ylim(-1, +2)\n",
    "plt.xlabel('g-r')\n",
    "plt.ylabel('r-i')\n",
    "plt.title('Color-color diagram of all objects in one tract of object catalogs')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "4527a23f-081f-4b46-88cf-407d5935fd2d"
    }
   },
   "source": [
    "This plot includes all objects, including lots of low SNR points and so is a bit of a mess, but maybe we're starting to see a bit of some structure resolving.  Let's do the traditional simple 2D histogram version:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's overplot a stellar locus from [Davenport et al., 2014, MNRAS, 440, 3430](http://adsabs.harvard.edu/abs/2014MNRAS.440.3430D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def simple_stellar_locus_gmr_rmi(color='red', linestyle='--', linewidth=2.5):\n",
    "    \"\"\"Provide a simple two-slope definition of a stellar locus in (g-r, r-i).\n",
    "    \n",
    "    This simply defines one vertical line for M stars\n",
    "    and the general slope for the bluer stars.\"\"\"\n",
    "    m_stars_gmr = [+1.4, +1.4]  # mag\n",
    "    m_stars_rmi = [+0.5, +1.5]\n",
    "\n",
    "    other_stars_gmr = [-0.5, +1.4]\n",
    "    other_stars_rmi = [-0.4, +0.5]\n",
    "\n",
    "    model_gmr = other_stars_gmr + m_stars_gmr\n",
    "    model_rmi = other_stars_rmi + m_stars_rmi\n",
    "    \n",
    "    return model_gmr, model_rmi\n",
    "\n",
    "    \n",
    "def get_stellar_locus_davenport(color1='gmr', color2='rmi',\n",
    "                                datafile='assets/Davenport_2014_MNRAS_440_3430_table1.txt'):\n",
    "    data = pd.read_table(datafile, sep='\\s+', header=1)\n",
    "    return data[color1], data[color2]\n",
    "    \n",
    "    \n",
    "def plot_stellar_locus(color1='gmr', color2='rmi',\n",
    "                       color='red', linestyle='--', linewidth=2.5):\n",
    "    model_gmr, model_rmi = get_stellar_locus_davenport(color1, color2)\n",
    "    plot_kwargs = {'linestyle': linestyle, 'linewidth': linewidth, 'color': color,\n",
    "                   'scalex': False, 'scaley': False}\n",
    "    plt.plot(model_gmr, model_rmi, **plot_kwargs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_color_color(z, color1, color2, range1=(-1, +2), range2=(-1, +2), bins=31):\n",
    "    \"\"\"Plot a color-color diagram.  Overlay stellar locus\"\"\"\n",
    "    band1, band2 = color1[0], color1[-1]\n",
    "    band3, band4 = color2[0], color2[-1]\n",
    "    H, xedges, yedges = np.histogram2d(\n",
    "        z['%s_mag' % band1] - z['%s_mag' % band2],\n",
    "        z['%s_mag' % band3] - z['%s_mag' % band4],\n",
    "        range=(range1, range2), bins=bins)\n",
    "        \n",
    "    zi = H.T\n",
    "    xi = (xedges[1:] + xedges[:-1])/2\n",
    "    yi = (yedges[1:] + yedges[:-1])/2\n",
    "\n",
    "    plt.figure(figsize=(8, 8))\n",
    "    plt.pcolormesh(xi, yi, zi, cmap=cmap)\n",
    "    plt.contour(xi, yi, zi)\n",
    "    plt.xlabel('%s-%s' % (band1, band2))\n",
    "    plt.ylabel('%s-%s' % (band3, band4))\n",
    "\n",
    "    plot_stellar_locus(color1, color2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "3a89fe4d-718f-4748-9ca7-28e0b805e04d"
    }
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8,6.5))\n",
    "plt.hist2d(df['g_mag']-df['r_mag'], df['r_mag']-df['i_mag'],\n",
    "           cmap=cmap, range=((-1, +2), (-1, +2)), bins=51)\n",
    "plt.colorbar()\n",
    "plt.xlabel('g-r')\n",
    "plt.ylabel('r-i')\n",
    "\n",
    "plot_stellar_locus('gmr', 'rmi')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(Don't worry about the runtime error.  There are some NaNs in the catalog magnitudes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Restrict to \"good\" sources\n",
    "Clean this up a bit to define\n",
    " 1. objects with SNR > 25 in g and r.\n",
    " 2. stars from this SNR > 25 sample.\n",
    " 3. bright-star set consistent with the usage in the LSST SRD of stars with SNR~100 (although that threshold may arguably be meant to be on the SNR in individual images).\n",
    " \n",
    "The DM Science Pipeline calculates an 'extendness', which is the degree to which an object is classified as 'extended', i.e., not a point source (star).  To select stars, we select for objects with an extendeness classification less than 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "snr_threshold = 25\n",
    "mag_err_threshold = 1/snr_threshold\n",
    "good_snr = df[(df['g_mag_err'] < mag_err_threshold) & (df['r_mag_err'] < mag_err_threshold)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bright_snr_threshold = 100\n",
    "mag_err_threshold = 1/bright_snr_threshold\n",
    "bright_stars = df[(df['g_mag_err'] < mag_err_threshold) & (df['r_mag_err'] < mag_err_threshold)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "3a0b17cb-981e-49a3-a872-1f5ef2ab59cd"
    }
   },
   "outputs": [],
   "source": [
    "safe_max_extended = 1.0\n",
    "stars = good_snr[good_snr['base_ClassificationExtendedness_value'] < safe_max_extended]\n",
    "galaxies = good_snr[good_snr['base_ClassificationExtendedness_value'] >= safe_max_extended]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"%d stars (SNR > %.0f)\" % (len(stars), snr_threshold))\n",
    "print(\"%d bright stars (SNR > %.0f)\" % (len(bright_stars), bright_snr_threshold))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "663149ff-4774-43d3-be5a-0a3bf8ed8952"
    }
   },
   "outputs": [],
   "source": [
    "hist_kwargs = {'bins': np.linspace(-2, 3, 51),\n",
    "               'range': (-2, +3),\n",
    "               'linewidth': 4,\n",
    "               'histtype': 'step'}\n",
    "          \n",
    "plt.hist(df['g_mag'] - df['r_mag'], label='all', color='grey', **hist_kwargs)\n",
    "plt.hist(good_snr['g_mag'] - good_snr['r_mag'], label='SNR > %.0f' % snr_threshold, color='black',\n",
    "         **hist_kwargs)\n",
    "plt.hist(galaxies['g_mag'] - galaxies['r_mag'], label='galaxies', **hist_kwargs)\n",
    "plt.hist(stars['g_mag'] - stars['r_mag'], label='stars', **hist_kwargs)\n",
    "plt.xlabel('g-r')\n",
    "plt.ylabel('# / 0.5 mag bin')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "dacd2cd2-204b-43b9-9f68-447b52b581fd"
    }
   },
   "outputs": [],
   "source": [
    "plot_color_color(galaxies, 'gmr', 'rmi')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_color_color(stars, 'gmr', 'rmi')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. You can see that we still have some galaxies, but otherwise we're starting to see the stars following the red line.\n",
    "2. There's a slight offset with the red line, but we haven't really worried exactly about dust, metallicity, or filter transmission functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_color_color(stars, 'rmi', 'imz')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's briefly look at the (inverse) magnitude error vs. magnitude to see if the behavior seems reasonable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(df['r_mag'], 1/df['r_mag_err'], label='1/r_mag_err')\n",
    "plt.xlabel('r')\n",
    "plt.ylabel('1/r_err')\n",
    "plt.axhline(bright_snr_threshold, linestyle='--', color='red', label='Bright star')\n",
    "plt.axhline(snr_threshold, linestyle='--', color='orange', label='Good SNR')\n",
    "plt.legend()\n",
    "plt.xlim(15, 28)\n",
    "# plt.ylim(0, 1.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I don't fully understand this above plot."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Color-Magnitude plot\n",
    "\n",
    "Stars and galaxies also occupy different, but overlapping spaces in color-magnitude plots.  Here are SDSS color-magnitude and color-color plots for g, r, i.  Above we looked at the lower, color-color plots.\n",
    "\n",
    "http://www.astroml.org/examples/datasets/plot_sdss_imaging.html\n",
    "![](http://www.astroml.org/_images/plot_sdss_imaging_1.png)\n",
    "\n",
    "Let's see what we get from our simulations->images->DM science pipeline processing for the color-magnitude plots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z = galaxies\n",
    "plt.scatter(z['g_mag'] - z['r_mag'], z['g_mag'], marker='.', label='galaxies')\n",
    "plt.xlabel('g-r')\n",
    "plt.ylabel('g')\n",
    "plt.xlim(-1, +2)\n",
    "plt.ylim(25, 15)\n",
    "plt.title('Galaxies')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare to the upper-left plot from the AtroML visualization of the SDSS data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z = stars\n",
    "H, xedges, yedges = np.histogram2d(z['g_mag']-z['r_mag'], z['g_mag'],\n",
    "                                    range=((-1, +2), (15, 25)), bins=31)\n",
    "zi = H.T\n",
    "xi = (xedges[1:] + xedges[:-1])/2\n",
    "yi = (yedges[1:] + yedges[:-1])/2\n",
    "\n",
    "plt.pcolormesh(xi, yi, zi, cmap=cmap)\n",
    "plt.contour(xi, yi, zi)\n",
    "plt.xlabel('g-r')\n",
    "plt.ylabel('g')\n",
    "plt.ylim(25, 15)\n",
    "plt.title('Stars')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare to the upper-right plot from the AtroML visualization of the SDSS data.  It looks reasonable, but we still have significant galaxy contamination, as we similarly saw in the color-color plots above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "4791c937-9435-405d-8490-a4f9f6d80aa6"
    }
   },
   "source": [
    "## Appendix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generic Catalog Reader (GCR)\n",
    "\n",
    "After a little bit you'll find yourself wanting to ask interesting questions, like how well does the DC2 analysis recover the input sources?  Or easily loading in all of the patches of a tract.  A convenient framework for generalizing analyses of DESC-DC related catalogs is the Generic Catalog Reader:\n",
    "\n",
    "* GCR is a general utility\n",
    "https://github.com/yymao/generic-catalog-reader\n",
    "* The catalogs supporting DC2 are in:\n",
    "https://github.com/LSSTDESC/gcr-catalogs\n",
    "\n",
    "See the \"DC2 object Run1.1p GCR tutorial\" series!\n",
    "For an example of using the GCR to do work similar to the agove:\n",
    "[DC2 object Run1.1p GCR tutorial -- Part I: GCR Access.ipynb](object_gcr_1_intro.ipynb)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "desc-python",
   "language": "python",
   "name": "desc-python"
  },
  "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.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}