fiberassign/FiberAssignDECaLS.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# fiberassign on DECaLS (DR8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os, sys, subprocess, fitsio\n",
    "from collections import Counter\n",
    "\n",
    "import numpy as np\n",
    "from astropy.table import Table\n",
    "from astropy.io import fits\n",
    "\n",
    "import desimodel.io\n",
    "import desimodel.focalplane\n",
    "import desimodel.footprint\n",
    "from desitarget.targetmask import desi_mask, obsconditions\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Learning goals\n",
    "\n",
    "In this notebook you will learn how to:\n",
    "\n",
    "* Use observational data from DR8 to prepare it for fiberassign.\n",
    "* Run fiberassign using the files generated above.\n",
    "* Explore the outputs of fiberassign.\n",
    "\n",
    "\n",
    "See https://desi.lbl.gov/trac/wiki/Computing/JupyterAtNERSC for instructions on configuring jupyter kernels with pre-installed DESI software at NERSC.  This tutorial was last tested with the 22.2 kernel on April 29 2022.\n",
    "\n",
    "This notebook reuses code in the [main fiberassign notebook](https://github.com/desihub/tutorials/blob/master/FiberAssign.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# You can set this to a local path if you have the software installed and a copy of the data.\n",
    "# workdir = os.path.join(os.environ['HOME'], 'scratch', 'desi', 'tutorials', 'fiberassign_obs')\n",
    "workdir = os.path.join(os.environ['SCRATCH'], 'desi', 'test', 'fiberassign_obs')\n",
    "os.makedirs(workdir, exist_ok=True)\n",
    "os.chdir(workdir)\n",
    "outdir = os.path.join(workdir, 'output')\n",
    "os.makedirs(outdir, exist_ok=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Paths and filenames with the observational data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Change this if you are running locally.\n",
    "# paths = {\n",
    "#     \"targets\": \"/home/kisner/scratch/desi/tutorials\", \n",
    "#     \"skies\": \"/home/kisner/scratch/desi/tutorials\", \n",
    "#     \"gfas\": \"/home/kisner/scratch/desi/tutorials\",\n",
    "# }\n",
    "\n",
    "paths = {\"targets\": \"/global/cfs/cdirs/desi/target/catalogs/dr8/0.31.1/targets/main/resolve/\", \n",
    "         \"skies\": \"/global/cfs/cdirs/desi/target/catalogs/dr8/0.32.0/skies/\", \n",
    "         \"gfas\": \"/global/cfs/cdirs/desi/target/catalogs/dr8/0.32.0/gfas/\",\n",
    "}\n",
    "\n",
    "names = {\"targets\": \"dr8-hp-10,66,69,75,87,105.fits\", \n",
    "         \"skies\":\"dr8-0.32.0.fits\", \n",
    "         \"gfas\": \"dr8-0.32.0.fits\"}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some more files needed to run fiberassign"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mtlfile = os.path.join(workdir, 'mtl.fits')\n",
    "truthfile = os.path.join(workdir, 'truth.fits')\n",
    "targetcutfile = os.path.join(workdir, 'targets.fits') \n",
    "skycutfile = os.path.join(workdir, 'sky.fits') \n",
    "targetfile = os.path.join(paths[\"targets\"], \"targets-{}\".format(names[\"targets\"]))\n",
    "skyfile = os.path.join(paths[\"skies\"], \"skies-{}\".format(names[\"skies\"]))\n",
    "gfafile = os.path.join(paths[\"gfas\"], \"gfas-{}\".format(names[\"gfas\"]))\n",
    "tilefile = os.path.join(workdir, \"tiles.fits\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# tile selection\n",
    "program = \"dark\"\n",
    "\n",
    "tiles = desimodel.io.load_tiles()\n",
    "bright = tiles['PROGRAM']=='BRIGHT'\n",
    "    \n",
    "small = ((tiles['RA']>32) & (tiles['RA']<37) & (tiles['DEC']<3) & (tiles['DEC']>-3))\n",
    "\n",
    "if program==\"bright\":\n",
    "    Table(tiles[(bright)&(small)]).write(tilefile, overwrite=True)\n",
    "else:\n",
    "    Table(tiles[(~bright) & (small)]).write(tilefile, overwrite=True)\n",
    "\n",
    "print(\"Wrote tiles to {}\".format(tilefile))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# target selection results\n",
    "print('Started reading {}'.format(targetfile))\n",
    "targetdata = fitsio.read(targetfile, 'TARGETS')\n",
    "ii = (targetdata['RA']>30) &  (targetdata['RA']<39) & (targetdata['DEC']<5) & (targetdata['DEC']>-5)\n",
    "targetdata = targetdata[ii]\n",
    "print('Done reading target data')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(targetdata['RA'], targetdata['DEC'],s=0.1, alpha=0.1)\n",
    "plt.xlabel('RA [deg]')\n",
    "plt.ylabel('DEC [deg]')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# write down the targets to disk\n",
    "\n",
    "if (not os.path.exists(targetcutfile)):\n",
    "    Table(targetdata).write(targetcutfile, overwrite=True)\n",
    "    print('Done writing target cut data')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# compute sky file\n",
    "if (not os.path.exists(skycutfile)):\n",
    "    skydata = fitsio.read(skyfile)\n",
    "    ii = (skydata['RA']>30) &  (skydata['RA']<39) & (skydata['DEC']<5) & (skydata['DEC']>-5)\n",
    "    jj = (skydata['DESI_TARGET'] & desi_mask.SKY)!=0\n",
    "    skydata = skydata[ii&jj]\n",
    "    Table(skydata).write(skycutfile, overwrite=True)\n",
    "    print('Done writing sky cut data')\n",
    "skydata = Table.read(skycutfile)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(skydata['RA'], skydata['DEC'],s=0.1, alpha=0.1)\n",
    "plt.xlabel('RA [deg]')\n",
    "plt.ylabel('DEC [deg]')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#compute MTL\n",
    "if not os.path.exists(mtlfile):\n",
    "    print('computing mtl')\n",
    "    import desitarget.mtl\n",
    "    mtl = desitarget.mtl.make_mtl(targetdata, 'DARK|GRAY')\n",
    "\n",
    "    mtl.meta['EXTNAME'] = 'MTL'\n",
    "    mtl.write(mtlfile)\n",
    "    \n",
    "    #print some stats\n",
    "    print('MWS_TARGETS: {}'.format(np.count_nonzero(mtl['MWS_TARGET']!=0)))\n",
    "    print('BGS_TARGETS: {}'.format(np.count_nonzero(mtl['BGS_TARGET']!=0)))\n",
    "    print('DESI_TARGETS: {}'.format(np.count_nonzero(mtl['DESI_TARGET']!=0)))\n",
    "    print('finished computing mtl')\n",
    "mtl = Table.read(mtlfile)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(mtl['RA'], mtl['DEC'], s=0.1, alpha=0.1)\n",
    "plt.xlabel('RA [deg]')\n",
    "plt.ylabel('DEC [deg]')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By default, fiberassign uses the real focalplane geometry and properties at the current time.  During operations, if a fiber broke yesterday and you run fiberassign today with an up-to-date desimodel data checkout, then that fiber will not be assigned.  For this tutorial, we will run fiberassign with a *rundate* set to a time in the past before the start of commissioning.  This will give us a nominal focalplane layout with all positioners working."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assign_date = \"2022-01-01T00:00:00+00:00\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cmd = 'fiberassign --overwrite --mtl mtl.fits --sky sky.fits'\n",
    "cmd += ' --rundate {}'.format(assign_date)\n",
    "cmd += ' --footprint {}'.format(tilefile)\n",
    "cmd += ' --outdir {}'.format(outdir)\n",
    "\n",
    "print(cmd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('RUNNING: '+cmd)\n",
    "try:\n",
    "    results = subprocess.check_output(cmd.split(), stderr=subprocess.STDOUT)\n",
    "    print(results.decode())\n",
    "    print('--- SUCCESS ---')\n",
    "except subprocess.CalledProcessError as ex:\n",
    "    print('--- ERROR {} ---'.format(ex.returncode))\n",
    "    print(ex.output.decode())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tiles = Table.read(\"tiles.fits\")\n",
    "targets = Table.read(\"mtl.fits\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Gather all outputs\n",
    "assignments = list()\n",
    "for tileid in tiles['TILEID']:\n",
    "    outtile = '{}/fiberassign-{:06d}.fits'.format(outdir, tileid)\n",
    "    assignments.append(Table.read(outtile, hdu=1))\n",
    "    \n",
    "assigned_targetids = np.concatenate([tmp['TARGETID'] for tmp in assignments])\n",
    "isAssigned = np.in1d(targets['TARGETID'], assigned_targetids)\n",
    "\n",
    "plt.figure(figsize=(12,8))\n",
    "plt.plot(targets['RA'][isAssigned], targets['DEC'][isAssigned], 'k,')\n",
    "plt.title('Targets assigned to fibers')\n",
    "plt.xlabel('RA [deg]')\n",
    "plt.ylabel('DEC [deg]')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12,8))\n",
    "plt.plot(targets['RA'][~isAssigned], targets['DEC'][~isAssigned], 'k,')\n",
    "plt.title('Targets left unassigned to fibers')\n",
    "plt.xlabel('RA [deg]')\n",
    "plt.ylabel('DEC [deg]')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now run some basic QA on these outputs.  We run this on the merged outputs so that we have access to more target information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cmd = 'fba_run_qa --prefix fiberassign-'\n",
    "cmd += ' --rundate {}'.format(assign_date)\n",
    "cmd += ' --footprint {}'.format(tilefile)\n",
    "cmd += ' --dir {}'.format(outdir)\n",
    "\n",
    "print(cmd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('RUNNING: '+cmd)\n",
    "try:\n",
    "    results = subprocess.check_output(cmd.split(), stderr=subprocess.STDOUT)\n",
    "    print(results.decode())\n",
    "    print('--- SUCCESS ---')\n",
    "except subprocess.CalledProcessError as ex:\n",
    "    print('--- ERROR {} ---'.format(ex.returncode))\n",
    "    print(ex.output.decode())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the results back in and print\n",
    "import json\n",
    "from pprint import PrettyPrinter\n",
    "\n",
    "qa_file = os.path.join(outdir, \"qa.json\")\n",
    "qa_data = None\n",
    "with open(qa_file, \"r\") as f:\n",
    "    qa_data = json.load(f)\n",
    "\n",
    "pp = PrettyPrinter(indent=2)\n",
    "pp.pprint(qa_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "DESI 22.2",
   "language": "python",
   "name": "desi-22.2"
  },
  "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": 4
}