Validation/DC2_refcat_loader_demo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0b1f24fe",
   "metadata": {},
   "source": [
    "# DC2 Refcat Loader Demo\n",
    "\n",
    "<!--- <br>Owner: **Keith Bechtol** ([@bechtol](https://github.com/LSSTScienceCollaborations/StackClub/issues/new?body=@bechtol)) --->\n",
    "<br>Owner: **Peter Ferguson** ([@psferguson](https://github.com/LSSTScienceCollaborations/StackClub/issues/new?body=@psferguson))\n",
    "<br>Last Verified to Run: **2021-10-03**\n",
    "<br>Verified Stack Release: **w_2021_33**\n",
    "Last verified to run on 2021-08-16 with LSST Science Pipelines release w_2021_33 <br>\n",
    "Contact authors: Peter Ferguson <br>\n",
    "Target audience: All DP0 delegates. <br>\n",
    "Container Size: medium <br>\n",
    "Questions welcome at <a href=\"https://community.lsst.org/c/support/dp0\">community.lsst.org/c/support/dp0</a> <br>\n",
    "Find DP0 documentation and resources at <a href=\"https://dp0-1.lsst.io\">dp0-1.lsst.io</a> <br>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f310a9e",
   "metadata": {},
   "source": [
    "**Credit:** This tutorial was originally developed by Keith Bechtol."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f77eff54",
   "metadata": {},
   "source": [
    "### Learning Objectives\n",
    "\n",
    "This notebook demonstrates how to: <br>\n",
    "1. Determine the reference catalog associated with a dataset \n",
    "2. Load this reference catalog \n",
    "3. Load a source catalog\n",
    "4. Load the reference catalog that overlaps"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1dc0cda5",
   "metadata": {},
   "source": [
    "### Set Up \n",
    "You can find the Stack version by using `eups list -s` on the terminal command line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6aa48917",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Site, host, and stack version\n",
    "! echo $EXTERNAL_INSTANCE_URL\n",
    "! echo $HOSTNAME\n",
    "! eups list -s | grep lsst_distrib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7e9091cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os, os.path\n",
    "import numpy as np\n",
    "from astropy.time import Time\n",
    "\n",
    "import lsst.geom\n",
    "from lsst.pipe.tasks.loadReferenceCatalog import LoadReferenceCatalogConfig, LoadReferenceCatalogTask\n",
    "from lsst.meas.algorithms import ReferenceObjectLoader\n",
    "import lsst.daf.butler as dafButler\n",
    "from lsst.utils import getPackageDir\n",
    "\n",
    "from astropy.table import vstack\n",
    "import astropy.units as u\n",
    "import astropy.coordinates as coord\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "061aca1f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up some plotting defaults:\n",
    "\n",
    "params = {\n",
    "   'axes.labelsize': 28,\n",
    "   'font.size': 24,\n",
    "   'legend.fontsize': 14,\n",
    "   'xtick.major.width': 3,\n",
    "   'xtick.minor.width': 2,\n",
    "   'xtick.major.size': 12,\n",
    "   'xtick.minor.size': 6,\n",
    "   'xtick.direction': 'in',\n",
    "   'xtick.top': True,\n",
    "   'lines.linewidth':3,\n",
    "   'axes.linewidth':3,\n",
    "   'axes.labelweight':3,\n",
    "   'axes.titleweight':3,\n",
    "   'ytick.major.width':3,\n",
    "   'ytick.minor.width':2,\n",
    "   'ytick.major.size': 12,\n",
    "   'ytick.minor.size': 6,\n",
    "   'ytick.direction': 'in',\n",
    "   'ytick.right': True,\n",
    "   'figure.figsize': [9, 8]\n",
    "   }\n",
    "\n",
    "plt.rcParams.update(params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ace6082",
   "metadata": {},
   "outputs": [],
   "source": [
    "### "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b6c280a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Location of the DC2 Gen3 repository on this site\n",
    "URL = os.getenv('EXTERNAL_INSTANCE_URL')\n",
    "if URL.endswith('data.lsst.cloud'): # IDF\n",
    "    repo = \"s3://butler-us-central1-dp01\"\n",
    "elif URL.endswith('ncsa.illinois.edu'): # NCSA\n",
    "    repo = \"/repo/dc2\"\n",
    "else:\n",
    "    raise Exception(f\"Unrecognized URL: {URL}\")\n",
    "\n",
    "collections=['2.2i/runs/DP0.1']\n",
    "\n",
    "\n",
    "\n",
    "config= os.path.join(repo,'butler.yaml')\n",
    "butler = dafButler.Butler(config=config)\n",
    "registry = butler.registry"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "592ced69",
   "metadata": {},
   "source": [
    "Given this collection we can list the associated reference catalogs.\n",
    "\n",
    "For DP0.1 there is just one: `cal_ref_cat_2_2`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b819fad6",
   "metadata": {},
   "outputs": [],
   "source": [
    "registry.getCollectionSummary('refcats').datasetTypes.names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63dfa976",
   "metadata": {},
   "outputs": [],
   "source": [
    "refDataset='cal_ref_cat_2_2'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a8004a0",
   "metadata": {},
   "source": [
    "For a given dataID we can see what reference datasets are available"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d6df84a",
   "metadata": {},
   "outputs": [],
   "source": [
    "dataId = {'visit': 192350, 'detector': 175, 'band': 'i', 'instrument':'LSSTCam-imSim'}\n",
    "visit=dataId['visit']\n",
    "detector=175"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57d12662",
   "metadata": {},
   "outputs": [],
   "source": [
    "refcatRefs = list(registry.queryDatasets(datasetType=refDataset,\n",
    "                                          collections=[\"refcats\"],\n",
    "                                          #**dataId))\n",
    "                                          instrument=dataId['instrument'],\n",
    "                                          where=f'visit={visit} AND detector={detector}').expanded())\n",
    "refDataIds=[_.dataId for _ in refcatRefs]\n",
    "refCatsDef = [butler.getDeferred(refDataset, __, collections=['refcats']) for __ in refDataIds]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da8b6f88",
   "metadata": {},
   "source": [
    "Then we can load the source catalog data as well as the refcat data, and convert them to astropy tables "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff52ad8b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the source catalog for this visit and convert to astropy table\n",
    "datasetRefs=list(registry.queryDatasets(datasetType='src',\n",
    "                                          collections=\"2.2i/runs/DP0.1\",\n",
    "                                          **dataId))\n",
    "sourceCat = butler.getDirect(datasetRefs[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da4b5c19",
   "metadata": {},
   "outputs": [],
   "source": [
    "#load the associated refcats explicitly \n",
    "refCats=[butler.getDirect(__) for __ in refcatRefs]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dba7f508",
   "metadata": {},
   "outputs": [],
   "source": [
    "#next we plot the two loaded datasets\n",
    "fit,ax=plt.subplots()\n",
    "for refCat in refCats:\n",
    "    ax.scatter(refCat[\"coord_ra\"], refCat[\"coord_dec\"], label=\"refcat\",s=1)\n",
    "plt.scatter(sourceCat[\"coord_ra\"], sourceCat[\"coord_dec\"], label=\"sourcecat\", s=1)\n",
    "plt.legend()\n",
    "plt.xlabel(\"RA\")\n",
    "plt.ylabel(\"DEC\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "89b321d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "refCatsDef"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2315c069",
   "metadata": {},
   "outputs": [],
   "source": [
    "#we can also load the refcat with a spatial query\n",
    "config = LoadReferenceCatalogConfig()\n",
    "config.refObjLoader.ref_dataset_name = refDataset\n",
    "\n",
    "config.refObjLoader.load(os.path.join(getPackageDir('obs_lsst'),\n",
    "                                          'config',\n",
    "                                          'filterMap.py'))\n",
    "config.doApplyColorTerms = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "647a12ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "loaderTask = LoadReferenceCatalogTask(config=config,\n",
    "                                      dataIds=refDataIds,\n",
    "                                      refCats=refCatsDef)\n",
    "\n",
    "# Define center relative to DC2 catalog\n",
    "center = lsst.geom.SpherePoint(np.median(sourceCat['coord_ra']),\n",
    "                               np.median(sourceCat['coord_dec']),\n",
    "                               lsst.geom.radians)\n",
    "# Alternatively, define center relative to reference catalog\n",
    "# center = lsst.geom.SpherePoint(refCats[0]['coord_ra'][0],\n",
    "#                                refCats[0]['coord_dec'][0],\n",
    "#                                lsst.geom.radians)\n",
    "print('Using center (RA, DEC) =', center)\n",
    "\n",
    "refCatSpatial = loaderTask.getSkyCircleCatalog(center,\n",
    "                                         1.0*lsst.geom.degrees,\n",
    "                                         ['i'])\n",
    "print('Found %i reference catalog objects'%(len(refCatSpatial)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "404143c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "fit,ax=plt.subplots()\n",
    "\n",
    "ax.scatter(refCatSpatial[\"ra\"], refCatSpatial[\"dec\"], label=\"refcat\",s=1)\n",
    "plt.scatter(sourceCat[\"coord_ra\"]*180/np.pi, sourceCat[\"coord_dec\"]*180/np.pi, label=\"sourcecat\", s=1)\n",
    "plt.legend()\n",
    "plt.xlabel(\"RA\")\n",
    "plt.ylabel(\"DEC\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e9a3e1a4",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "LSST",
   "language": "python",
   "name": "lsst"
  },
  "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}