tutorials/extragalactic_gcr_cluster_members.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plotting Positions of Galaxy Cluster Members in Extragalactic Catalogs\n",
    "\n",
    "In this example script we show how to select for and look at members of individual clusters.\n",
    "\n",
    "Owners: **Dan Korytov [@dkorytov](https://github.com/LSSTDESC/DC2-analysis/issues/new?body=@dkorytov)**, **Patricia Larsen**\n",
    "\n",
    "Last verified run: **Nov 30, 2018** by @yymao\n",
    "\n",
    "This notebook demonstrates how to access the extra galactic catalog through the Generic Catalog Reader (GCR, https://github.com/yymao/generic-catalog-reader) as well as how filter on galaxy features and cluster membership.\n",
    "\n",
    "__Objectives__:\n",
    "\n",
    "After working through and studying this Notebook you should be able to\n",
    "\n",
    "1. Access extragalactic catalogs (protoDC2, cosmoDC2) through the GCR.\n",
    "2. Select galaxy cluster centrals as a proxy for clusters.\n",
    "3. Select galaxies in individual clusters by using the host_id quantity.\n",
    "4. Plotting galaxy clustermembers positions on the sky as well as their comoving position in space.\n",
    "\n",
    "\n",
    "__Logistics__: This notebook is intended to be run through the JupyterHub NERSC interface available here: https://jupyter-dev.nersc.gov. To setup your NERSC environment, please follow the instructions available here: https://confluence.slac.stanford.edu/display/LSSTDESC/Using+Jupyter-dev+at+NERSC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import GCRCatalogs\n",
    "import numpy as np\n",
    "from astropy.table import Table\n",
    "from GCR import GCRQuery\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reading catalog\n",
    "We load in the catalog with the \"load_catalog\" command, and then the values with the \"get_quantities\" command using filters to select sub-samples of the catalog.  \n",
    "\n",
    "### Help for error messages:\n",
    "If this fails to find the appropriate quantities, check that the desc-python kernel is being used and if this is not available source the kernels by running the following command on a terminal at nersc: \"source /global/common/software/lsst/common/miniconda/setup_current_python.sh\"\n",
    "\n",
    "We are loading in a smaller version of the full cosmoDC2 catalog - this contains the same information as the full catalog but with a smaller sky area."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gc = GCRCatalogs.load_catalog('cosmoDC2_v1.1.4_small')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "galaxy_data = gc.get_quantities(['ra', 'dec', 'mag_r', 'halo_id'], filters=['mag_r < 19'])\n",
    "cluster_data = gc.get_quantities(['ra','dec', 'halo_mass', 'halo_id'], \n",
    "                                 filters=['is_central', 'halo_mass > 1e14', 'redshift < 0.2'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot individual clusters\n",
    "The cluster_data we loaded selects for central galaxies in clusters of mass higher than 10^14 solar masses/h below redshift 0.2. We will plot three such clusters in the proto-DC2 region. \n",
    "\n",
    "For each of these clusters we can match the galaxies by selecting on 'halo_id' and make a scatter plot of the projected positions around the central galaxy. Note that the size of the circles corresponds to the r-band magnitude of the galaxies. \n",
    "\n",
    "The cluster data is cast into a familiar astropy table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cluster_data = Table(cluster_data)\n",
    "for i, cluster in enumerate(cluster_data):\n",
    "    if (i >= 3):\n",
    "        break # plot only the first 3\n",
    "    members = GCRQuery('halo_id == {}'.format(cluster['halo_id'])).filter(galaxy_data)\n",
    "    plt.figure()\n",
    "    plt.scatter(\n",
    "        members['ra'], \n",
    "        members['dec'], \n",
    "        s=(19-members['mag_r'])*8, \n",
    "        label='Galaxy Members [{}]'.format(len(members['ra']))\n",
    "    )\n",
    "    plt.plot(cluster['ra'], cluster['dec'], 'xr', label='Cluster Center')\n",
    "    plt.legend(loc='best', framealpha=0.3)\n",
    "    plt.xlabel(r'ra [deg]')\n",
    "    plt.ylabel(r'dec [deg]')\n",
    "    plt.title('Halo ID:  {}\\nHalo Mass:  {:.2e} h^-1 Msun'.format(cluster['halo_id'], cluster['halo_mass']))\n",
    "        \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extensions:\n",
    "We can load further information on the cluster members. For instance the second cluster looks a little odd in projected space, so we re-make these plots in comoving cartesian coordinates x and y in the example below. We also map the colours to the x-direction velocities. \n",
    "\n",
    "To do this you need to load the required quantities from the catalog before using them. A simple way to double check the quantity names is the command \"gc.list_all_quantities()\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "galaxy_data = gc.get_quantities(['ra', 'dec', 'mag_r', 'halo_id', 'position_x', 'position_y', 'velocity_x', 'velocity_y'], filters=['mag_r < 19'])\n",
    "cluster_data = gc.get_quantities(['ra','dec', 'halo_mass', 'halo_id', 'position_x', 'position_y'], \n",
    "                                 filters=['is_central', 'halo_mass > 1e14', 'redshift < 0.2'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cluster_data = Table(cluster_data)\n",
    "for i, cluster in enumerate(cluster_data):\n",
    "    if (i >= 3):\n",
    "        break # plot only the first 3\n",
    "    members = GCRQuery('halo_id == {}'.format(cluster['halo_id'])).filter(galaxy_data)\n",
    "    plt.figure()\n",
    "    plt.scatter(\n",
    "        members['position_x'],\n",
    "        members['position_y'],\n",
    "        s=(19-members['mag_r'])*8, \n",
    "        label='Galaxy Members [{}]'.format(len(members['ra'])),\n",
    "        c=members['velocity_x'],\n",
    "        cmap='viridis')\n",
    "    plt.plot(cluster['position_x'], cluster['position_y'],'xr',label='Cluster Center', alpha=0.8)\n",
    "    plt.legend(loc='best',framealpha=0.3)\n",
    "    plt.xlabel(r'x [Mpc/h]')\n",
    "    plt.ylabel(r'y [Mpc/h]')\n",
    "    plt.title('Halo ID: {}\\nHalo Mass: {:.2e}'.format(cluster['halo_id'], cluster['halo_mass']))\n",
    "    plt.colorbar(label='km/s')\n",
    "plt.show()"
   ]
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}