FC1-sat-pos-vel.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from itertools import combinations\n",
    "from collections import Counter, defaultdict\n",
    "from astropy.table import Table\n",
    "from astropy.coordinates import SkyCoord, Distance\n",
    "from SAGA.utils import add_skycoord\n",
    "from easyquery import Query, QueryMaker\n",
    "from astropy.wcs import WCS\n",
    "from astropy.io import fits\n",
    "import SAGA\n",
    "from SAGA.database import FitsTable\n",
    "from SAGA import ObjectCuts as C\n",
    "from SAGA.utils import add_skycoord, fill_values_by_query\n",
    "from SAGA.utils.distance import z2v\n",
    "from matplotlib.patches import Ellipse\n",
    "\n",
    "print(SAGA.__version__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "saga = SAGA.QuickStart()\n",
    "saga.set_default_base_version(\"paper2\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "base = SAGA.database.FitsTable(\"/home/yymao/Documents/Research/SAGA/PaperII/data-archive/saga_base_all.fits\").read()\n",
    "base = saga.host_catalog.construct_host_query(\"paper2_complete\").filter(base)\n",
    "assert len(np.unique(base[\"HOSTID\"])) == 36\n",
    "\n",
    "sats = Query(C.is_sat).filter(base)\n",
    "del base"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hostlist = saga.host_catalog.load(include_stats=True, query=\"paper2_complete\")[\"HOSTID\", \"DIST\", \"RA\", \"DEC\", \"COMMON_NAME\", \"coord\"]\n",
    "hostlist.sort(\"RA\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use ESO DSS tool to get the host background images\n",
    "# print(*[\" \".join([host[\"HOSTID\"], host[\"coord\"].to_string(style=\"hmsdms\", sep=\" \"), \"60 60\"]) for host in hostlist], sep=\"\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "from astropy.wcs import FITSFixedWarning\n",
    "from astropy.io.fits.verify import VerifyWarning\n",
    "warnings.filterwarnings(\"ignore\", category=FITSFixedWarning)\n",
    "warnings.filterwarnings(\"ignore\", category=VerifyWarning)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ncols = 6\n",
    "nrows = 6\n",
    "assert ncols*nrows >= len(hostlist)\n",
    "bg_keep = {}\n",
    "\n",
    "fig = plt.figure(figsize=(14, 14*nrows/ncols*1.1))\n",
    "gs = plt.GridSpec(nrows+1, ncols, fig, top=1, left=0, right=1, bottom=0, wspace=0.25, hspace=0.25, height_ratios=tuple([1]*nrows + [0.05]))\n",
    "\n",
    "for i, host in enumerate(hostlist):\n",
    "    \n",
    "    host_rvir = np.rad2deg(np.arcsin(0.3 / host[\"DIST\"]))\n",
    "    host_ra_corr = np.abs(np.cos(np.deg2rad(host[\"DEC\"])))\n",
    "    host_center = (host[\"RA\"], host[\"DEC\"])\n",
    "    host_name = host[\"COMMON_NAME\"]\n",
    "    fp = fits.open(\"/home/yymao/Documents/Research/SAGA/PaperII/host-dss-images/{}.fits\".format(host[\"HOSTID\"]), output_verify=\"silentfix\")\n",
    "    host_background = fp[0]\n",
    "    sats_this = QueryMaker.equals(\"HOSTID\", host[\"HOSTID\"]).filter(sats)\n",
    "    ec = plt.cm.cool((z2v(sats_this[\"SPEC_Z\"])-sats_this[\"HOST_VHOST\"])/275/2 + 0.5)\n",
    "\n",
    "    if host[\"COMMON_NAME\"] in [\"NGC3976\", \"NGC4158\", \"NGC5602\"]:\n",
    "        bg_keep[host[\"COMMON_NAME\"]] = host_background.data.copy()\n",
    "        k = host_background.data.ravel()\n",
    "        k[k > 17500] = 0\n",
    "        host_background.data[:] = k.reshape(host_background.data.shape)\n",
    "    \n",
    "    ax = plt.subplot(gs[i], projection=WCS(host_background.header), label=host[\"HOSTID\"])\n",
    "    ax.imshow(host_background.data, cmap='Greys')\n",
    "    ax.scatter(sats_this[\"RA\"], sats_this[\"DEC\"], transform=ax.get_transform('fk5'), s=35, c='none', lw=2, edgecolors=ec)\n",
    "    ax.add_patch(Ellipse((host[\"RA\"], host[\"DEC\"]), 2.0*host_rvir/host_ra_corr, 2.0*host_rvir, \n",
    "                 edgecolor='#999999', facecolor='none', transform=ax.get_transform('fk5'))) \n",
    "\n",
    "    ax.set_xlim(-0.5, host_background.data.shape[1] - 0.5)\n",
    "    ax.set_ylim(-0.5, host_background.data.shape[0] - 0.5)\n",
    "    ax.text(0.03, 1.03, host[\"COMMON_NAME\"], transform=ax.transAxes, va='bottom', fontweight=\"bold\", fontsize=18)\n",
    "    for c in ax.coords:\n",
    "        c.set_ticklabel(exclude_overlapping=True, size=9, pad=0)\n",
    "        c.display_minor_ticks(True)\n",
    "        c.set_axislabel('', size=0)\n",
    "    #ax.coords[1].set_major_formatter('dd:mm')\n",
    "    #ax.coords[0].set_major_formatter('mm')\n",
    "\n",
    "        \n",
    "    ax.tick_params(labelsize=9)\n",
    "    ax.set_xlabel(' ')\n",
    "    ax.set_ylabel(' ')\n",
    "    del host_background.data, host_background\n",
    "    fp.close()\n",
    "    \n",
    "ax = plt.subplot(gs[38:40])\n",
    "v = np.linspace(-275, 275, 512)\n",
    "ax.imshow(np.atleast_2d(v), aspect='auto', extent=(-275, 275, 0, 1), cmap=\"cool\", label=\"colorbar\")\n",
    "ax.tick_params(labelsize=14)\n",
    "ax.set_xlabel(r\"$\\Delta \\mathrm{HRV}\\quad[\\mathrm{km}\\,\\mathrm{s}^{-1}]$\", fontsize=16)\n",
    "ax.yaxis.set_visible(False)\n",
    "\n",
    "plt.savefig(\"/home/yymao/Downloads/host_sat_config.pdf\", bbox_inches=\"tight\", dpi=150, pad_inches=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for k, d in bg_keep.items():\n",
    "    if k == \"NGC3976\":\n",
    "        plt.plot(d[0], label=k)\n",
    "    else:\n",
    "        plt.plot(d[:,0], label=k)\n",
    "    \n",
    "plt.legend()\n",
    "plt.axhline(17500, c=\"C3\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}