database/spectroscopic-production-database.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Spectroscopic Production Database Tutorial\n",
    "\n",
    "## Abstract\n",
    "\n",
    "This tutorial will cover the basics of using the spectroscopic production database, which is loaded from the outputs of the DESI pipeline.  Currently, this is based on software \"release\" `23.1`, and uses a [PostgreSQL](https://www.postgresql.org/) database. We use [SQLAlchemy](http://www.sqlalchemy.org/) to abstract away the details of the database.\n",
    "\n",
    "## Requirements\n",
    "\n",
    "This tutorial uses data from the `fuji` production (`/global/cfs/cdirs/desi/public/edr/spectro/redux/fuji`), and the **DESI EDR 23.1** kernel.\n",
    "\n",
    "## Getting Help\n",
    "\n",
    "If you find a problem using this notebook, please [fill out an new issue on GitHub describing your situation](https://github.com/desihub/tutorials/issues/new/).\n",
    "\n",
    "Additional details can be found in the [DESI data documentation](https://data.desi.lbl.gov/doc/access/database/).\n",
    "\n",
    "## Initial Setup\n",
    "\n",
    "This just imports everything we need and sets up paths and environment variables so we can find things."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "INFO:iers.py:82:freeze_iers: Freezing IERS table used by astropy time, coordinates.\n",
      "INFO:ephem.py:80:get_ephem: Restored ephemerides for (2019-01-01,2027-12-31) from /global/cfs/cdirs/desi/public/epo/example_files/ephem_2019-01-01_2027-12-31.fits.\n",
      "sqlalchemy==1.4.46\n",
      "specprodDB==1.1.0\n",
      "desisurvey==0.19.0\n"
     ]
    }
   ],
   "source": [
    "#\n",
    "# Imports\n",
    "#\n",
    "import os\n",
    "import sys\n",
    "import itertools\n",
    "from types import MethodType\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.font_manager import fontManager, FontProperties\n",
    "from sqlalchemy import __version__ as sqlalchemy_version\n",
    "from sqlalchemy import inspect, and_\n",
    "from sqlalchemy.sql import func\n",
    "from sqlalchemy.orm import aliased\n",
    "import astropy.units as u\n",
    "#\n",
    "# DESI software\n",
    "#\n",
    "from desiutil.log import get_logger, DEBUG\n",
    "from desitarget.targetmask import (desi_mask, mws_mask, bgs_mask)\n",
    "# from desisim.spec_qa import redshifts as dsq_z\n",
    "from desisurvey import __version__ as desisurvey_version\n",
    "from desisurvey.ephem import get_ephem, get_object_interpolator\n",
    "from desisurvey.utils import get_observer\n",
    "from specprodDB import __version__ as specprodDB_version\n",
    "import specprodDB.load as db\n",
    "#\n",
    "# Set the spectroscopic production run.\n",
    "#\n",
    "specprod = os.environ['SPECPROD'] = 'fuji'  # Change this to 'guadalupe' if needed.\n",
    "#\n",
    "# Initialize ephemerides, to find Moon, etc.\n",
    "#\n",
    "os.environ['DESISURVEY_OUTPUT'] = '/global/cfs/cdirs/desi/public/epo/example_files'\n",
    "ephem = get_ephem()\n",
    "#\n",
    "# get_ephem() will run freeze_iers(), so we import these after that.\n",
    "#\n",
    "from astropy.time import Time\n",
    "from astropy.coordinates import ICRS\n",
    "#\n",
    "# Working directory.\n",
    "#\n",
    "workingdir = os.getcwd()\n",
    "print(f'sqlalchemy=={sqlalchemy_version}')\n",
    "print(f'specprodDB=={specprodDB_version}')\n",
    "print(f'desisurvey=={desisurvey_version}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This function will compute various Moon paramters needed below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def moon(self, mjd, ra, dec):\n",
    "    \"\"\"Compute relative location of the Moon.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    mjd : float\n",
    "        Time of observation\n",
    "    ra : float\n",
    "        Right Ascension\n",
    "    dec : float\n",
    "        Declination\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    tuple\n",
    "        Moon separation, Moon altitude, Moon illumination fraction\n",
    "    \"\"\"\n",
    "    observation_time = Time(mjd, format='mjd')\n",
    "    position = ICRS(ra=ra*u.deg, dec=dec*u.deg)\n",
    "    zenith = get_observer(observation_time, alt=90 * u.deg, az=0 * u.deg).transform_to(ICRS)\n",
    "    alt = 90 * u.deg - position.separation(zenith)\n",
    "    moon_dec, moon_ra = get_object_interpolator(self.get_night(observation_time), 'moon', altaz=False)(observation_time.mjd)\n",
    "    moon_position = ICRS(ra=moon_ra*u.deg, dec=moon_dec*u.deg)\n",
    "    moon_sep = position.separation(moon_position).to(u.deg).value\n",
    "    moon_alt = (90 * u.deg - moon_position.separation(zenith)).to(u.deg).value\n",
    "    moon_frac = ephem.get_moon_illuminated_fraction(observation_time.mjd).tolist()\n",
    "    return (moon_sep, moon_alt, moon_frac)\n",
    "\n",
    "ephem.moon = MethodType(moon, ephem)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contents of the Database\n",
    "\n",
    "### Schema\n",
    "\n",
    "All tables are grouped into a database *schema* and that schema is named for the production run, (*e.g.* `fuji`).  When writing \"raw\" SQL, table names need to be schema-qualified, for example, `fuji.target`.  However, the SQLAlchemy abstraction layer is designed to take care of this for you.\n",
    "\n",
    "### Important notes\n",
    "\n",
    "* This database does not contain any sky spectra. Both deliberately-targeted (`targetid & 2**59 != 0`) and negative targetid (`targetid < 0`) sky spectra are excluded.\n",
    "* Only quantities derived from *cumulative* tile-based spectra are included at the present time.\n",
    "\n",
    "### The tables\n",
    "\n",
    "* `photometry`. This contains the pure photometric data. Usually this is derived from the LS DR9 Tractor data, but not every *targeted* object has Tractor photometry.\n",
    "  - Loaded from `tractorphot` files in the [lsdr9-photometry VAC](https://data.desi.lbl.gov/doc/releases/edr/vac/lsdr9-photometry/): `/global/cfs/cdirs/desi/public/edr/vac/edr/lsdr9-photometry/fuji/v2.1/potential-targets/tractorphot/tractorphot-potential-*-fuji.fits`.\n",
    "  - SQLAlchemy object: `db.Photometry`.\n",
    "  - Primary key: `targetid`.\n",
    "  - [Column descriptions](https://data.desi.lbl.gov/doc/access/database/#photometry).\n",
    "* `target`. This contains the targeting bits and other data generated by `desitarget`.\n",
    "  - Loaded from `targetphot` file in the [lsdr9-photometry VAC](https://data.desi.lbl.gov/doc/releases/edr/vac/lsdr9-photometry/): `/global/cfs/cdirs/desi/public/edr/vac/edr/lsdr9-photometry/fuji/v2.1/potential-targets/targetphot-potential-fuji.fits`.\n",
    "  - SQLAlchemy object: `db.Target`.\n",
    "  - Unique identifier: (`targetid`, `survey`, `tileid`).\n",
    "  - Primary key: `id`, a unique, arbitrary integer composed from (`targetid`, `survey`, `tileid`).\n",
    "  - [Column descriptions](https://data.desi.lbl.gov/doc/access/database/#target).\n",
    "* `tile`. This contains information about observations grouped by tile.\n",
    "  - Loaded from top-level `tiles-fuji.fits`.\n",
    "  - SQLAlchemy object: `db.Tile`.\n",
    "  - Primary key: `tileid`.\n",
    "  - [Column descriptions](https://data.desi.lbl.gov/doc/access/database/#tile).\n",
    "* `exposure`. This contains information about individual exposures.\n",
    "  - Loaded from top-level `exposures-fuji.fits`, `EXPOSURES` HDU.\n",
    "  - SQLAlchemy object: `db.Exposure`.\n",
    "  - Primary key: `expid`.\n",
    "  - [Column descriptions](https://data.desi.lbl.gov/doc/access/database/#exposure).\n",
    "* `frame`. This contains information about individual exposures, but broken down by camera.  There will usually, but not always, be 30 frames per exposure.\n",
    "  - Loaded from top-level `exposures-fuji.fits`, `FRAMES` HDU.\n",
    "  - SQLAlchemy object: `db.Frame`.\n",
    "  - Unique identifier: (`expid`, `camera`).\n",
    "  - Primary key: `frameid`, composed from `expid` and a mapping of `camera` to an arbitrary integer.\n",
    "  - [Column descriptions](https://data.desi.lbl.gov/doc/access/database/#frame).\n",
    "* `fiberassign`. This contains information about fiber positions.\n",
    "  - Loaded from fiberassign files in the [tiles product](https://data.desi.lbl.gov/public/edr/target/fiberassign/tiles/tags/0.5/).  All fiberassign files corresponding to tiles in the `tile` table are loaded.\n",
    "  - SQLAlchemy object: `db.Fiberassign`.\n",
    "  - Unique identifier: (`tileid`, `targetid`, `location`).\n",
    "  - Primary key: `id`, a unique, arbitrary integer composed from (`tileid`, `targetid`, `location`).\n",
    "  - [Column descriptions](https://data.desi.lbl.gov/doc/access/database/#fiberassign).\n",
    "* `potential`. This contains a list of `targetid`s that *could* have been targeted on a given tile.\n",
    "  - Loaded from the `POTENTIAL_ASSIGNMENTS` HDU in the same fiberassign files mentioned above.\n",
    "  - SQLAlchemy object: `db.Potential`.\n",
    "  - Unique identifier: (`tileid`, `targetid`, `location`).\n",
    "  - Primary key: `id`, a unique, arbitrary integer composed from (`tileid`, `targetid`, `location`).\n",
    "  - [Column descriptions](https://data.desi.lbl.gov/doc/access/database/#potential).\n",
    "* `zpix`. This contains the pipeline redshifts grouped by HEALPixel.\n",
    "  - Loaded from the [zcat VAC](https://data.desi.lbl.gov/doc/releases/edr/vac/zcat/): `/global/cfs/cdirs/desi/public/edr/vac/edr/zcat/fuji/v1.0/zall-pix-edr-vac.fits`.\n",
    "  - Based on the `zall-pix-fuji.fits` file in the `zcatalog/` directory.\n",
    "  - SQLAlchemy object: `db.Zpix`.\n",
    "  - Unique identifier: (`targetid`, `survey`, `program`).\n",
    "  - Primary key: `id`, a unique, arbitrary integer composed from (`targetid`, `survey`, `program`).\n",
    "  - [Column descriptions](https://data.desi.lbl.gov/doc/access/database/#zpix).\n",
    "* `ztile`. This contains the pipeline redshifts grouped by tile in a variety of ways.\n",
    "  - Loaded from the [zcat VAC](https://data.desi.lbl.gov/doc/releases/edr/vac/zcat/): `/global/cfs/cdirs/desi/public/edr/vac/edr/zcat/fuji/v1.0/zall-tilecumulative-edr-vac.fits`.\n",
    "  - Based on the `zall-tilecumulative-fuji.fits` file in the `zcatalog/` directory.\n",
    "  - SQLAlchemy object: `db.Ztile`.\n",
    "  - Unique identifier: (`targetid`, `spgrp`, `spgrpval`, `tileid`).\n",
    "  - Primary key: `id`, a unique, arbitrary integer composed from (`targetid`, `spgrp`, `spgrpval`, `tileid`).\n",
    "  - [Column descriptions](https://data.desi.lbl.gov/doc/access/database/#ztile).\n",
    "* `version`. This contains metadata related to loading the database schema associated with the production run. It has no relationships to other tables.\n",
    "  - SQLAlchemy object: `db.Version`.\n",
    "  - Primary key: `id`, a sequential integer.\n",
    "\n",
    "### Foreign key relationships\n",
    "\n",
    "* `fuji.target.targetid` -> `fuji.photometry.targetid`.\n",
    "* `fuji.target.tileid` -> `fuji.tile.tileid`.\n",
    "* `fuji.exposure.tileid` -> `fuji.tile.tileid`.\n",
    "* `fuji.frame.expid` -> `fuji.exposure.expid`.\n",
    "* `fuji.fiberassign.targetid` -> `fuji.photometry.targetid`.\n",
    "* `fuji.fiberassign.tileid` -> `fuji.tile.tileid`.\n",
    "* `fuji.potential.targetid` -> `fuji.photometry.targetid`.\n",
    "* `fuji.potential.tileid` -> `fuji.tile.tileid`.\n",
    "* `fuji.zpix.targetid` -> `fuji.photometry.targetid`.\n",
    "* `fuji.ztile.targetphotid` -> `fuji.target.id`.\n",
    "* `fuji.ztile.targetid` -> `fuji.photometry.targetid`.\n",
    "* `fuji.ztile.tileid` -> `fuji.tile.tileid`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initial Database Connection\n",
    "\n",
    "This connection uses a `~/.pgpass` file to store connection credentials.\n",
    "The following code will attempt to connect to the database.  If it fails, it prints the\n",
    "commands to update your `~/.pgpass` file (but it won't do that automatically in case you\n",
    "did a Jupyter \"run all cells\")."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# For much more output, use DEBUG/verbose mode in the commands below, e.g:\n",
    "#\n",
    "# db.log = get_logger(DEBUG)\n",
    "# postgresql = db.setup_db(schema=specprod, hostname='specprod-db.desi.lbl.gov', username='desi_public', verbose=True)\n",
    "#\n",
    "db.log = get_logger()\n",
    "try:\n",
    "    postgresql = db.setup_db(schema=specprod, hostname='specprod-db.desi.lbl.gov', username='desi_public')\n",
    "except RuntimeError:\n",
    "    print(\"\\nIt looks like you need to add access credentials to your $HOME/.pgpass file.\")\n",
    "    print(\"Try running this from the command line and then rerun this cell.\")\n",
    "    print()\n",
    "    print(\"cat /global/common/software/desi/desi_public.pgpass >> ~/.pgpass; chmod 600 ~/.pgpass\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Learning About the Tables\n",
    "\n",
    "The tables in the database are listed above.  To inspect an individual table, you can use the `__table__` attribute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# Print the table columns and their types.\n",
    "#\n",
    "[(c.name, c.type) for c in db.Zpix.__table__.columns]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also `inspect()` the database.  For details see [here](http://docs.sqlalchemy.org/en/latest/core/inspection.html?highlight=inspect#module-sqlalchemy.inspection)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "inspector = inspect(db.engine)\n",
    "for table_name in inspector.get_table_names(schema=specprod):\n",
    "    print(table_name)\n",
    "    for column in inspector.get_columns(table_name, schema=specprod):\n",
    "        print(\"Column: {name} {type}\".format(**column))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercises\n",
    "\n",
    "* What is the type of the `night` column of the `exposures` table?\n",
    "* What is the primary key of the `ztile` table?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Simple Queries\n",
    "\n",
    "Queries are set up with the `.query()` method on Session objects.  In this case, there's a prepared Session object called `db.dbSession`.  `.filter()` corresponds to a `WHERE` clause in SQL.\n",
    "\n",
    "In most of the examples below, we include the equivalent raw SQL command that corresponds to the query."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Versions\n",
    "\n",
    "Each schema contains a `version` table that describes the software and other tagged data sets used to load that schema. This is mainly intended for internal metadata tracking purposes, but it is still visible.\n",
    "\n",
    "```SQL\n",
    "SELECT * FROM fuji.version;\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q = db.dbSession.query(db.Version).all()\n",
    "q"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exposures, Nights, Tiles\n",
    "\n",
    "Here are some simple queries that demonstrate simple connections between nights, exposures and tiles.\n",
    "\n",
    "#### How many tiles are there?\n",
    "\n",
    "```SQL\n",
    "SELECT COUNT(tileid) FROM fuji.tile;\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q = db.dbSession.query(db.Tile).count()\n",
    "q"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### On which nights were a particular tile observed?\n",
    "\n",
    "```SQL\n",
    "SELECT night, expid FROM fuji.exposure WHERE tileid = 100;\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q = db.dbSession.query(db.Exposure.night, db.Exposure.expid).filter(db.Exposure.tileid == 100).all()\n",
    "q"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Which tiles were observed on a night?\n",
    "\n",
    "```SQL\n",
    "SELECT tileid, survey, program FROM fuji.exposure WHERE night = 20210115;\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q = db.dbSession.query(db.Exposure.tileid, db.Exposure.survey, db.Exposure.program).filter(db.Exposure.night == 20210115).all()\n",
    "q"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Select ELG Targets\n",
    "\n",
    "Note the special way we obtain the bitwise AND operator (`desi_mask.ELG == 2**1`).\n",
    "\n",
    "```SQL\n",
    "SELECT * from fuji.target WHERE (desi_target & 2) != 0;\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q = db.dbSession.query(db.Target).filter(db.Target.desi_target.op('&')(desi_mask.ELG) != 0).all()\n",
    "q[:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "[(row.targetid, row.desi_target, row.survey, row.program, row.tileid) for row in q[:10]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise\n",
    "\n",
    "* How many objects in the `zpix` table have `spectype` 'GALAXY'?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Redshift and Classification\n",
    "\n",
    "Simple query filtering on string values. Note that the slice notation `[:20]` is equivalent to `LIMIT 20` in SQL.\n",
    "\n",
    "```SQL\n",
    "SELECT spectype, subtype, z FROM fuji.zpix WHERE spectype = 'STAR' AND subtype != '' LIMIT 20;\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q = db.dbSession.query(db.Zpix.spectype, db.Zpix.subtype, db.Zpix.z).filter(db.Zpix.spectype == 'STAR').filter(db.Zpix.subtype != '')[:20]\n",
    "q"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Joining Tables\n",
    "\n",
    "### A Simple Join\n",
    "\n",
    "Let's look at the nights and exposures on which a particular `TARGETID` was observed.\n",
    "\n",
    "```SQL\n",
    "SELECT f.tileid, e.expid, e.night FROM fuji.fiberassign AS f JOIN fuji.exposure AS e ON f.tileid = e.tileid WHERE f.targetid = 933811403620352;\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q = db.dbSession.query(db.Fiberassign.tileid, db.Exposure.expid, db.Exposure.night).filter(db.Fiberassign.tileid == db.Exposure.tileid).filter(db.Fiberassign.targetid == 933811403620352).all()\n",
    "q"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Another Simple Join\n",
    "\n",
    "In this case, we'll look at photometric flux and measured redshift. We'll `ORDER`(`BY`) the results and `LIMIT` the query with slice notation.\n",
    "\n",
    "```SQL\n",
    "SELECT p.*, z.* FROM fuji.photometry AS p JOIN fuji.zpix AS z ON p.targetid = z.targetid ORDER BY z.z, p.flux_g LIMIT 50;\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q = db.dbSession.query(db.Photometry, db.Zpix).filter(db.Photometry.targetid == db.Zpix.targetid).order_by(db.Zpix.z, db.Photometry.flux_g)[:50]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "[(row.Photometry.flux_g, row.Photometry.flux_r, row.Photometry.flux_z, row.Zpix.z) for row in q]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "flux_g = np.array([row.Photometry.flux_g for row in q])\n",
    "flux_r = np.array([row.Photometry.flux_r for row in q])\n",
    "flux_z = np.array([row.Photometry.flux_z for row in q])\n",
    "g_minus_r = np.log10(flux_r/flux_g)\n",
    "r_minus_z = np.log10(flux_z/flux_r)\n",
    "redshift = np.array([row.Zpix.z for row in q])\n",
    "fig, axes = plt.subplots(1, 1, figsize=(5, 5), dpi=100)\n",
    "p = axes.plot(g_minus_r, r_minus_z, 'k.')\n",
    "foo = axes.set_xlim([-0.2, 1.0])\n",
    "foo = axes.set_ylim([-0.2, 1.0])\n",
    "foo = axes.set_aspect('equal')\n",
    "foo = axes.set_xlabel('$g - r$')\n",
    "foo = axes.set_ylabel('$r - z$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise\n",
    "\n",
    "* Create a color-color plot for objects targeted as QSOs, and spectroscopically confirmed as such."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A More Complicated Join\n",
    "\n",
    "Let's look at objects that appear on more than one tile. For each of those tiles, how many exposures where there?\n",
    "\n",
    "In this example, we're using `sqlalchemy.sql.func` to get the equivalent of `COUNT(*)` and a subquery that itself is a multi-table join.\n",
    "\n",
    "```SQL\n",
    "SELECT t.nexp, f.tileid, q1.targetid, q1.n_assign FROM fuji.fiberassign AS f\n",
    "    JOIN (SELECT ff.targetid, COUNT(*) AS n_assign FROM fuji.fiberassign AS ff GROUP BY ff.targetid) AS q1 ON f.targetid = q1.targetid\n",
    "    JOIN fuji.tile AS t ON f.tileid = t.tileid LIMIT 100;\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# db.dbSession.rollback()\n",
    "q1 = db.dbSession.query(db.Fiberassign.targetid, func.count('*').label('n_assign')).group_by(db.Fiberassign.targetid).subquery()\n",
    "q2 = db.dbSession.query(db.Tile.nexp, db.Fiberassign.tileid, q1.c.targetid, q1.c.n_assign).filter(q1.c.n_assign>2).filter(db.Fiberassign.targetid == q1.c.targetid).filter(db.Tile.tileid == db.Fiberassign.tileid)[:100]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's see what the redshift table thinks are the number of exposures for these objects.\n",
    "\n",
    "```SQL\n",
    "SELECT z.* FROM fuji.zpix AS z WHERE z.targetid IN (6432023904256, 6448025174016[, ...]);\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q3 = db.dbSession.query(db.Zpix).filter(db.Zpix.targetid.in_([row[2] for row in q2])).all()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "[(row.coadd_numexp, row.zcat_nspec) for row in q3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise\n",
    "\n",
    "* What is the distribution of number of exposures?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Connecting Redshifts to Targeting\n",
    "\n",
    "Matching redshift results back to targeting is an important task for many DESI science projects. Below we briefly describe the fundamentals of this matching using the database.\n",
    "\n",
    "### Tile-based\n",
    "\n",
    "The `ztile` table contains a foreign key pointing to the primary key of the target table. This means that for any row of `ztile`, there is a pre-matched entry in the `target` table, and it can be accessed with a simple join, such as:\n",
    "```sql\n",
    "SELECT z.targetid, z.spgrp, z.spgrpval, z.tileid, z.z, z.zwarn, z.spectype,\n",
    "    t.sv1_desi_target, t.sv1_bgs_target, t.sv1_mws_target, t.sv1_scnd_target\n",
    "    p.ra, p.dec\n",
    "FROM fuji.ztile AS z\n",
    "JOIN fuji.target AS t ON z.targetphotid = t.id\n",
    "JOIN fuji.photometry AS p ON z.targetid = p.targetid\n",
    "LIMIT 50;\n",
    "```\n",
    "Note the `z.targetphotid = t.id`. This query also demonstrates a `JOIN` with the photometry table, where `targetid` is unique."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q = db.dbSession.query(db.Ztile.targetid, db.Ztile.spgrp, db.Ztile.spgrpval, db.Ztile.tileid, db.Ztile.z, db.Ztile.zwarn, db.Ztile.spectype,\n",
    "                       db.Target.sv1_desi_target, db.Target.sv1_bgs_target, db.Target.sv1_mws_target, db.Target.sv1_scnd_target,\n",
    "                       db.Photometry.ra, db.Photometry.dec).join(db.Target, db.Ztile.targetphotid == db.Target.id).join(db.Photometry)\n",
    "q[:50]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### HEALPixel-based\n",
    "\n",
    "The redshift results from HEALPixel-based coadds in the `zpix` table may be made up of several tiles, so it is not as simple to match these with targeting data. In the majority of cases, the targeting information is the same for a `targetid` across multiple tiles, as long as the tiles are part of the same survey (*e.g.* sv1). However it is not *always* the case, so we have instead collected the best values of the targeting information and placed those values directly in the `zpix` table. Here is the analogous query to the query above:\n",
    "```sql\n",
    "SELECT z.targetid, z.survey, z.program, z.healpix, z.z, z.zwarn, z.spectype,\n",
    "    z.sv1_desi_target, z.sv1_bgs_target, z.sv1_mws_target, z.sv1_scnd_target\n",
    "    p.ra, p.dec\n",
    "FROM fuji.zpix AS z\n",
    "JOIN fuji.photometry AS p ON z.targetid = p.targetid\n",
    "LIMIT 50;\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q = db.dbSession.query(db.Zpix.targetid, db.Zpix.survey, db.Zpix.program, db.Zpix.healpix, db.Zpix.z, db.Zpix.zwarn, db.Zpix.spectype,\n",
    "                       db.Zpix.sv1_desi_target, db.Zpix.sv1_bgs_target, db.Zpix.sv1_mws_target, db.Zpix.sv1_scnd_target,\n",
    "                       db.Photometry.ra, db.Photometry.dec).join(db.Photometry)\n",
    "q[:50]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "* Determine which tiles and exposures were involved in the determination of a HEALPixel-based redshift. *Hint*: the `fiberassign` table can connect `targetid` to `tileid`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fly me to the Moon\n",
    "\n",
    "How does the Moon affect redshifts?  First, let's find exposures that exposures that had the Moon above the horizon."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "moon_up = [e.expid for e in db.dbSession.query(db.Exposure).all() if ephem.moon(e.mjd, e.tilera, e.tiledec)[1] > 0]\n",
    "len(moon_up)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "moon_up"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So there are a few.  But there is a subtle issue: redshifts are based on *all* exposures, and the exposures are deliberately designed to enforce the bright/dark dichotomy in targeting. There are execptions though: certain LRGs also get targeted in the BGS & MWS, so that's not hard to capture. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "various_lrgs = (desi_mask.LRG | desi_mask.BGS_ANY | desi_mask.MWS_ANY)\n",
    "various_lrgs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```SQL\n",
    "SELECT z.targetid, z.z, z.zerr, z.zwarn\n",
    "    FROM fuji.ztile AS z\n",
    "    JOIN fuji.target AS t ON z.targetphotid = t.id\n",
    "    JOIN fuji.fiberassign AS f ON t.targetid = f.targetid\n",
    "    JOIN fuji.exposure AS e ON f.tileid == e.tileid\n",
    "    WHERE z.spgrp = 'cumulative'\n",
    "    AND t.desi_target & 3458764513820540929 != 0\n",
    "    AND e.expid IN (90250, 87505, 87382[...]);\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q_up = db.dbSession.query(db.Ztile.targetid, db.Ztile.z, db.Ztile.zerr, db.Ztile.zwarn)\\\n",
    "                   .join(db.Target, db.Target.id==db.Ztile.targetphotid)\\\n",
    "                   .join(db.Fiberassign, db.Target.targetid==db.Fiberassign.targetid)\\\n",
    "                   .join(db.Exposure, db.Fiberassign.tileid==db.Exposure.tileid)\\\n",
    "                   .filter(db.Ztile.spgrp=='cumulative')\\\n",
    "                   .filter(db.Target.desi_target.op('&')(various_lrgs) != 0)\\\n",
    "                   .filter(db.Exposure.expid.in_(moon_up)).all()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```SQL\n",
    "SELECT z.targetid, z.z, z.zerr, z.zwarn\n",
    "    FROM fuji.ztile AS z\n",
    "    JOIN fuji.target AS t ON z.targetphotid = t.id\n",
    "    JOIN fuji.fiberassign AS f ON t.targetid = f.targetid\n",
    "    JOIN fuji.exposure AS e ON f.tileid == e.tileid\n",
    "    WHERE z.spgrp = 'cumulative'\n",
    "    AND t.desi_target & 3458764513820540929 != 0\n",
    "    AND e.expid NOT IN (90250, 87505, 87382[...]);\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q_dn = db.dbSession.query(db.Ztile.targetid, db.Ztile.z, db.Ztile.zerr, db.Ztile.zwarn)\\\n",
    "                   .join(db.Target, db.Target.id==db.Ztile.targetphotid)\\\n",
    "                   .join(db.Fiberassign, db.Target.targetid==db.Fiberassign.targetid)\\\n",
    "                   .join(db.Exposure, db.Fiberassign.tileid==db.Exposure.tileid)\\\n",
    "                   .filter(db.Ztile.spgrp=='cumulative')\\\n",
    "                   .filter(db.Target.desi_target.op('&')(various_lrgs) != 0)\\\n",
    "                   .filter(~db.Exposure.expid.in_(moon_up)).all()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unfortunately however, the database currently only contains cumulative tile redshifts, not per-exposure redshifts, so it's not really meaningful to say whether the Moon was up or not. We'll just call this a work in progress."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q_up"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q_dn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Survey Progress\n",
    "\n",
    "Let's see which nights have data, and count the number of exposures per night.\n",
    "\n",
    "```SQL\n",
    "SELECT e.night, COUNT(e.expid) AS n_exp FROM fuji.exposure AS e GROUP BY e.night ORDER BY e.night;\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q = db.dbSession.query(db.Exposure.night, func.count(db.Exposure.expid).label('n_exp')).group_by(db.Exposure.night).order_by(db.Exposure.night).all()\n",
    "q"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Observation timestamp for a given night.  Note how we have both MJD and a corresponding `datetime.datetime` object in the database.\n",
    "\n",
    "```SQL\n",
    "SELECT e.expid, e.mjd, e.date_obs FROM fuji.exposure AS e WHERE e.night = 20210428 ORDER BY e.expid;\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q = db.dbSession.query(db.Exposure.expid, db.Exposure.mjd, db.Exposure.date_obs).filter(db.Exposure.night == 20210428).order_by(db.Exposure.expid).all()\n",
    "q"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, for a given target in the `target` table, when was the observation completed?  In other words, if a target has multiple observations, we want the date of the *last* observation.  First, how many targets are there?\n",
    "\n",
    "```SQL\n",
    "SELECT COUNT(*) FROM fuji.target;\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_targets = db.dbSession.query(db.Target).count()\n",
    "N_targets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we look for targets that have observations and find the MJD of the observation.\n",
    "\n",
    "```SQL\n",
    "SELECT f.targetid, e.expid, e.mjd FROM fuji.fiberassign AS f\n",
    "    JOIN (SELECT tt.targetid FROM fuji.target AS tt JOIN fuji.fiberassign AS ff ON tt.targetid = ff.targetid\n",
    "              JOIN fuji.exposure AS ee ON ff.tileid = ee.tileid GROUP BY tt.targetid) AS q1 ON f.targetid = q1.targetid\n",
    "    JOIN fuji.exposure AS e ON f.tileid = e.tileid ORDER BY q1.targetid, e.expid;\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# Find all targetids that have observations.\n",
    "#\n",
    "q1 = db.dbSession.query(db.Target.targetid).filter(db.Target.targetid == db.Fiberassign.targetid).filter(db.Fiberassign.tileid == db.Exposure.tileid).group_by(db.Target.targetid).subquery()\n",
    "#\n",
    "# Find the exposure times for the targetids that have been observed\n",
    "#\n",
    "q2 = db.dbSession.query(db.Fiberassign.targetid, db.Exposure.expid, db.Exposure.mjd).filter(db.Fiberassign.targetid == q1.c.targetid).filter(db.Fiberassign.tileid == db.Exposure.tileid).order_by(q1.c.targetid, db.Exposure.expid).all()\n",
    "targetid, expid, mjd = zip(*q2)\n",
    "targetid = np.array(targetid)\n",
    "expid = np.array(expid)\n",
    "mjd = np.array(mjd)\n",
    "#\n",
    "# Use the counts to give the *last* observation.\n",
    "#\n",
    "unique_targetid, i, j, c = np.unique(targetid, return_index=True, return_inverse=True, return_counts=True)\n",
    "unique_expid = expid[i + (c-1)]\n",
    "unique_mjd = mjd[i + (c-1)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have the targets and the date of last observation.  But it's sorted by `targetid`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ii = unique_expid.argsort()\n",
    "unique_targetid, i3, j3, c3 = np.unique(unique_expid[ii], return_index=True, return_inverse=True, return_counts=True)\n",
    "N_completed = np.cumsum(c3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "min_mjd = 10*(int(mjd.min())//10)\n",
    "fig, axes = plt.subplots(1, 1, figsize=(8, 8), dpi=100)\n",
    "p1 = axes.plot(unique_mjd[ii][i3] - min_mjd, N_completed/N_targets, 'k-')\n",
    "foo = axes.set_xlabel(f'MJD - {min_mjd:d}')\n",
    "foo = axes.set_ylabel('Fraction completed')\n",
    "foo = axes.grid(True)\n",
    "# foo = axes.legend(loc=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "* Break down the progress by target class, target bit, etc."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using Relationships in SQLAlchemy\n",
    "\n",
    "Here we demonstrate how table relationships can simplify certain queries. First we grab a single `db.Exposure` object.\n",
    "\n",
    "```SQL\n",
    "SELECT * FROM fuji.exposure where expid = 86507;\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "exposure = db.dbSession.query(db.Exposure).filter(db.Exposure.expid == 86507).one()\n",
    "exposure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How do we get the `db.Frame` objects associated with this exposure?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "exposure.frames"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What tile is associated with this exposure?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "exposure.tile"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What fiberassignments were made on this tile?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "exposure.tile.fiberassign[:20]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What redshifts were measured on this tile?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "exposure.tile.ztile_redshifts[:20]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "* Rewrite example queries above using relationships wherever possible."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using q3c in SQLAlchemy\n",
    "\n",
    "[q3c](https://github.com/segasai/q3c) ([Koposov & Bartunov 2006](https://ui.adsabs.harvard.edu/abs/2006ASPC..351..735K/abstract)) is a popular library that provides spatial indexing and searching in astronomical databases. Here we'll demonstrate how to access this functionality in SQLAlchemy.  Any database function is accessible with `sqlalchemy.sql.func`.  This is a radial (\"cone\") search on an arbitrary point in the DESI footprint:\n",
    "\n",
    "```SQL\n",
    "SELECT p.*, z.*, q3c_dist(p.ra, p.dec, 180.0, 0.0) AS radial_distance\n",
    "    FROM fuji.photometry AS p JOIN fuji.zpix AS z ON p.targetid = z.targetid\n",
    "    WHERE q3c_radial_query(p.ra, p.dec, 180.0, 0.0, 1.0/60.0); -- 1 arcmin\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q = db.dbSession.query(db.Photometry, db.Zpix, func.q3c_dist(db.Photometry.ra, db.Photometry.dec, 180.0, 0.0).label(\"radial_distance\")).join(db.Zpix).filter(func.q3c_radial_query(db.Photometry.ra, db.Photometry.dec, 180.0, 0.0, 1.0/60.0)).all()  # 1 arcmin\n",
    "q"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "* What spectra are near your favourite object?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "DESI 23.1",
   "language": "python",
   "name": "desi-23.1"
  },
  "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.10.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}