{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "from urllib import parse\n",
    "from xml import etree\n",
    "import zipfile\n",
    "\n",
    "import requests\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "from astropy import units as u\n",
    "from astropy import table\n",
    "from astropy import wcs\n",
    "from astropy.coordinates import Distance\n",
    "\n",
    "from astropy.io import fits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import glue_jupyter as gj\n",
    "from glue import core as gcore"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load the Pleiades "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start by downloading the Gaia DR2 dataset for the 3 degree area around the Pleiades using `astroquery`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astroquery import gaia"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The below will either download the dataset and save it locally, or it will load the saved copy if you've already downloaded."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tab = gaia.Gaia.query_object_async('Messier 45', 90*u.arcmin)\n",
    "tab = table.QTable(tab)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now lets select a sub-set of interesting columns, and sub-sample on the part of the dataset that has a valid distance and is within a reasonable range that could be anywhere near us."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subtab = tab['designation', 'ra', 'dec', 'parallax', 'parallax_error', 'pmra', 'pmra_error', 'pmdec', 'pmdec_error', \n",
    "             'phot_g_mean_mag' ,'bp_rp']\n",
    "\n",
    "subtab['distance'] = u.Quantity(subtab['parallax']).to(u.pc, u.parallax())\n",
    "subtab['distance_error'] = u.Quantity(subtab['parallax_error']).to(u.pc, u.parallax())\n",
    "\n",
    "subtab = subtab[np.isfinite(subtab['distance'])&(u.Quantity(subtab['distance'])<300*u.pc)&(subtab['distance']>0)]\n",
    "len(subtab)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Set up Glupyter "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we set up a data object for Glupyter, and create a Glupyter app for further examination."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gaia_data = gcore.Data('Gaia_Pleiades', **subtab)\n",
    "app = gj.jglue(gaia_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scats = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "view3d = app.scatter3d(x='ra', y='dec', z='distance')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scats.append(app.scatter2d(x='distance', y='phot_g_mean_mag'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "scats.append(app.scatter2d(x='bp_rp', y='phot_g_mean_mag'))\n",
    "\n",
    "state = scats[-1].state\n",
    "state.y_max, state.y_min = state.y_min, state.y_max"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scats.append(app.scatter2d(x='pmra', y='pmdec'))\n",
    "scats[-1].scale_x.min = -100\n",
    "scats[-1].scale_x.max = 100\n",
    "scats[-1].scale_y.min = -100\n",
    "scats[-1].scale_y.max = 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scats.append(app.scatter2d(x='pmra', y='pmra_error'))\n",
    "scats[-1].state.x_min, scats[-1].state.x_max = -50, 50\n",
    "\n",
    "scats.append(app.scatter2d(x='pmdec', y='pmdec_error'))\n",
    "scats[-1].state.x_min, scats[-1].state.x_max = -100, 50"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now go ahead a select a subset using all the plots above..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if len(gaia_data.subsets) == 0:\n",
    "    print(\"You didn't select anything before running this cell, \"\n",
    "          \"so I'll select the brightest stars for you!\")\n",
    "    app.subset('Bright stars', gaia_data.id['phot_g_mean_mag'] < 8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "seltab = subtab[gaia_data.subsets[0].to_mask()]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And compute the distance to the Pleiades!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.mean(seltab['distance'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or some percentiles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.percentile(seltab['distance'], [10, 32, 50, 68, 90])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And while we're at it - lets see what the mean proper motion is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.mean(seltab['pmra']), np.mean(seltab['pmdec'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.mean(np.hypot(seltab['pmra']*np.cos(seltab['dec']), seltab['pmdec']))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Now overlay HST image "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we look at an HST image of a small part of the Pleiades, and try overplotting some of the Gaia stars"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astroquery.mast import Observations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "mast_qry = Observations.query_criteria(target_name='PLEIADESFIELD1B', obs_collection='hst')\n",
    "mast_qry"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now lets download just the final drizzled image for the observation with the bluer filter (highest resolution): F475W"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result = Observations.download_products(mast_qry['obsid'][mast_qry['filters']=='F475W'], productSubGroupDescription='DRZ')\n",
    "result['Local Path']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert len(result) == 1\n",
    "\n",
    "f475w_image = fits.open(result['Local Path'][0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we add the image data to glue and display it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sciim = app.add_data(f475w=f475w_image)[0]\n",
    "imview = app.imshow(data=sciim)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we glue together the image and the Gaia data. (Note that it might take some time for the view to catch up after you make the link.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "app.add_link(gaia_data, 'ra', sciim, 'Right Ascension')\n",
    "app.add_link(gaia_data,'dec', sciim, 'Declination')\n",
    "\n",
    "imview.add_data(gaia_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Markers should appear in the image viewer. You may need to fiddle with the color a bit to see the markers.  Since you've already made your subset selection: are any of the Pleiades members in the HST field?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You'll probably see them: they are there, but clearly offset from where the HST image has them!  This might at first appear to be a bug, but it's not.  Why might this be?"
   ]
  }
 ],
 "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}