tutorials/Example-KDA_info_table.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example: Create a new KDA info table\n",
    "\n",
    "The BD wrapper contains a new prior that uses literature distance results to inform the P_far prior that helps resolve the kinematic distance ambiguity (KDA). In this notebook we show how to create a KDA info table that can be added to the already existing tables (which can be found in the `KDA_info` directory).\n",
    "\n",
    "The current KDA info tables mostly cover regions in the first Galactic quadrant. See Appendix A in Riener et al. (2020b) for a discussion about the KDA prior and currently available KDA info data sets.  \n",
    "\n",
    "In the following we give the code we used to produce a KDA info table from Table 1 of [Roman-Duval et al. 2009](https://ui.adsabs.harvard.edu/abs/2009ApJ...699.1153R/abstract) and Table 2 of [Rathborne et al. 2009](https://ui.adsabs.harvard.edu/abs/2009ApJS..182..131R/abstract) for a sample of molecular clouds from the Galactic Ring Survey."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "from astropy.table import Table, Column, join\n",
    "from astropy import units as u\n",
    "\n",
    "\n",
    "dirpath_kda_info = 'results'  # change this to 'os.path.join('..', 'KDA_info') for your own new tables'\n",
    "kda_table_name = 'Roman-Duval+09'\n",
    "\n",
    "#  check if dirpath_kda_info already exists, if not create it\n",
    "if not os.path.exists(dirpath_kda_info):\n",
    "    os.makedirs(dirpath_kda_info)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we load the original tables:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#  Table 1 from Roman-Duval+ 2009\n",
    "\n",
    "t1 = Table.read(\n",
    "    'https://cdsarc.unistra.fr/ftp/J/ApJ/699/1153/table1.dat',\n",
    "    readme='https://cdsarc.unistra.fr/ftp/J/ApJ/699/1153/ReadMe',\n",
    "    format=\"ascii.cds\")\n",
    "\n",
    "#  Table 2 from Rathborne+ 2009\n",
    "\n",
    "t2 = Table.read(\n",
    "    'https://cdsarc.unistra.fr/ftp/J/ApJS/182/131/table2.dat',\n",
    "    readme='https://cdsarc.unistra.fr/ftp/J/ApJS/182/131/ReadMe',\n",
    "    format=\"ascii.cds\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Table 1 from [Roman-Duval et al. 2009](https://ui.adsabs.harvard.edu/abs/2009ApJ...699.1153R/abstract) contains the final kinematic distances to the Galactic clouds but does not explicitly state whether the near or far distance was chosen. In the following, we use classes and functions contained in the module `BD_wrapper.kinematic_distance` to infer the chosen kinematic distance solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from BD_wrapper.kinematic_distance import KinematicDistance, infer_kinematic_distances\n",
    "\n",
    "\n",
    "#  initialize the KinematicDistance object\n",
    "kd = KinematicDistance()\n",
    "kd.bdc_version = '2.4'\n",
    "kd.initialize()\n",
    "\n",
    "#  the f_avg column indicates whether the associated velocity dispersion value was determined\n",
    "#  or measured for the individual source (0) or an average value was chosen (1). For the GRS\n",
    "#  table the dispersion along the velocity axis was determined for each cloud, so we set the\n",
    "#  entire f_avg column to 1.\n",
    "f_avg = np.zeros(len(t1))\n",
    "t1.add_column(Column(data=f_avg.astype('int'), name='f_avg'))\n",
    "\n",
    "#  the in_grs columns just indicates whether the source is within the PPV region covered by the GRS.\n",
    "#  This column can be ignored for other data sets.\n",
    "within_grs = np.ones(len(t1))\n",
    "t1.add_column(Column(data=within_grs.astype('int'), name='in_grs'))\n",
    "\n",
    "#  infer the kinematic distance solutions and the chosen KDA resolution\n",
    "#  Roman-Duval et al. (2009) assumed a distance to the Galactic center of 8.5 kpc, se we set R_0=8.5\n",
    "dist_n, dist_f, kda = infer_kinematic_distances(t1['GLON'], t1['Dist'], R_0=8.5)\n",
    "\n",
    "#  add the three new columns to Table 1 from Roman-Duval et al. (2009)\n",
    "t1.add_column(Column(data=dist_n, name='dist_kd_n', format=\"{0:.2f}\"))\n",
    "t1.add_column(Column(data=dist_f, name='dist_kd_f', format=\"{0:.2f}\"))\n",
    "t1.add_column(Column(data=kda, name='KDA'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we combine Table 1 of [Roman-Duval et al. 2009](https://ui.adsabs.harvard.edu/abs/2009ApJ...699.1153R/abstract) and Table 2 of [Rathborne et al. 2009](https://ui.adsabs.harvard.edu/abs/2009ApJS..182..131R/abstract) as we need information from both catalogues. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "table = join(t1, t2, keys=['GRSMC'], join_type='left')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can go on to produce the KDA info table. There are few things to consider here. To keep the structure to the other KDA info tables similar we suggest renaming the columns containing the names of the sources `'ID'` and the column containing the distances `'Dist_lit'`.\n",
    "\n",
    "The spatial and spectral extent of the sources should be given in terms of their FWHM extent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from BD_wrapper.KDA_tables import kda_info_table, kda_info_table_ini\n",
    "\n",
    "table.rename_column('GRSMC', 'ID')\n",
    "table.rename_column('Dist', 'Dist_lit')\n",
    "\n",
    "#  a: Major axis size of ellipse in deg {erroneously called semi-major axis in Rathborne+ 2009}\n",
    "#  b: Minor axis size of ellipse in deg {erroneously called semi-minor axis in Rathborne+ 2009}\n",
    "#  DelV: FWHM extent in velocity [km/s]\n",
    "#  pa: Position angle measured w.r.t. the positive GLON axis -> offset_pa = 0\n",
    "\n",
    "#  in the following we specify the corresponding column name in our merged table\n",
    "\n",
    "#  Galactic longitude value in [deg]\n",
    "glon = 'GLONc'\n",
    "#  Galactic latitude value in [deg]\n",
    "glat = 'GLATc'\n",
    "#  VLSR position of the source in [km/s]\n",
    "vlsr = 'Vlsr_1'\n",
    "#  semi-major axis of source, the unit can be specified with unit_a\n",
    "a = 'a'\n",
    "#  semi-minor axis of source, the unit can be specified with unit_b\n",
    "b = 'b'\n",
    "#  position angle of the source; if pa is measured east of Galactic north \n",
    "pa = 'pa'\n",
    "#  spectral extent of the source, typically given as the full width at half maximum (FWHM)\n",
    "dvlsr = 'DelV'\n",
    "#  column specifying whether the chosen literature distance correponds to the near \n",
    "#  or far kinematic distance solution (given as 'N' or 'F')\n",
    "kda = 'KDA'\n",
    "\n",
    "#  the following factors ease the homogenization for different catalogues as entries\n",
    "#  can be for example given in terms of the FWHM extent or standard deviation\n",
    "\n",
    "fwhm_factor = 2.354820045\n",
    "#  if a is given as major axis -> factor_a = 0.5\n",
    "#  if a is given as semi-major axis -> factor_a = 1.0\n",
    "factor_a = 0.5\n",
    "#  if b is given as minor axis -> factor_b = 0.5\n",
    "#  if a is given as semi-minor axis -> factor_b = 1.0\n",
    "factor_b = 0.5\n",
    "#  if dvlsr is given as the spectral FWHM extent -> factor_dvlsr = 1\n",
    "#  if dvlsr is given as velocity dispersion -> factor_dvlsr = fwhm_factor\n",
    "factor_dvlsr = 1\n",
    "#  if position angle is measured w.r.t. the positive GLON axis -> offset_pa = 0\n",
    "#  if position angle is measured east of north -> offset_pa = 90\n",
    "offset_pa = 0\n",
    "\n",
    "#  specify the units of a, b, and dvlsr \n",
    "unit_a = u.deg\n",
    "unit_b = u.deg\n",
    "unit_dvlsr = u.km/u.s\n",
    "\n",
    "#  make the KDA info table\n",
    "table = kda_info_table(\n",
    "    table, \n",
    "    glon=glon, \n",
    "    glat=glat, \n",
    "    vlsr=vlsr,\n",
    "    a=a, \n",
    "    b=b, \n",
    "    pa=pa, \n",
    "    dvlsr=dvlsr, \n",
    "    kda=kda,\n",
    "    factor_a=factor_a, \n",
    "    factor_b=factor_b, \n",
    "    factor_dvlsr=factor_dvlsr, \n",
    "    offset_pa=offset_pa,\n",
    "    unit_a=unit_a, \n",
    "    unit_b=unit_b, \n",
    "    unit_dvlsr=unit_dvlsr, \n",
    "    add_cols=['Dist_lit', 'KDA', 'f_avg', 'ID', 'in_grs'])\n",
    "\n",
    "table.write(\n",
    "    os.path.join(dirpath_kda_info, kda_table_name + '.dat'),\n",
    "    format='ascii', \n",
    "    overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we need to produce a .ini-file giving specifications for the KDA info table. See Appendix A in Riener et al. (2020b) for more information about the parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reference = 'RD+09'\n",
    "weight_cat = 0.5\n",
    "threshold_spatial = 0.9\n",
    "threshold_spectral = 0.125\n",
    "size = 'fwhm'\n",
    "linewidth = 'fwhm'\n",
    "\n",
    "text_ini = kda_info_table_ini(\n",
    "    reference, \n",
    "    weight_cat=weight_cat, \n",
    "    threshold_spatial=threshold_spatial, \n",
    "    threshold_spectral=threshold_spectral, \n",
    "    size=size, \n",
    "    linewidth=linewidth)\n",
    "with open(os.path.join(dirpath_kda_info, kda_table_name + '.ini'), 'w') as wfile:\n",
    "    wfile.write(text_ini)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [default]",
   "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}