query_dl_nicer.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# query_dl_nicer.ipynb\n",
    "This notebook:\n",
    "* queries the NICER observation tables for a source\n",
    "* saves the meta information on observations for a source to a table\n",
    "* downloads the data files from the NASA servers\n",
    "* decrypts GPG files if needed\n",
    "* applies the barycenter correction\n",
    "\n",
    "You must have heainit running in the shell you're running this from in order to run `nicerl2` and `barycorr`! If it's not running, close the notebook, run heainit, then re-open the notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import subprocess\n",
    "import requests  ## to get and read in the website for BeautifulSoup to parse\n",
    "from bs4 import BeautifulSoup  ## to scrape and parse the website html\n",
    "from astropy.table import Table, Column  ## to use astropy tables as our data storage and interaction format\n",
    "from astropy import units as u\n",
    "import urllib.request ## to download the data files via ftp\n",
    "import sys\n",
    "import numpy as np\n",
    "from datetime import datetime"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Credentials and decryption passphrase"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Team credentials, to log into the site with the segment summary table.\n",
    "cred_user = \n",
    "cred_pass = \n",
    "# Decryption passphrase for unvalidated gpg files.\n",
    "gpg_pass = "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Methods that parse the observation segment summary table website"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def scrape_site_for_seg_table():\n",
    "    \"\"\" Scrapes the NICER target segment website for the body and header \n",
    "        of the big table listing all the observations. \"\"\"\n",
    "    nicer_tab_url=\"https://heasarc.gsfc.nasa.gov/docs/nicer/team_schedule/nicer_seg_team.html\"\n",
    "    response = requests.get(nicer_tab_url, auth=(cred_user, cred_pass))  # using the NICER team credentials\n",
    "#     print(response)\n",
    "    soup = BeautifulSoup(response.text, 'lxml')\n",
    "#     print(soup)\n",
    "    t = soup.find('table', attrs={'id':'nicer_observation_segment_summary'})\n",
    "    t_header = t.find('thead')\n",
    "    t_body = t.find('tbody')\n",
    "    return t_header, t_body\n",
    "\n",
    "def get_column_labels(t_head):\n",
    "    \"\"\" Gets the text for the column headers from the table on the website, \n",
    "        to be used as the 'description' for our output table. \"\"\"\n",
    "    label_list = []\n",
    "    h_rows = t_head.find_all('th')\n",
    "    for row in h_rows:\n",
    "        label=row.text\n",
    "        label_list.append(label)\n",
    "    return label_list\n",
    "\n",
    "def get_object_observation_table(obj_name):\n",
    "    \"\"\" Gets and makes a table of all the observations for a specified object, \n",
    "        and returns the information as an astropy table. \"\"\"\n",
    "    tab_header, tab_body = scrape_site_for_seg_table()\n",
    "    col_labels = get_column_labels(tab_header)\n",
    "\n",
    "    obsID_list = []\n",
    "    start_time_list = []\n",
    "    stop_time_list = []\n",
    "    ra_list = []\n",
    "    dec_list = []\n",
    "    expo_ontarg_list = []\n",
    "    expo_good_list = []\n",
    "    proc_date_list = []\n",
    "    proc_state_list = []\n",
    "    proc_status_list = []\n",
    "    proc_ver_list = []\n",
    "    dir_link_list = []\n",
    "\n",
    "    # rows = tab_body.find_all('tr')[0:20]\n",
    "    rows = tab_body.find_all('tr')\n",
    "    for row in rows:\n",
    "        cells = row.select(\"td\")\n",
    "#         print(cells)\n",
    "        target_id = int(cells[3].string)\n",
    "        if target_id != 0:  # the Null target name, where there's no data\n",
    "            target_name = cells[4].string\n",
    "            if target_name == obj_name: \n",
    "                obsID = int(cells[0].string)\n",
    "                start_time = datetime.strptime(cells[1].string, \"%Y-%m-%dT%H:%M:%S\")\n",
    "                stop_time = datetime.strptime(cells[2].string, \"%Y-%m-%dT%H:%M:%S\")\n",
    "                ra = float(cells[5].string)\n",
    "                dec = float(cells[6].string)\n",
    "                expo_ontarg = cells[7].string\n",
    "                expo_good = cells[8].string\n",
    "                proc_date = datetime.strptime(cells[9].string, \"%Y-%m-%dT%H:%M:%S\")\n",
    "                proc_state = cells[10].string\n",
    "                proc_status = cells[11].string\n",
    "                proc_ver = cells[12].string\n",
    "                archive_links = cells[13].find_all('a', href=True)\n",
    "                if len(str(archive_links)) > 20:\n",
    "                    dir_link = archive_links[-1]['href']\n",
    "                    obsID_list.append(obsID)\n",
    "                    start_time_list.append(start_time)\n",
    "                    stop_time_list.append(stop_time)\n",
    "                    ra_list.append(ra)\n",
    "                    dec_list.append(dec)\n",
    "                    expo_ontarg_list.append(expo_ontarg)\n",
    "                    expo_good_list.append(expo_good)\n",
    "                    proc_date_list.append(proc_date)\n",
    "                    proc_state_list.append(proc_state)\n",
    "                    proc_status_list.append(proc_status)\n",
    "                    proc_ver_list.append(proc_ver)\n",
    "                    dir_link_list.append(dir_link)\n",
    "\n",
    "    obsID_col = Column(obsID_list, name='obs id', dtype='i', description=col_labels[0])\n",
    "    start_time_col = Column(start_time_list, name='start time', dtype=datetime, unit='UTC', description=col_labels[1])\n",
    "    stop_time_col = Column(stop_time_list, name='stop time', dtype=datetime, unit='UTC', description=col_labels[2])\n",
    "    ra_col = Column(ra_list, name='ra', dtype='f', unit=u.deg, description=col_labels[5])\n",
    "    dec_col = Column(dec_list, name='dec', dtype='f', unit=u.deg, description=col_labels[6])\n",
    "    expo_ontarg_col = Column(expo_ontarg_list, name='on-targ expo', dtype='f', unit=u.s, description=col_labels[7])\n",
    "    expo_good_col = Column(expo_good_list, name='good expo', dtype='f', unit=u.s, description=col_labels[8])\n",
    "    proc_date_col = Column(proc_date_list, name='proc date', dtype=datetime, description=col_labels[9])\n",
    "    proc_state_col = Column(proc_state_list, name='proc state', dtype='U', description=col_labels[10])\n",
    "    proc_status_col = Column(proc_status_list, name='proc status', dtype='i1', description=col_labels[11])\n",
    "    proc_ver_col = Column(proc_ver_list, name='proc ver', dtype='U', description=col_labels[12])\n",
    "    dir_link_col = Column(dir_link_list, name='dir link', dtype='U', description='Archive link, directory')\n",
    "    ## more on datatypes in python here: https://docs.scipy.org/doc/numpy/reference/arrays.dtypes.html\n",
    "\n",
    "    obj_tab = Table([obsID_col, start_time_col, stop_time_col, ra_col, \n",
    "                    dec_col, expo_ontarg_col, expo_good_col, proc_date_col, \n",
    "                    proc_state_col, proc_status_col, proc_ver_col,\n",
    "                    dir_link_col])\n",
    "    return obj_tab"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get a table of all the NICER observations for a specified source"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "obj_name = \"Swift_J1728.9-3613\"\n",
    "obj_prefix = \"SwiftJ1728\"\n",
    "obj_tab = get_object_observation_table(obj_name)\n",
    "print(obj_tab.info)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(min(obj_tab['start time']))\n",
    "print(max(obj_tab['start time']))\n",
    "# print(obj_tab['obs id','start time', 'good expo'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp_tab = obj_tab\n",
    "## Certain obsIDs\n",
    "# obsID_want = [1050390122, 1050390134]\n",
    "# obsID_mask = [x in obsID_want for x in tmp_tab['obs id']]\n",
    "# tmp_tab = tmp_tab[obsID_mask]\n",
    "## Exposure longer than 100 seconds\n",
    "exp_mask = [x > 100 for x in tmp_tab['on-targ expo']]\n",
    "tmp_tab = tmp_tab[exp_mask]\n",
    "## obsIDs greater than one (proprietary GOF ones have negative numbers)\n",
    "obsid_mask = [x > 1 for x in tmp_tab['obs id']]\n",
    "tmp_tab = tmp_tab[obsid_mask]\n",
    "## Dates earlier than Sept 1 2018\n",
    "# date_mask = [(x >= datetime(2018, 6, 1)) and (x < datetime(2018, 8, 1)) for x in tmp_tab['start time']]\n",
    "\n",
    "# date_mask = [(x >= datetime(2021, 2, 1)) for x in tmp_tab['start time']]\n",
    "# tmp_tab = tmp_tab[date_mask]\n",
    "print(tmp_tab['obs id','start time', 'on-targ expo', 'good expo'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you like what you see, download those observations!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "obj_tab = tmp_tab\n",
    "print(obj_tab)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Save obsID list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "homedir = os.path.expanduser(\"~\")\n",
    "exe_dir = os.getcwd()\n",
    "listdir = \"%s/Documents/Research/%s/in\" % (homedir, obj_prefix)\n",
    "if not os.path.exists(listdir):\n",
    "    os.makedirs(listdir)  # Recursive mkdir\n",
    "obsID_file = \"%s/%s_obsIDs.txt\" % (listdir, obj_prefix)\n",
    "np.savetxt(obsID_file, obj_tab['obs id'], fmt=\"%d\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Specify local directories for downloading the data files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_dir = homedir+\"/Reduced_data/\"+obj_name+\"/\"\n",
    "if not os.path.exists(data_dir):\n",
    "    os.makedirs(data_dir) # Recursive mkdir\n",
    "dl_script = data_dir+\"download.sh\"\n",
    "dl_log = data_dir+\"download.log\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Downloading the data files via FTP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(obj_tab['dir link'][0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(dl_script, 'w') as f:\n",
    "    f.write(\"rm %s \\n\" % dl_log)\n",
    "    for tab_row in obj_tab:\n",
    "        obsID =  str(tab_row['obs id'])\n",
    "        obsID_dir = data_dir + \"/\" + obsID\n",
    "        ## Downloading the whole thing\n",
    "        if not os.path.exists(obsID):\n",
    "            f.write(\"wget -r --directory-prefix=%s --append-output=%s -nH --cut-dirs=5 -l0 -c -N -np =R 'index*' -e robots=off --retr-symlinks --show-progress --progress=bar %s \\n\" % \\\n",
    "                     (data_dir, dl_log, tab_row['dir link']))\n",
    "\n",
    "print(\"Run these things at the command line:\")\n",
    "print(\"bash\")\n",
    "print(\"cd %s\" % data_dir)\n",
    "print(\"chmod u+x %s\" % os.path.basename(dl_script))\n",
    "print(\"./%s\" % os.path.basename(dl_script))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Running nicerl2 (if you didn't download the event_cl files) \n",
    "(nicercal, niextract-events, nimaketime, and nicer-mergeclean with min_fpm=7) to get cleaned event lists"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "red_script = data_dir+\"reduce.sh\"\n",
    "red_log = data_dir+\"reduce.log\"\n",
    "\n",
    "with open(red_script, 'w') as f:\n",
    "    f.write(\"rm %s \\n\" % red_log)\n",
    "    for obsid in obj_tab['obs id']:\n",
    "        f.write(\"nicerl2 indir=%s clobber=yes \\n\" % str(obsid)) \n",
    "print(\"Run these things at the command line:\")\n",
    "print(\"bash\")\n",
    "print(\"heainit\")\n",
    "print(\"cd %s\" % data_dir)\n",
    "print(\"chmod u+x %s\" % os.path.basename(red_script))\n",
    "print(\"./%s\" % os.path.basename(red_script))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## If files are encrypted, decrypt here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Make a list of the event files (in their local directory)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cl_list = []\n",
    "orb_list = []\n",
    "for obsid in obj_tab['obs id']:\n",
    "    cl_file = \"%s/%s/xti/event_cl/ni%s_0mpu7_cl.evt.gz\" % (data_dir, str(obsid), str(obsid))\n",
    "    orb_file = \"%s/%s/auxil/ni%s.orb.gz\" % (data_dir, str(obsid), str(obsid))\n",
    "    if os.path.exists(cl_file):\n",
    "        cl_list.append(cl_file)\n",
    "    else:\n",
    "        print(\"CL does not exist for obsid=%s\" % str(obsid))\n",
    "    if os.path.exists(orb_file):\n",
    "        orb_list.append(orb_file)\n",
    "    else:\n",
    "        print(\"ORB does not exist for obsid=%s\" % str(obsid))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## If the files were encrypted, use these two lines\n",
    "# evt_list = [filename.replace('.gpg','') for filename in evt_list]\n",
    "# orb_list = [filename.replace('.gpg','') for filename in orb_list]\n",
    "\n",
    "evt_list = [os.path.relpath(filename,start=data_dir) for filename in cl_list]\n",
    "print(evt_list)\n",
    "orb_list = [os.path.relpath(filename,start=data_dir) for filename in orb_list]\n",
    "print(orb_list)\n",
    "assert len(evt_list) == len(orb_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Save event list file names to a text file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "evt_out_file = listdir + \"/in/\" + obj_prefix + \"_evtlists.txt\"\n",
    "with open(evt_out_file, 'w') as f:\n",
    "    [f.write(\"%s\\n\" % evt_name) for evt_name in evt_list]\n",
    "print(\"Event list printed to: %s\" % evt_out_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Apply barycenter correction\n",
    "It doesn't like doing this to zipped files, so it will probably give a segmentation fault."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir(data_dir)\n",
    "bary_list = []\n",
    "bary_script = obj_name+\"_barycorr.sh\"\n",
    "print(\"chmod u+x %s\" % bary_script)\n",
    "with open(bary_script, mode='w') as out:\n",
    "    for (evt_file, orb_file) in zip(evt_list, orb_list):\n",
    "        bary_file = evt_file.replace('.evt', '_bary.evt')\n",
    "    #     barycorr infile=evt_file outfile=bary_file orbitfiles=orb_file refframe=ICRS clobber=yes\n",
    "        bary_cmd = \"barycorr infile=./%s outfile=./%s orbitfiles=./%s refframe=ICRS clobber=yes\" % (evt_file, bary_file, orb_file)\n",
    "        out.write(bary_cmd+\"\\n\")\n",
    "    #     bary_cmd = ['barycorr', 'infile=%s' % evt_file, 'outfile=%s' % bary_file, 'orbitfiles=%s' % orb_file, 'refframe=ICRS', 'clobber=yes']\n",
    "#         p = subprocess.Popen(bary_cmd, stdout=subprocess.PIPE, stderr=subprocess.STDOUT, shell=True)\n",
    "#         output = p.communicate()\n",
    "#         print(output)\n",
    "#         normal_output = output[0]\n",
    "#         error = output[1]\n",
    "#         exitCode = p.returncode\n",
    "        bary_list.append(bary_file)\n",
    "# os.chmod(bary_script, 754)\n",
    "# p = subprocess.Popen(bary_script, stdout=subprocess.PIPE, stderr=subprocess.STDOUT, shell=True)        \n",
    "print(\"./\"+bary_script)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bary_out_file = listdir + \"/in/\" + obj_prefix + \"_bary_evtlists.txt\"\n",
    "with open(bary_out_file, 'w') as f:\n",
    "    [f.write(\"%s\\n\" % bary_name) for bary_name in bary_list]\n",
    "print(\"Barycentered event list printed to: %s\" % bary_out_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "barycorr infile=ni1108030106_0mpu7_cl.evt.gz outfile=ni1108030106_0mpu7_cl_bary.evt.gz orbitfiles=ni1108030106.orb.gz refframe=ICRS clobber=yes\n",
    "\n"
   ]
  }
 ],
 "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
}