Gaia-extinctions

Gaia sprint - extinction.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Gaia extinctions\n",
    "\n",
    "Greg Green, Doug Finkbeiner\n",
    "\n",
    "Gaia Sprint, NYC, October 2016"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we include the standard boilerplate:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from __future__ import print_function, division\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we copy+paste the function from http://argonaut.skymaps.info/usage#function-call, which allows us to query the Bayestar dust map remotely."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import json, requests\n",
    "\n",
    "def query(lon, lat, coordsys='gal', mode='full'):\n",
    "    '''\n",
    "    Send a line-of-sight reddening query to the Argonaut web server.\n",
    "    \n",
    "    Inputs:\n",
    "      lon, lat: longitude and latitude, in degrees.\n",
    "      coordsys: 'gal' for Galactic, 'equ' for Equatorial (J2000).\n",
    "      mode: 'full', 'lite' or 'sfd'\n",
    "    \n",
    "    In 'full' mode, outputs a dictionary containing, among other things:\n",
    "      'distmod':    The distance moduli that define the distance bins.\n",
    "      'best':       The best-fit (maximum proability density)\n",
    "                    line-of-sight reddening, in units of SFD-equivalent\n",
    "                    E(B-V), to each distance modulus in 'distmod.' See\n",
    "                    Schlafly & Finkbeiner (2011) for a definition of the\n",
    "                    reddening vector (use R_V = 3.1).\n",
    "      'samples':    Samples of the line-of-sight reddening, drawn from\n",
    "                    the probability density on reddening profiles.\n",
    "      'success':    1 if the query succeeded, and 0 otherwise.\n",
    "      'converged':  1 if the line-of-sight reddening fit converged, and\n",
    "                    0 otherwise.\n",
    "      'n_stars':    # of stars used to fit the line-of-sight reddening.\n",
    "      'DM_reliable_min':  Minimum reliable distance modulus in pixel.\n",
    "      'DM_reliable_max':  Maximum reliable distance modulus in pixel.\n",
    "    \n",
    "    Less information is returned in 'lite' mode, while in 'sfd' mode,\n",
    "    the Schlegel, Finkbeiner & Davis (1998) E(B-V) is returned.\n",
    "    '''\n",
    "    \n",
    "    url = 'http://argonaut.skymaps.info/gal-lb-query-light'\n",
    "    \n",
    "    payload = {'mode': mode}\n",
    "    \n",
    "    if coordsys.lower() in ['gal', 'g']:\n",
    "        payload['l'] = lon\n",
    "        payload['b'] = lat\n",
    "    elif coordsys.lower() in ['equ', 'e']:\n",
    "        payload['ra'] = lon\n",
    "        payload['dec'] = lat\n",
    "    else:\n",
    "        raise ValueError(\"coordsys '{0}' not understood.\".format(coordsys))\n",
    "    \n",
    "    headers = {'content-type': 'application/json'}\n",
    "    \n",
    "    r = requests.post(url, data=json.dumps(payload), headers=headers)\n",
    "    \n",
    "    try:\n",
    "        r.raise_for_status()\n",
    "    except requests.exceptions.HTTPError as e:\n",
    "        print('Response received from Argonaut:')\n",
    "        print(r.text)\n",
    "        raise e\n",
    "    \n",
    "    return json.loads(r.text)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we write our own function to interpolate the output of the Argonaut server in the distance axis:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def get_ebv(lon, lat, dist, coordsys='gal'):\n",
    "    \"\"\"\n",
    "    Returns samples of the Bayestar E(B-V) at a location in space, or\n",
    "    at an array of locations.\n",
    "    \n",
    "    Args:\n",
    "        lon (float or list/array of floats): Longitude (either l or RA),\n",
    "            in deg.\n",
    "        lat (float or list/array of floats): Latitude (either b or Dec),\n",
    "            in deg.\n",
    "        dist (float or list/array of floats): Distance, in pc.\n",
    "        coordsys (Optional[str]): Either 'gal' (for Galactic) or 'equ'\n",
    "            (for equatorial). Defaults to 'gal'.\n",
    "    \n",
    "    Returns:\n",
    "        Bayestar E(B-V), as a float or array of floats.\n",
    "    \"\"\"\n",
    "    \n",
    "    array_input = hasattr(lon, '__len__')\n",
    "    if not array_input:\n",
    "        lon = [lon]\n",
    "        lat = [lat]\n",
    "        dist = [dist]\n",
    "    if isinstance(lon, np.ndarray):\n",
    "        lon = lon.tolist()\n",
    "        lat = lat.tolist()\n",
    "        dist = dist.tolist()\n",
    "    \n",
    "    q = query(lon, lat, coordsys=coordsys)\n",
    "    samples = np.array(q['samples'])\n",
    "    dm_bin_edges = np.array(q['distmod'])\n",
    "    \n",
    "    dm = 5. * np.log10(dist) - 5.\n",
    "    bin_idx_ceil = np.searchsorted(dm_bin_edges, dm)\n",
    "    \n",
    "    ebv = np.zeros(samples.shape[:-1], dtype='f4')\n",
    "    \n",
    "    idx_near = (bin_idx_ceil == 0)\n",
    "    if np.any(idx_near):\n",
    "        a = 10.**(0.2 * (dm[idx_near] - dm_bin_edges[0]))\n",
    "        ebv[idx_near] = a[:,None] * samples[idx_near, :, 0]\n",
    "    \n",
    "    idx_far = (bin_idx_ceil == dm_bin_edges.size)\n",
    "    if np.any(idx_far):\n",
    "        ebv[idx_far] = samples[idx_far, :, -1]\n",
    "    \n",
    "    idx_btw = ~idx_near & ~idx_far\n",
    "    if np.any(idx_btw):\n",
    "        dm_ceil = dm_bin_edges[bin_idx_ceil[idx_btw]]\n",
    "        dm_floor = dm_bin_edges[bin_idx_ceil[idx_btw]-1]\n",
    "        a = (dm_ceil - dm[idx_btw]) / (dm_ceil - dm_floor)\n",
    "        ebv[idx_btw] = (\n",
    "            (1.-a[:,None]) * samples[idx_btw, :, bin_idx_ceil[idx_btw]]\n",
    "            +    a[:,None] * samples[idx_btw, :, bin_idx_ceil[idx_btw]-1]\n",
    "        )\n",
    "    \n",
    "    if not array_input:\n",
    "        ebv = ebv[0]\n",
    "    \n",
    "    return ebv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output is samples of the E(B-V) reddening to the given location:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.13038647,  0.13011999,  0.14299667,  0.12550372,  0.09960373,\n",
       "        0.12699078,  0.10722706,  0.11520647,  0.17113589,  0.12061216,\n",
       "        0.1029247 ,  0.14841001,  0.10138353,  0.13360509,  0.11327431,\n",
       "        0.14808059,  0.08403745,  0.12453725,  0.09821726,  0.14949039], dtype=float32)"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "get_ebv(90., 10., 700.) # l, b, distance (in pc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also query multiple locations at once:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.10133666,  0.11019059,  0.09275803,  0.11112501,  0.07423601,\n",
       "         0.09352152,  0.09788472,  0.10909561,  0.09985435,  0.09880251,\n",
       "         0.10847612,  0.11925913,  0.09274809,  0.11662331,  0.08620881,\n",
       "         0.1036787 ,  0.12739938,  0.12906626,  0.08649748,  0.12032939],\n",
       "       [ 0.07150812,  0.07159445,  0.08678027,  0.0912026 ,  0.1039279 ,\n",
       "         0.08121216,  0.0905378 ,  0.08752158,  0.09822951,  0.0979644 ,\n",
       "         0.09067247,  0.09841058,  0.1070139 ,  0.09406628,  0.09593894,\n",
       "         0.09060942,  0.12192534,  0.09859475,  0.09126664,  0.0935045 ],\n",
       "       [ 0.00386094,  0.01582209,  0.01418209,  0.01608209,  0.01594421,\n",
       "         0.01164722,  0.00893419,  0.02178   ,  0.01103047,  0.01536419,\n",
       "         0.01349   ,  0.01518047,  0.01013628,  0.01101628,  0.00959837,\n",
       "         0.01911209,  0.00968209,  0.01826209,  0.02244628,  0.01194209]], dtype=float32)"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ra = np.random.random(3) * 360.\n",
    "dec = np.random.random(3) * 120. - 30.\n",
    "d = np.random.random(3) * 2000.\n",
    "get_ebv(ra, dec, d, coordsys='equ') # Notice that we're using Equatorial coordinates here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to convert these measurements to G-band extinctions, we need to calibrate the Bayestar E(B-V) to G-band extinction relationship, which should be linear (to good approximation)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}