diff --git a/Basics/DataInventory.ipynb b/Basics/DataInventory.ipynb
index 394a9a12..a103461d 100644
--- a/Basics/DataInventory.ipynb
+++ b/Basics/DataInventory.ipynb
@@ -8,10 +8,10 @@
}
},
"source": [
- "# Exploring the Shared Datasets in the LSST Science Platform\n",
- "
Owner(s): **Phil Marshall** ([@drphilmarshall](https://github.com/LSSTScienceCollaborations/StackClub/issues/new?body=@drphilmarshall)), \n",
- "
Last Verified to Run: **2018-08-05**\n",
- "
Verified Stack Release: **16.0**\n",
+ "# An Inventory of the Shared Datasets in the LSST Science Platform\n",
+ "
Owner(s): **Phil Marshall** ([@drphilmarshall](https://github.com/LSSTScienceCollaborations/StackClub/issues/new?body=@drphilmarshall)), **Rob Morgan** ([@rmorgan10](https://github.com/LSSTScienceCollaborations/StackClub/issues/new?body=@rmorgan10))\n",
+ "
Last Verified to Run: **2019-08-13**\n",
+ "
Verified Stack Release: **18.1**\n",
"\n",
"In this notebook we'll take a look at some of the datasets available on the LSST Science Platform. \n",
"\n",
@@ -20,12 +20,14 @@
"After working through this tutorial you should be able to: \n",
"1. Start figuring out which of the available datasets is going to be of most use to you in any given project; \n",
"\n",
- "When it is finished, you should be able to:\n",
- "2. Plot the patches and tracts in a given dataset on the sky;\n",
- "3. List the available catalogs in a given dataset.\n",
+ "When it is finished, you should be able to use the `stackclub.Taster` to:\n",
+ "2. Report on the available data in a given dataset;\n",
+ "3. Plot the patches and tracts in a given dataset on the sky.\n",
+ "\n",
+ "**Outstanding Issue:** The `Taster` augments the functionality of the Gen-2 butler, which provides limited capabilities for discovery what data *actually* exist. Specifically, the `Taster` is relying heavily on the `queryMetadata` functionality of the Gen-2 butler, which is limited to a small number of datasets and does not actually guarentee that those datasets exist. The user should beware of over-interpreting the true *existence* of datasets queried by the `Taster`. This should be improved greatly with the Gen-3 butler.\n",
"\n",
"### Logistics\n",
- "This notebook is intended to be runnable on `lsst-lspdev.ncsa.illinois.edu` from a local git clone of https://github.com/LSSTScienceCollaborations/StackClub.\n",
+ "This notebook is intended to be runnable on `lsst-lsp-stable.ncsa.illinois.edu` from a local git clone of https://github.com/LSSTScienceCollaborations/StackClub.\n",
"\n",
"## Set-up"
]
@@ -54,7 +56,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "When editing the `stackclub` package files, we want the latest version to be imported when we re-run the import command. To enable this, we need the %autoreload magic command."
+ "You may need to restart the kernel after doing this. When editing the `stackclub` package files, we want the latest version to be imported when we re-run the import command. To enable this, we need the `%autoreload` magic command."
]
},
{
@@ -75,7 +77,7 @@
}
},
"source": [
- "For accessing the datasets using the Butler, and then visualizing the results, we'll need the following modules:"
+ "To just get a taste of the data that the Butler will deliver for a chosen dataset, we have added a `taster` class to the `stackclub` library. All needed imports are contained in that file, so we only need to import the `stackclub` library to work through this notebook."
]
},
{
@@ -84,37 +86,12 @@
"metadata": {},
"outputs": [],
"source": [
- "%matplotlib inline\n",
- "# %matplotlib ipympl\n",
- "\n",
- "import os, glob\n",
"import numpy as np\n",
- "import matplotlib as mpl\n",
- "import matplotlib.pyplot as plt\n",
- "from IPython.display import IFrame, display, Markdown\n",
+ "%matplotlib inline\n",
"\n",
"import stackclub"
]
},
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": [
- "import lsst.daf.persistence as dafPersist\n",
- "import lsst.daf.base as dafBase\n",
- "\n",
- "import lsst.afw.math as afwMath\n",
- "import lsst.afw.geom as afwGeom\n",
- "\n",
- "import lsst.afw.detection as afwDetect\n",
- "import lsst.afw.image as afwImage\n",
- "import lsst.afw.table as afwTable\n",
- "\n",
- "import lsst.afw.display as afwDisplay"
- ]
- },
{
"cell_type": "markdown",
"metadata": {},
@@ -130,7 +107,7 @@
"source": [
"# What version of the Stack am I using?\n",
"! echo $HOSTNAME\n",
- "! eups list -s | grep lsst_distrib"
+ "! eups list -s lsst_distrib"
]
},
{
@@ -142,14 +119,14 @@
},
"source": [
"## Listing the Available Datasets\n",
- "First, let's look at what is currently available. There are two primary shared dataset folders in the LSP, the read-only `/datasets` folder, and the group-writeable folder `/projects/shared/datasets`. Let's see what's in there"
+ "First, let's look at what is currently available. There are several shared data folders in the LSP, the read-only `/datasets` folder, the project-group-writeable folder `/project/shared/data`, and the Stack Club shared directory `/project/stack-club`. Let's take a look at what's in `/project/shared/data`. Specifically, we want to see butler-friendly data _repositories_, distinguished by their containing a file called `_mapper`, or `repositoryCfg.yaml` in their top level."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "**`/projects/shared/data`:** These datasets are designed to be small test sets, ideal for tutorials."
+ "**`/project/shared/data`:** These datasets are designed to be small test sets, ideal for tutorials."
]
},
{
@@ -158,8 +135,11 @@
"metadata": {},
"outputs": [],
"source": [
- "shared_datasets = ! ls -d /project/shared/data/* | grep -v README\n",
- "shared_datasets"
+ "shared_repos_with_mappers = ! ls -d /project/shared/data/*/_mapper | grep -v README | cut -d'/' -f1-5 | sort | uniq\n",
+ "shared_repos_with_yaml_files = ! ls -d /project/shared/data/*/repositoryCfg.yaml | grep -v README | cut -d'/' -f1-5 | sort | uniq\n",
+ "shared_repos = np.unique(shared_repos_with_mappers + shared_repos_with_yaml_files)\n",
+ "\n",
+ "shared_repos"
]
},
{
@@ -168,11 +148,8 @@
"metadata": {},
"outputs": [],
"source": [
- "%%bash\n",
- "shared_datasets=$( ls -d /project/shared/data/* | grep -v README )\n",
- "for dataset in $shared_datasets; do\n",
- " du -sh $dataset\n",
- "done"
+ "for repo in shared_repos:\n",
+ " ! du -sh $repo"
]
},
{
@@ -189,8 +166,11 @@
"metadata": {},
"outputs": [],
"source": [
- "datasets = ! ls -d /datasets/* | grep -v USAGE | grep -v html\n",
- "datasets"
+ "repos_with_mappers = ! ls -d /datasets/*/repo/_mapper |& grep -v \"No such\" | cut -d'/' -f1-4 | sort | uniq\n",
+ "repos_with_yaml_files = ! ls -d /datasets/*/repo/repositoryCfg.yaml |& grep -v \"No such\" | cut -d'/' -f1-4 | sort | uniq\n",
+ "repos = np.unique(repos_with_mappers + repos_with_yaml_files)\n",
+ "\n",
+ "repos"
]
},
{
@@ -199,19 +179,21 @@
"metadata": {},
"outputs": [],
"source": [
- "# %%bash\n",
- "# datasets=$( ls -d /datasets/* | grep -v USAGE | grep -v html )\n",
- "# for dataset in $datasets; do\n",
- "# du -h $dataset\n",
- "# done"
+ "\"\"\"\n",
+ "for repo in repos:\n",
+ " ! du -sh $repo\n",
+ "\"\"\";"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "## Visualizing Sky Coverage\n",
- "In this section, we'll plot the available patches and tracts in a given dataset on the sky, following the LSST DESC tutorial [dm_butler_skymap.ipynb](https://github.com/LSSTDESC/DC2-analysis/blob/master/tutorials/dm_butler_skymap.ipynb). In fact, we will _import_ this notebook, so that we can re-use its functions. This operation is handled by the `stackclub.wimport` function."
+ "## Exploring the Data Repo with the Stack Club `Taster`\n",
+ "\n",
+ "The `stackclub` library provides a `Taster` class, to explore the datasets in a given repo. As an example, let's take a look at some HSC data using the `Taster`. When instantiating the `Taster`, if you plan to use it for visualizing sky coverage, you can provide it with a path to the tracts from the main repo.\n",
+ "\n",
+ "### Initializing the `Taster`"
]
},
{
@@ -220,16 +202,68 @@
"metadata": {},
"outputs": [],
"source": [
- "dm_butler_skymap_notebook = \"https://github.com/LSSTDESC/DC2-analysis/raw/master/tutorials/dm_butler_skymap.ipynb\"\n",
+ "# Parent repo\n",
+ "repo = '/datasets/hsc/repo/'\n",
"\n",
- "skymapper = stackclub.wimport(dm_butler_skymap_notebook, vb=True)"
+ "#Location of tracts for a particular rerun and depth relative to main repo\n",
+ "rerun = 'DM-13666' # DM-13666, DM-10404 \n",
+ "depth = 'WIDE' # WIDE, DEEP, UDEEP\n",
+ "tract_location = 'rerun/' + rerun + '/' + depth"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Execute one of the following two cells. The latter will make `tarquin` aware of the tracts for the dataset while the former will just look at the repo as a whole and not visualize any sky area."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "tarquin = stackclub.Taster(repo, vb=True)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "tarquin = stackclub.Taster(repo, vb=True, path_to_tracts=tract_location)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Properties of the `Taster`"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The taster, `tarquin`, carries a butler around with it:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "type(tarquin.butler)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "> BUG: remote notebooks are not yet `wimport`-able. A workaround could be to import the downloaded file explicitly. This is not yet working, hence the commented out failed attempt below."
+ "If we ask the taster to investigate a folder that is not a repo, its butler will be `None`"
]
},
{
@@ -238,9 +272,7 @@
"metadata": {},
"outputs": [],
"source": [
- "# import sys, os\n",
- "# import stackclub\n",
- "# sys.path.append(os.getcwd() + '/.downloads')"
+ "failed = stackclub.Taster('not-a-repo', vb=True)"
]
},
{
@@ -249,14 +281,14 @@
"metadata": {},
"outputs": [],
"source": [
- "# import dm_butler_skymap"
+ "print(failed.butler)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
- "Now we can attempt to plot the available tracts, using the `plot_skymap_tract()` function."
+ "The taster uses its butler to query the metadata of the repo for datasets, skymaps etc. "
]
},
{
@@ -265,17 +297,119 @@
"metadata": {},
"outputs": [],
"source": [
- "# repo = \"/project/shared/data/Twinkles_subset/output_data_v2\"\n",
- "repo = \"/datasets/hsc/repo/rerun/DM-13666/WIDE\"\n",
- "butler = dafPersist.Butler(repo)\n",
+ "tarquin.look_for_datasets_of_type(['raw', 'calexp', 'deepCoadd_calexp', 'deepCoadd_mergeDet'])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "> **PROBLEM: these last two datatypes are not listed in the repo metadata. This is one of the issues with the Gen-2 butler and the`Taster` is not smart enough to search the tract folders for catalog files. This should be updated with Gen-3.**"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "tarquin.look_for_skymap()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The `what_exists` method searches for everything \"interesting\". In the `taster.py` class, interesting currently consists of \n",
+ "* `'raw'`\n",
+ "* `'calexp'` \n",
+ "* `'src'`\n",
+ "* `'deepCoadd_calexp'`\n",
+ "* `'deepCoadd_meas'` \n",
"\n",
- "# Glob the merged coadd folder for the tracts that have data. Unfortunately, this information is not\n",
- "# directly accessible from the data butler.\n",
- "tracts = sorted([int(os.path.basename(x)) for x in\n",
- " glob.glob(os.path.join(repo, 'deepCoadd-results', 'merged', '*'))])\n",
+ "but this method can easily be updated to include more dataset types."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "tarquin.what_exists()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "If one wishes to check the existance of all dataset types, you can use the `all` parameter of the `what_exists()` method to do exactly that. Checking all dataset types may take a minute or so (while the `Taster` does a lot of database queries)."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "tarquin.what_exists(all=True)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "A dictionary with existence information is stored in the `exists` attribute:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "tarquin.exists"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The `Taster` can report on the data available, counting the number of visits, sources, etc, according to what's in the repo. It uses methods like this one:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "tarquin.estimate_sky_area()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "and this one:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "tarquin.count_things()\n",
"\n",
- "# How many tracts do we have?\n",
- "print(\"Found {} tracts\".format(len(tracts)))"
+ "print(tarquin.counts)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "When the `estimate_sky_area` method runs, `tarquin` collects all the tracts associated with the repo. A list of the tracts is stored in the attribute `tarquin.tracts`."
]
},
{
@@ -284,15 +418,108 @@
"metadata": {},
"outputs": [],
"source": [
- "\"\"\"\n",
- "Uncomment this cell when the `wimport` bug is fixed (or avoided).\n",
+ "tarquin.tracts"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Using the tracts, we can get a rough estimate for what parts of the sky have been targeted in the dataset. The method for doing this is `tarquin.plot_sky_coverage`, and follows the example code given in [Exploring_A_Data_Repo.ipynb](Exploring_A_Data_Repo.ipynb)."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "tarquin.plot_sky_coverage()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "To have your `Taster` do all the above, and just report on what it finds, do:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "tarquin.report()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "If you are interested in learning which fields, filters, visits, etc. have been counted by `tarquin`, remember that `tarquin` carries an instance of the `Butler` with it, so you can run typical `Butler` methods. For example, if you found the number of filters being 13 odd, you can look at the filters like this:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "tarquin.butler.queryMetadata('calexp', ['filter'])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "For more on the `Taster`'s methods, do, for example:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# help(tarquin)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Example Tastings"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Let's compare the WIDE, DEEP and UDEEP parts of the HSC dataset."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "repo = '/datasets/hsc/repo/'\n",
+ "rerun = 'DM-13666'\n",
"\n",
- "# Now, loop over all the tracts, plotting them as gray, numbered, rectangles:\n",
- "ax = None\n",
- "for tract in tracts:\n",
- " skyMap = butler.get('deepCoadd_skyMap')\n",
- " ax = skymapper.plot_skymap_tract(skyMap, tract=tract, title='', ax=ax)\n",
- "\"\"\";"
+ "for depth in ['WIDE', 'DEEP', 'UDEEP']:\n",
+ " tract_location = 'rerun/' + rerun + '/' + depth\n",
+ " \n",
+ " taster = stackclub.Taster(repo, path_to_tracts=tract_location)\n",
+ " taster.report()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "You may notice that all **Metadata Characteristics** beginning with \"Number of\" are the same for the three depths. This is a result of `tarquin`'s `Butler` getting this information from the repo as a whole, rather than the specific depth we specified for the tracts. There is more information on why the `Butler` works in this way in the [Exploring_A_Data_Repo.ipynb](https://github.com/LSSTScienceCollaborations/StackClub/blob/project/data_inventory/drphilmarshall/Basics/Exploring_A_Data_Repo.ipynb) notebook."
]
},
{
@@ -301,7 +528,57 @@
"source": [
"## Summary\n",
"\n",
- "In this notebook we took a first look at the datasets available to us in two shared directories in the LSST science platform filesystem."
+ "In this notebook we took a first look at the datasets available to us in two shared directories in the LSST science platform filesystem, and used the `stackclub.Taster` class to report on their basic properties, and their sky coverage. Details on the methods used by the `Taster` can be found in the [Exploring_A_Data_Repo.ipynb](https://github.com/LSSTScienceCollaborations/StackClub/blob/project/data_inventory/drphilmarshall/Basics/Exploring_A_Data_Repo.ipynb) notebook, or by executing the follwoing cell:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "help(tarquin)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# STILL TODO\n",
+ "* Build defensiveness into the `Taster` so that it can handle a wider variety of datasets.\n",
+ "* Update `Taster` to use Gen-3 butler\n",
+ "\n",
+ "### Looking at other shared datasets and repos\n",
+ "\n",
+ "The following loops over all shared datasets fails in interesting ways: some folders don't seem to be `Butler`-friendly. We need to do a bit more work to identify the actual repos available to us, and then use the `Taster` to provide a guide to all of them."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "for repo in shared_repos:\n",
+ " try:\n",
+ " taster = stackclub.Taster(repo)\n",
+ " taster.report()\n",
+ " except:\n",
+ " print(\"Taster failed to explore repo \",repo)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "for repo in repos:\n",
+ " try:\n",
+ " taster = stackclub.Taster(repo)\n",
+ " taster.report()\n",
+ " except:\n",
+ " print(\"Taster failed to explore repo \",repo)"
]
}
],
@@ -322,7 +599,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
- "version": "3.6.2"
+ "version": "3.7.2"
},
"livereveal": {
"scroll": true,
diff --git a/Basics/Exploring_A_Data_Repo.ipynb b/Basics/Exploring_A_Data_Repo.ipynb
new file mode 100644
index 00000000..df203aab
--- /dev/null
+++ b/Basics/Exploring_A_Data_Repo.ipynb
@@ -0,0 +1,827 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Exploring a Data Repository\n",
+ "\n",
+ "
Owner: **Rob Morgan** ([@rmorgan10](https://github.com/LSSTScienceCollaborations/StackClub/issues/new?body=@rmorgan10)), **Phil Marshall** ([@drphilmarshall](https://github.com/LSSTScienceCollaborations/StackClub/issues/new?body=@drphilmarshall)), **Alex Drlica-Wagner** ([@kadrlica](https://github.com/LSSTScienceCollaborations/StackClub/issues/new?body=@kadrlica))\n",
+ "
Last Verified to Run: **2019-08-13**\n",
+ "
Verified Stack Release: **18.1**\n",
+ "\n",
+ "This notebook examines the content of a data repository and shows how to determine the inputs for each component. \n",
+ "\n",
+ "### Learning Objectives:\n",
+ "After working through and studying this notebook you should be able to understand how to use the Butler to figure out: \n",
+ " 1. What a data repo is;\n",
+ " 2. Which data types are present in a data repository;\n",
+ " 3. If coadds have been made, what the available tracts are;\n",
+ " 4. Which parts of the sky those tracts cover.\n",
+ " \n",
+ "### Logistics\n",
+ "This notebook is intended to be runnable on `lsst-lsp-stable.ncsa.illinois.edu` from a local git clone of https://github.com/LSSTScienceCollaborations/StackClub.\n",
+ "\n",
+ "\n",
+ "## Set Up"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import os\n",
+ "import sys\n",
+ "import warnings\n",
+ "import matplotlib.pyplot as plt\n",
+ "from IPython.display import display, Markdown\n",
+ "import numpy as np\n",
+ "import os, glob\n",
+ "%matplotlib inline\n",
+ "\n",
+ "# Filter some warnings printed by v16.0 of the stack\n",
+ "# warnings.simplefilter(\"ignore\", category=FutureWarning)\n",
+ "# warnings.simplefilter(\"ignore\", category=UserWarning)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## What is a Data Repo?\n",
+ "\n",
+ "A data repo is a directory containing raw images, calibration files, metadata and configuration information, defining an LSST-format dataset. Data repositories contain either a `_mapper` file or a `repositoryCfg.yaml` file, which record the \"obs package\" that was used to organize the data. The obs package gives the repository more structure and organization than an ordinary data directory. Let's take a look at this file structure in the HSC data repo."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### The HSC Data Repo: What's in there?"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "repo = '/datasets/hsc/repo'"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We'll use the `hsc` data repository as our testing ground, and start by figuring out what it contains. In the `hsc` case, the `_mapper` file is in the top level folder, while the data repo for each field is a few levels down."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "! ls /datasets/hsc/repo/"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We can see the `_mapper` file here, and at contains one line giving the name of the `Mapper` object for the HSC repo:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "! cat /datasets/hsc/repo/_mapper\n",
+ "\n",
+ "# Import the Mapper object once you know its name\n",
+ "from lsst.obs.hsc import HscMapper"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "You can get some more information on this object like this:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# help(HscMapper)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The mapper defines a (large) number of different \"dataset types\". Some of these are specific to this particular data repo, others are more general. Even filtering out some intermediate dataset types, we are still left with a long list. But, once we figure out which dataset types we are interested in, we can start querying for information about those datasets."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "mapper = HscMapper(root=repo)\n",
+ "all_dataset_types = mapper.getDatasetTypes()\n",
+ "\n",
+ "remove = ['_config', '_filename', '_md', '_sub', '_len', '_schema', '_metadata']\n",
+ "\n",
+ "shortlist = []\n",
+ "for dataset_type in all_dataset_types:\n",
+ " keep = True\n",
+ " for word in remove:\n",
+ " if word in dataset_type:\n",
+ " keep = False\n",
+ " if keep:\n",
+ " shortlist.append(dataset_type)\n",
+ "\n",
+ "print(shortlist)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The `Butler`, directed by the `Mapper`, will have access to all the above dataset types. \n",
+ "\n",
+ "Another important file in the repo parent folder is `registry.sqlite3`. This database contains metadata for the HSC **raw** images."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Great, but where is the actual data, and how was it processed?"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The raw visit images are stored by field. In the HSC dataset the fields have names like `COSMOS` and `DEEPE09`. Within those field folders, there is a directory structure that eventually gets down to visit image FITS files whose names and paths contain the date/time and filter for that exposure. For example:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "! \\ls /datasets/hsc/repo/COSMOS/2015-01-18/01113/HSC-Y/HSC-0018476-00?.fits"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "How do the pipeline tasks know which raw data to process? This information is captured in the \"configs\". In the HSC repo there is no config folder or files in the top level directory - in fact the only two files are `_mapper` and `registry.sqlite3`. So what's going on? \n",
+ "\n",
+ "It turns out that the provenance of the stack processing of the HSC raw images is captured in \"rerun\" folders, one for each time the science pipelines were run on the data. Let's do some detective work to find out what happened to the HSC data. "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "! ls /datasets/hsc/repo/rerun"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "First, `DM-10404` looks like a run ID. What's in that folder? "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "! ls /datasets/hsc/repo/rerun/DM-10404"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "`DEEP`, `WIDE`, and `UDEEP` are the names of the sub-surveys of the HSC survey. We might expect each to contain results from the processing of that sub-survey's images."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "! ls /datasets/hsc/repo/rerun/DM-10404/UDEEP"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The numerically-named folders contain the generated catalog files, organized by sky tract. Note that a `config` folder is present, and also a `repositoryCfg.yaml` file - which means that this folder is itself a `repo`, from the Butler's point of view."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "! cat /datasets/hsc/repo/rerun/DM-10404/WIDE/repositoryCfg.yaml"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "What does it mean that `SFM` is the \"parent\" of this repo? Let's see what _that_ folder contains:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "! ls /datasets/hsc/repo/rerun/DM-10404/SFM"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "`SFM` seems to contain _all_ the tracts that have been produced - so is the entire HSC survey. It's `repositoryCfg.yaml` file shows that it's \"parent\" is the top level folder, `/datasets/hsc/repo/`."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "! cat /datasets/hsc/repo/rerun/DM-10404/SFM/repositoryCfg.yaml"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Let's see what the `DM-10404/UDEEP` repo's `config` folder contains:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "! ls -a /datasets/hsc/repo/rerun/DM-10404/UDEEP/config"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "These are the configuration files that were used when the science pipelines were run on these data. While we don't know which versions of the software were used, we at least know which tasks were run. \n",
+ "\n",
+ "Here's what a config file looks like (ignoring the many import statements and just looking at a few example lines):"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "! cat /datasets/hsc/repo/rerun/DM-10404/UDEEP/config/forcedPhotCcd.py | grep -v import | head -15"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The next level to dig into here is the tract folders within one of the repos in this rerun. For example:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "! ls /datasets/hsc/repo/rerun/DM-10404/UDEEP/00814/HSC-Y/tract9570 | head -10"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Those FITS files contain the forced source tables."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Instantiating the Butler and looking for Dataset Types\n",
+ "\n",
+ "Now that we have an idea of the structure of the repo itself, let's use the Butler to explore the data within the repo. Here we will demonstrate a few useful `Butler` methods for learning about the data in a repo. Let's choose one of the rerun repos, and investigate its properties. We'll summon two butlers, one that is pointed at the parent repo, and another (an \"under butler\") that is asked to focus on a particular sub-survey in a particular re-run."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "parent_repo = '/datasets/hsc/repo'\n",
+ "\n",
+ "# Choose a re-run repo:\n",
+ "rerun_id = 'DM-10404'\n",
+ "depth = 'UDEEP'\n",
+ "# Try a different one:\n",
+ "# rerun = 'DM-13666'\n",
+ "# depth = 'WIDE'\n",
+ "\n",
+ "repo = parent_repo + '/rerun/' + rerun_id + '/' + depth\n",
+ "print(repo)\n",
+ "\n",
+ "from lsst.daf.persistence import Butler\n",
+ "\n",
+ "butler = Butler(parent_repo)\n",
+ "under_butler = Butler(repo)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The `butler` can check whether a datatype (like the source catalogs) actually exists or not, but it needs a specific dataset ID to check whether that specific part of the dataset exists. \n",
+ "\n",
+ "Note that the metadata being queried here is in the `registry.sqlite3` database in the _parent_ repo - and so refers to the _initial_ processing run, not the most recent rerun. We'll need to work carefully around this below."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Obtaining Basic Dataset Properties Using the Butler\n",
+ "Now we can start using Butler methods to query metadata for the repo. For this dataset, we can look at the filters used, number of visits, number of pointings, etc. by examining the Butler's keys and metadata. For these basic properties, we will look at the `calexp` and `src` tables. The contents of these tables are derived from the processing of individual sensors, and exist in the parent folder. (That means that we can use either of our two butlers to query for them.)\n",
+ "\n",
+ "Note that the metadata is created from the raw exposures loaded into the sqlite registry. The fact that we can get metadata for a specific datasetType and dataId **does not** imply that the data exist on disk (we will check this in a subsequent step)."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# This would be faster if only one query were issued...\n",
+ "visits = butler.queryMetadata('calexp', ['visit'])\n",
+ "pointings = butler.queryMetadata('calexp', ['pointing'])\n",
+ "ccds = butler.queryMetadata('calexp', ['ccd'])\n",
+ "fields = butler.queryMetadata('calexp', ['field'])\n",
+ "filters = butler.queryMetadata('calexp', ['filter'])\n",
+ "sources = butler.queryMetadata('src', ['id'])\n",
+ "\n",
+ "# It is possible to specify multiple formats -- i.e., butler.queryMetadata('calexp', ['visit','ccd'])\n",
+ "metadata = butler.queryMetadata('calexp', ['visit','pointing','ccd','field','filter'])"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "num_visits = len(visits)\n",
+ "num_pointings = len(pointings)\n",
+ "num_ccds = len(ccds)\n",
+ "num_fields = len(fields)\n",
+ "num_filters = len(filters)\n",
+ "num_sources = len(sources)\n",
+ "num_metadata = len(metadata)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "print(\"The HSC {}/{} rerun contains {} visits.\".format(rerun_id,depth,num_visits))\n",
+ "print(\"The HSC {}/{} rerun contains {} pointings.\".format(rerun_id,depth,num_pointings))\n",
+ "print(\"The HSC {}/{} rerun contains {} ccds.\".format(rerun_id,depth,num_ccds))\n",
+ "print(\"The HSC {}/{} rerun contains {} fields.\".format(rerun_id,depth,num_fields))\n",
+ "print(\"The HSC {}/{} rerun contains {} filters.\".format(rerun_id,depth,num_filters))\n",
+ "print(\"The HSC {}/{} rerun contains {} sources.\".format(rerun_id,depth,num_sources))\n",
+ "# print(\"The HSC {}/{} rerun contains {} coadd sources.\".format(rerun_id,depth,num_coadd_sources))\n",
+ "# print(\"The HSC {}/{} rerun contains {} forced sources.\".format(rerun_id,depth,num_forced_sources))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "As a check, let's compare what our two butlers find when asked for the number of sources:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "alt_num_sources = len(under_butler.queryMetadata('src', ['id']))\n",
+ "\n",
+ "print(\"The butler says that we have {:d} input sources.\".format(num_sources))\n",
+ "print(\"The under butler says we have {:d} input sources.\".format(alt_num_sources))\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "So, we could have done our query with the `under_butler` as well. In practice, it's best to specify a Butler for the rerun repo, because that Butler will also have access to the parent repo.\n",
+ "\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### So where are the Calexps and Source Catalogs?\n",
+ "\n",
+ "Notice that while we were able to get metadata for the `processCcd` outputs (the `calexp` and `src`), that **does not** guarantee that these products are on disk. The metadata is created from the raw inputs stored in the registry and a template for the derived data products. In this case, we have a reason to be suspicious: remember that when we examined the directory structure we did not see a directory for the `processCcd` products. \n",
+ "\n",
+ "To check the existence of the data requires the use of the `datasetExists` method of the `butler`. Let's give this a try."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "First we select the metadata for a specific calexp"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "dataId={'filter':'HSC-I','visit':872,'ccd':10}\n",
+ "butler.queryMetadata('calexp', ['visit','ccd','filter','field','pointing'], dataId=dataId)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Ok, so we have the metadata, let's try to get the calexp..."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "try: \n",
+ " butler.get('calexp',dataId=dataId)\n",
+ "except Exception as e:\n",
+ " print(str(e))\n",
+ " print(\"\\nWhat? Does the data exist?\")\n",
+ "\n",
+ "# Explicitly check for existence\n",
+ "exists = butler.datasetExists('calexp',dataId=dataId)\n",
+ "print(\"\\nbutler.datasetExists: \" + str(exists))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "So it appears that we don't have these intermediate products stored in this repo."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Coadd Sky Area\n",
+ "\n",
+ "One may also be interested in the total sky area imaged for a particular coadd rerun/depth. We can estimate and visualize this from the coadd tract info that neither our `under_butler` nor our `butler` has access to. To collect all the tracts, we have to get them via the file structure. This operation will hopefully be `Butler`-ized with the Gen3 Butler."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Collect tract indices from file names\n",
+ "tracts = sorted([int(os.path.basename(x)) for x in\n",
+ " glob.glob(os.path.join(repo, 'deepCoadd-results', 'merged', '*'))])\n",
+ "num_tracts = len(tracts)\n",
+ "\n",
+ "print(\"Found {} merged tracts in repo {}\".format(num_tracts, repo))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "A quick way of extimating the sky area covered is to sum the areas of the inner boxes of all the tracts. For more information on the properties of tracts, you can look at the [Documentation](http://doxygen.lsst.codes/stack/doxygen/x_masterDoxyDoc/classlsst_1_1skymap_1_1tract_info_1_1_tract_info.html)."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "As a quick note, the file structure only tells us the names of the tracts in the particular rerun/depth to look at. The actual `TractInfo` objects are obtained by selecting the tracts we want from the `deepCoadd_skyMap` dataset in our particular rerun repo. Therefore, we will have to ask the `under_butler` to bring us this dataset for the particular rerun/depth. "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Calculate area from all tracts\n",
+ "skyMap = under_butler.get('deepCoadd_skyMap')\n",
+ "total_area = 0.0 #deg^2\n",
+ "plotting_vertices = []\n",
+ "for test_tract in tracts:\n",
+ " # Get inner vertices for tract\n",
+ " tractInfo = skyMap[test_tract]\n",
+ " vertices = tractInfo._vertexCoordList\n",
+ " plotting_vertices.append(vertices)\n",
+ " \n",
+ " #calculate area of box\n",
+ " av_dec = 0.5 * (vertices[2][1] + vertices[0][1])\n",
+ " av_dec = av_dec.asRadians()\n",
+ " delta_ra_raw = vertices[0][0] - vertices[1][0] \n",
+ " delta_ra = delta_ra_raw.asDegrees() * np.cos(av_dec)\n",
+ " delta_dec= vertices[2][1] - vertices[0][1]\n",
+ " area = delta_ra * delta_dec.asDegrees()\n",
+ " \n",
+ " #combine areas\n",
+ " total_area += area\n",
+ " \n",
+ "\n",
+ "# Round off the total area for presentation purposes\n",
+ "rounded_total_area = round(total_area, 2)\n",
+ "\n",
+ "print(\"Total area imaged (sq deg): \",rounded_total_area)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Displaying Dataset Characteristics"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Now let's print out a report of all the characteristcs we have found. We'll use the sky area from the rerun we chose, and the numbers common to all reruns."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "dataset_name = 'HSC'\n",
+ "display(Markdown('### %s' % repo))\n",
+ "\n",
+ "# Make a table of the collected metadata\n",
+ "collected_data = [num_visits, num_pointings, num_ccds, num_fields, num_filters, num_sources, \n",
+ " num_tracts, rounded_total_area]\n",
+ "data_names = (\"Number of Visits\", \"Number of Pointings\", \"Number of CCDs\", \"Number of Fields\", \n",
+ " \"Number of Filters\", \"Number of Sources\", \"Number of Tracts\", \"Total Sky Area (deg$^2$)\")\n",
+ "# TODO: include coadd sources and forced sources\n",
+ "\n",
+ "output_table = \"| Metadata Characteristic | Value | \\n | ---: | ---: | \\n \"\n",
+ "counter = 0\n",
+ "while counter < len(collected_data):\n",
+ " output_table += \"| %s | %s | \\n\" %(data_names[counter], collected_data[counter])\n",
+ " counter += 1\n",
+ "display(Markdown(output_table))\n",
+ "\n",
+ "# Show which fields and filters we're talking about:\n",
+ "display(Markdown('Fields: (%i total)' %num_fields))\n",
+ "print(fields)\n",
+ "display(Markdown('Filters: (%i total)' %num_filters))\n",
+ "print(filters)\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Plotting the sky coverage\n",
+ "\n",
+ "For this we will need our list of merged `tracts` from above, and also the `skyMap` object. We can then extract the sky coordinates of the corners of each tract, and use them to draw a set of rectangles to illustrate the sky coverage, following Jim Chiang's LSST DESC tutorial [dm_butler_skymap.ipynb](https://github.com/LSSTDESC/DC2-analysis/blob/master/tutorials/dm_butler_skymap.ipynb).\n",
+ "\n",
+ "In the future, we could imagine overlaying the focal plane and color the individual visits, using more of the code from Jim's notebook. Let's wait to see what functionality the Gen3 Butler provides first, and maybe return to visualization later."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.figure()\n",
+ "\n",
+ "for tract in tracts:\n",
+ " tractInfo = skyMap[tract]\n",
+ " \n",
+ " corners = [(x[0].asDegrees(), x[1].asDegrees()) for x in tractInfo.getVertexList()]\n",
+ " x = [k[0] for k in corners] + [corners[0][0]]\n",
+ " y = [k[1] for k in corners] + [corners[0][1]]\n",
+ " \n",
+ " \n",
+ " plt.plot(x,y, color='b')\n",
+ " \n",
+ "plt.xlabel('RA (deg)')\n",
+ "plt.ylabel('Dec (deg)')\n",
+ "plt.title('2D Projection of Sky Coverage')\n",
+ "\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We could imagine plotting the patches as well, to show which tracts were incomplete - but this gives us a rough idea of where our data is on the sky."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Finding a Coadd Object Catalog\n",
+ "\n",
+ "Since we know we have `deepCoadd` images, we might want to find the sources detected in those images. These are the precursors of `Objects`. If forced photometry has been run on these sources, there should be `deepCoadd_forced_src` catalogs present. The `config` tells us which tasks have been run in the rerun, but metadata describing the results of the run is not propagated back into the registry database in the parent folder. So, we need to get data from the `under_butler`.\n",
+ "\n",
+ "Since the registry only knows about the raw exposures (and products that can be directly derived from templates using the raws), we can't get valid metadata from `queryMetadata` directly. Instead, we query the directory structure for a valid filter, tract, patch combination."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "dirname=os.path.join(repo,'deepCoadd-results')\n",
+ "print('filters: '+dirname)\n",
+ "!ls $dirname\n",
+ "filter='HSC-I'\n",
+ "\n",
+ "dirname = os.path.join(dirname, filter)\n",
+ "print('\\ntracts: '+dirname)\n",
+ "!ls $dirname\n",
+ "tract=8523\n",
+ "\n",
+ "dirname = os.path.join(dirname,str(tract))\n",
+ "print('\\npatches: '+dirname)\n",
+ "!ls $dirname\n",
+ "patch='1,3'\n",
+ "\n",
+ "dirname = os.path.join(dirname,patch)\n",
+ "print('\\ndata products: '+dirname)\n",
+ "!ls $dirname"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We can also get some available tract, patches from the `skyMap` object. We need to feed it a list of coordinates defining a polygon to search for viable tract,patch info."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from lsst.geom import SpherePoint, Angle\n",
+ "# Create a list of corners of the polygon to search.\n",
+ "coordList = [SpherePoint(Angle(np.radians(0)),Angle(np.radians(-6))),\n",
+ " SpherePoint(Angle(np.radians(0)),Angle(np.radians(-2))),\n",
+ " SpherePoint(Angle(np.radians(40)),Angle(np.radians(-2))),\n",
+ " SpherePoint(Angle(np.radians(40)),Angle(np.radians(-6)))\n",
+ " ]\n",
+ "# Print the tract,patch info\n",
+ "tractInfo = skyMap.findTractPatchList(coordList)\n",
+ "for _tract in tractInfo:\n",
+ " print(_tract[0])\n",
+ " for _patch in _tract[1]:\n",
+ " print(' ',_patch)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Here we grab the coadd source catalog for a specific filter\n",
+ "dataId={'filter':'HSC-I','tract':tract,'patch':patch}\n",
+ "coadd_sources = under_butler.get('deepCoadd_forced_src',dataId=dataId)\n",
+ "coadd_sources.asAstropy()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# We can also get the merged detections\n",
+ "dataId={'tract':8523,'patch':'1,8'}\n",
+ "merged_sources = under_butler.get('deepCoadd_mergeDet',dataId=dataId)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Visualizing Simple Measurements\n",
+ "\n",
+ "We'd like to know something about the sources that have been detected and measured, beyond a simple total number. We create a simple histogram of the Kron flux of sources in our HSC-I band tile. You can find more details on merging catalogs across bands, etc on the [Science Pipelines documentation](https://pipelines.lsst.io/getting-started/multiband-analysis.html)."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "plt.figure()\n",
+ "plt.hist(coadd_sources['ext_photometryKron_KronFlux_flux'],bins=np.linspace(0,100,50))\n",
+ "plt.xlabel(\"log10(Kron Flux)\")\n",
+ "plt.title(\"{filter} {tract} {patch}\".format(filter=filter,tract=tract,patch=patch))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Summary\n",
+ "\n",
+ "We have shown a few techniques for exploring a data repo. To make this process straightforward, we are implementing all these techniques into methods of a `Taster` class, which is now a part of the `stackclub` library. The `Taster` will give you a taste of what the `Butler` delivers. we demonstrate the use of this class in the [DataInventory.ipynb](https://github.com/LSSTScienceCollaborations/StackClub/blob/project/data_inventory/drphilmarshall/Basics/DataInventory.ipynb) notebook.\n"
+ ]
+ }
+ ],
+ "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.7.2"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/Basics/dm_butler_skymap.ipynb b/Basics/dm_butler_skymap.ipynb
new file mode 100644
index 00000000..fb573e9e
--- /dev/null
+++ b/Basics/dm_butler_skymap.ipynb
@@ -0,0 +1,415 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Plotting the DC2 Run1.1p skyMap\n",
+ "
Owner: **Jim Chiang** ([@jchiang87](https://github.com/LSSTDESC/DC2-analysis/issues/new?body=@jchiang87))\n",
+ "
Last Verified to Run: **2018-10-26** (by @yymao)\n",
+ "\n",
+ "In this notebook, we show how to use the data butler to obtain information on the skyMap used in the coadd analyses performed by the DRP pipeline. These skyMaps are composed of tracts and patches on the sky. Each tract is a rectangular region of the sky with a common map projection; each tract is further divided into rectangular patches, which use the same tract coordinate system and which are a convenient size for processing the coadd data. A more complete description of the skyMap geometry is given in the HSC Software Pipeline paper ([Bosch et al. 2017](https://arxiv.org/abs/1705.06766)).\n",
+ "\n",
+ "Equipped with the info from the skyMap, we plot the tracts and patches that were used with the Run1.1p processing and overlay the WFD and uDDF simulation regions. We also use the butler to access the visit-level data and show how one can access the calexp (calibrated exposure) image data to obtain the PSF, zero-point, etc. as measured by the Stack for a given exposure. Finally, we show how to plot the sky region imaged on the focal plane for a given visit in two ways: the first using the CCD coordinates available from the calexps and the other using the lsst_sims code to compute those coordinates from the pointing information for the visit.\n",
+ "\n",
+ "## Set Up"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import os\n",
+ "import glob\n",
+ "import warnings\n",
+ "import sqlite3\n",
+ "import re\n",
+ "\n",
+ "import numpy as np\n",
+ "import matplotlib.pyplot as plt\n",
+ "from matplotlib.path import Path\n",
+ "import matplotlib.patches as patches\n",
+ "%matplotlib inline\n",
+ "\n",
+ "import lsst.afw.geom as afw_geom\n",
+ "import lsst.afw.cameraGeom as cameraGeom\n",
+ "import lsst.daf.persistence as dp\n",
+ "# The lsst_sims code issues some ignorable warnings regarding ids used for querying the object\n",
+ "# databases.\n",
+ "with warnings.catch_warnings():\n",
+ " warnings.simplefilter(\"ignore\")\n",
+ " import lsst.sims.coordUtils\n",
+ " from lsst.sims.catUtils.utils import ObservationMetaDataGenerator\n",
+ " from lsst.sims.utils import getRotSkyPos\n",
+ " \n",
+ "from desc_dc2_dm_data import REPOS"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Plotting Functions\n",
+ "First, we define some functions to plot the tract, patch, and CCD regions on the sky. These are copied from [example code](https://github.com/yalsayyad/dm_notebooks/blob/master/desc-ssim/DESC-SSim%20Patch%20Geometry.ipynb) that Yusra AlSayyad presented at the [2017-06-29 SSim meeting](https://confluence.slac.stanford.edu/pages/viewpage.action?pageId=224461017)."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def make_patch(vertexList, wcs=None):\n",
+ " \"\"\"\n",
+ " Return a Path in sky coords from vertex list in pixel coords.\n",
+ " \n",
+ " Parameters\n",
+ " ----------\n",
+ " vertexList: list of coordinates\n",
+ " These are the corners of the region to be plotted either in pixel coordinates or\n",
+ " sky coordinates.\n",
+ " wcs: lsst.afw.geom.skyWcs.skyWcs.SkyWcs [None]\n",
+ " The WCS object used to convert from pixel to sky coordinates.\n",
+ "\n",
+ " Returns\n",
+ " -------\n",
+ " matplotlib.path.Path: The encapsulation of the vertex info that matplotlib uses to\n",
+ " plot a patch.\n",
+ " \"\"\"\n",
+ " if wcs is not None:\n",
+ " skyPatchList = [wcs.pixelToSky(pos).getPosition(afw_geom.degrees)\n",
+ " for pos in vertexList]\n",
+ " else:\n",
+ " skyPatchList = vertexList\n",
+ " verts = [(coord[0], coord[1]) for coord in skyPatchList]\n",
+ " verts.append((0,0))\n",
+ " codes = [Path.MOVETO,\n",
+ " Path.LINETO,\n",
+ " Path.LINETO,\n",
+ " Path.LINETO,\n",
+ " Path.CLOSEPOLY,\n",
+ " ]\n",
+ " return Path(verts, codes)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def plot_skymap_tract(skyMap, tract=0, title=None, ax=None):\n",
+ " \"\"\"\n",
+ " Plot a tract from a skyMap.\n",
+ " \n",
+ " Parameters\n",
+ " ----------\n",
+ " skyMap: lsst.skyMap.SkyMap\n",
+ " The SkyMap object containing the tract and patch information.\n",
+ " tract: int [0]\n",
+ " The tract id of the desired tract to plot.\n",
+ " title: str [None]\n",
+ " Title of the tract plot. If None, the use `tract `.\n",
+ " ax: matplotlib.axes._subplots.AxesSubplot [None]\n",
+ " The subplot object to contain the tract plot. If None, then make a new one.\n",
+ "\n",
+ " Returns\n",
+ " -------\n",
+ " matplotlib.axes._subplots.AxesSubplot: The subplot containing the tract plot.\n",
+ " \"\"\"\n",
+ " if title is None:\n",
+ " title = 'tract {}'.format(tract)\n",
+ " tractInfo = skyMap[tract]\n",
+ " tractBox = afw_geom.Box2D(tractInfo.getBBox())\n",
+ " tractPosList = tractBox.getCorners()\n",
+ " wcs = tractInfo.getWcs()\n",
+ " xNum, yNum = tractInfo.getNumPatches()\n",
+ "\n",
+ " if ax is None:\n",
+ " fig = plt.figure(figsize=(12,8))\n",
+ " ax = fig.add_subplot(111)\n",
+ "\n",
+ " tract_center = wcs.pixelToSky(tractBox.getCenter())\\\n",
+ " .getPosition(afw_geom.degrees)\n",
+ " ax.text(tract_center[0], tract_center[1], '%d' % tract, size=16,\n",
+ " ha=\"center\", va=\"center\", color='blue')\n",
+ " for x in range(xNum):\n",
+ " for y in range(yNum):\n",
+ " patchInfo = tractInfo.getPatchInfo([x, y])\n",
+ " patchBox = afw_geom.Box2D(patchInfo.getOuterBBox())\n",
+ " pixelPatchList = patchBox.getCorners()\n",
+ " path = make_patch(pixelPatchList, wcs)\n",
+ " patch = patches.PathPatch(path, alpha=0.1, lw=1)\n",
+ " ax.add_patch(patch)\n",
+ " center = wcs.pixelToSky(patchBox.getCenter())\\\n",
+ " .getPosition(afw_geom.degrees)\n",
+ " ax.text(center[0], center[1], '%d,%d'%(x,y), size=6,\n",
+ " ha=\"center\", va=\"center\")\n",
+ "\n",
+ " skyPosList = [wcs.pixelToSky(pos).getPosition(afw_geom.degrees)\n",
+ " for pos in tractPosList]\n",
+ " ax.set_xlim(max(coord[0] for coord in skyPosList) + 1,\n",
+ " min(coord[0] for coord in skyPosList) - 1)\n",
+ " ax.set_ylim(min(coord[1] for coord in skyPosList) - 1,\n",
+ " max(coord[1] for coord in skyPosList) + 1)\n",
+ " ax.grid(ls=':',color='gray')\n",
+ " ax.set_xlabel(\"RA (deg.)\")\n",
+ " ax.set_ylabel(\"Dec (deg.)\")\n",
+ " ax.set_title(title)\n",
+ " return ax"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The following function loops over the available calexps as returned by the data butler to determine which CCDs to draw. Unfortunately, looping over those calexps using the butler is rather slow, but would be necessary if we wanted to access CCD-level information, like the PSF, from the calexps.\n",
+ "We include this function here for your edification, but then provide a faster function below."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def plot_focal_plane(butler, visit, ax, color='red'):\n",
+ " \"\"\"\n",
+ " Plot the CCDs in the LSST focal plane using the coordinate information in the calexps.\n",
+ " \n",
+ " Notes\n",
+ " -----\n",
+ " By looping over the available calexps, we only plot the CCDs for which image data\n",
+ " are available.\n",
+ " \n",
+ " Parameters\n",
+ " ----------\n",
+ " butler: lsst.daf.persistence.Butler\n",
+ " The data butler serving up data from the desired repo.\n",
+ " visit: int\n",
+ " The visit or obsHistID number.\n",
+ " ax: matplotlib.axes._subplots.AxesSubplot\n",
+ " The matplotlib subplot object onto which to plot the focal plane.\n",
+ " color: str ['red']\n",
+ " Color to use for plotting the individual CCDs.\n",
+ " \n",
+ " Returns\n",
+ " -------\n",
+ " matplotlib.axes._subplots.AxesSubplot: The subplot object used for plotting.\n",
+ " \"\"\"\n",
+ " # We use the `subset` method to obtain all of the `datarefs` (i.e., references to calexp\n",
+ " # data in this case) that satisfy an \"incomplete\" dataId. For visit-level calexp data,\n",
+ " # a unique dataset would specify visit, raft, and sensor. If we just give the visit, then\n",
+ " # references to the available data for all of the CCDs would be returned.\n",
+ " dataId = dict(visit=visit)\n",
+ " datarefs = list(butler.subset('calexp', dataId=dataid))\n",
+ " for i, dataref in enumerate(datarefs):\n",
+ " calexp = dataref.get('calexp')\n",
+ " # We're not going to do anything with it here, but we can get the PSF from the calexp\n",
+ " # like this:\n",
+ " # psf = calexp.getPsf()\n",
+ " # and we can get the zero-point (in ADU) like this\n",
+ " # zero_point = calexp.getCalib().getFluxMag0()\n",
+ " ccd_box = afw_geom.Box2D(calexp.getBBox())\n",
+ " wcs = calexp.getWcs()\n",
+ " path = make_patch(ccd_box.getCorners(), wcs)\n",
+ " ccd = patches.PathPatch(path, alpha=0.2, lw=1, color=color)\n",
+ " ax.add_patch(ccd)\n",
+ " center = wcs.pixelToSky(ccd_box.getCenter()).getPosition(afw_geom.degrees)\n",
+ " return ax"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The following fast version of the focal plane plotting code uses the lsst_sims package to obtain the location and orientation of the CCDs based on the pointing information for the desired visit. That pointing information is extracted from the dithered minion_1016 OpSim db that has been prepared for DC2. Since this code does not access the individual calexps for each CCD, it runs much faster. However, it assumes that the obs_lsstSim package was used in the analysis of the data, and it needs to use the inferred locations of the calexp files to determine if calexp data for a given CCD is available. "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def plot_focal_plane_fast(butler, visit, ax, color='red', opsimdb=None):\n",
+ " \"\"\"\n",
+ " Plot the CCDs in the LSST focal plane using CCD coordinates derived from the pointing\n",
+ " info using the lsst.sims code. \n",
+ " \n",
+ " Notes\n",
+ " -----\n",
+ " This function assumes that the obs_lsstSims package was used to define the camera geometry \n",
+ " for the analysis of the simulated image data.\n",
+ "\n",
+ " Parameters\n",
+ " ----------\n",
+ " butler: lsst.daf.persistence.Butler\n",
+ " The data butler serving up data from the desired repo.\n",
+ " visit: int\n",
+ " The visit or obsHistID number.\n",
+ " ax: matplotlib.axes._subplots.AxesSubplot\n",
+ " The matplotlib subplot object onto which to plot the focal plane.\n",
+ " color: str ['red']\n",
+ " Color to use for plotting the individual CCDs.\n",
+ " opsimDb: str [None]\n",
+ " Filename of the OpSim sqlite database. If None, then the dithered opsim db for Run1.1p\n",
+ " is used.\n",
+ "\n",
+ " Returns\n",
+ " -------\n",
+ " matplotlib.axes._subplots.AxesSubplot: The subplot object used for plotting.\n",
+ " \"\"\"\n",
+ " if opsimdb is None:\n",
+ " opsimdb = '/global/projecta/projectdirs/lsst/groups/SSim/DC2/minion_1016_desc_dithered_v4.db'\n",
+ " conn = sqlite3.connect(opsimdb)\n",
+ " obs_gen = ObservationMetaDataGenerator(database=opsimdb, driver='sqlite')\n",
+ "\n",
+ " # The dithered pointing info was added to the baseline minion_1016 db. We query for the values\n",
+ " # used for the desired visit.\n",
+ " curs = conn.execute('''select descDitheredRA, descDitheredDec, descDitheredRotTelPos\n",
+ " from summary where obshistid={}'''.format(visit))\n",
+ " ra, dec, rottelpos = [np.degrees(x) for x in curs][0]\n",
+ " \n",
+ " # An ObservationMetaData object used to pass the pointing info to the function in\n",
+ " # lsst.sims.coordUtils that provides the CCD coordinates.\n",
+ " obs_md = obs_gen.getObservationMetaData(obsHistID=visit, boundType='circle', boundLength=0.1)[0]\n",
+ " obs_md.pointingRA = ra\n",
+ " obs_md.pointingDec = dec\n",
+ " obs_md.OpsimMetaData['rotTelPos'] = rottelpos\n",
+ "\n",
+ " # Convert the rotation angle of the sky relative to the telescope to the sky angle relative to\n",
+ " # the camera.\n",
+ " obs_md.rotSkyPos = getRotSkyPos(ra, dec, obs_md, rottelpos)\n",
+ " \n",
+ " # Use the butler to get the camera appropriate for this observation. If the data were from a\n",
+ " # different camera, e.g., DECam or HSC, the corresponding camera objects with the associated\n",
+ " # CCD geometries would be returned.\n",
+ " camera = butler.get('camera')\n",
+ " \n",
+ " # Grab one of the calexps via its dataref so that we can ask for its filename and thereby infer\n",
+ " # the location on disk of all of the calexps for this visit.\n",
+ " dataref = list(butler.subset('calexp', visit=visit))[0]\n",
+ " calexp_path = os.path.dirname(os.path.dirname(dataref.get('calexp_filename')[0]))\n",
+ " \n",
+ " # The following code is specific to the obs_lsstSim package and how it names CCDs\n",
+ " # (e.g., \"R:2,2 S:1,1\") and formulates the path components for writing to disk. This\n",
+ " # code would not work for a different obs_ package/camera implementation.\n",
+ " \n",
+ " # Re-order the CCD vertex list returned by the lsst_sims code so that a rectangle is plotted.\n",
+ " corner_index = (np.array([0, 1, 3, 2]),)\n",
+ " for det in camera:\n",
+ " # Skip the WAVEFRONT and GUIDER CCDs\n",
+ " if det.getType() != cameraGeom.SCIENCE:\n",
+ " continue\n",
+ " detname = det.getName()\n",
+ " raft, sensor = re.match(r'R:?(\\d,?\\d)[_ ]S:?(\\d,?\\d)', detname).groups()\n",
+ " raft = 'R' + raft.replace(',', '')\n",
+ " sensor = 'S{}.fits'.format(sensor.replace(',', ''))\n",
+ " if os.path.isfile(os.path.join(calexp_path, raft, sensor)):\n",
+ " corners = np.array(lsst.sims.coordUtils.getCornerRaDec(detname, camera, obs_md))\n",
+ " path = make_patch(corners[corner_index])\n",
+ " ccd = patches.PathPatch(path, alpha=0.2, lw=1, color=color)\n",
+ " ax.add_patch(ccd)\n",
+ " \n",
+ " return ax"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "The following function just plots the boundaries of the Run1.1p regions as described in the [Run 1.1p Specifications document](https://docs.google.com/document/d/1aQOPL9smeDlhtlwDrp39Zuu2q8DKivDaHLQX3_omwOI/edit)."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def plot_Run1_1p_region(ax):\n",
+ " \"\"\"Plot the WFD and uDDF regions for Run1.1p.\"\"\"\n",
+ " uddf_ra = [53.764, 52.486, 52.479, 53.771, 53.764]\n",
+ " uddf_dec = [-27.533, -27.533, -28.667, -28.667, -27.533]\n",
+ "\n",
+ " wfd_ra = [52.25, 52.11, 58.02, 57.87, 52.25]\n",
+ " wfd_dec = [-27.25, -32.25, -32.25, -27.25, -27.25]\n",
+ "\n",
+ " ax.errorbar(wfd_ra, wfd_dec, fmt='-', color='green',\n",
+ " label='protoDC2 WFD boundary')\n",
+ " ax.errorbar(uddf_ra, uddf_dec, fmt='-', color='red',\n",
+ " label='uDDF boundary')\n",
+ " return ax"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Making the Plot\n",
+ "We now use the above functions to plot the tracts and patches on the sky. "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "repo = REPOS['1.2p']\n",
+ "butler = dp.Butler(repo)\n",
+ "\n",
+ "# Single visit to overlay:\n",
+ "visit = 219976\n",
+ "\n",
+ "# WARNING: this is a hack to search the coadd folder for the tracts that have data.\n",
+ "# Unfortunately, this information is not directly accessible from the data butler. \n",
+ "ref_path = butler.getUri('deepCoadd_forced_src', tract=4851, patch='0,0', filter='i')\n",
+ "ref_path, success, _ = ref_path.partition('/4851/0,0')\n",
+ "if not success:\n",
+ " ref_path = ref_path.partition('/i_t4851_p0,0')[0]\n",
+ "tracts = sorted(set((int(re.search(r'(\\d{4})', t).groups()[0]) for t in os.listdir(ref_path) if os.path.isdir(os.path.join(ref_path, t)))))\n",
+ "\n",
+ "# First, loop over all the tracts, plotting them as gray, numbered, rectangles:\n",
+ "ax = None\n",
+ "skyMap = butler.get('deepCoadd_skyMap')\n",
+ "for tract in tracts:\n",
+ " ax = plot_skymap_tract(skyMap, tract=tract, title='', ax=ax)\n",
+ "\n",
+ "# Now overlay a single focal plane, for the chosen visit, in violet:\n",
+ "plot_focal_plane_fast(butler, visit, ax, color=\"violet\")\n",
+ "\n",
+ "# Finally, overlay the Run 1.1p main survey and uDDF regions, and then adjust the axes:\n",
+ "plot_Run1_1p_region(ax)\n",
+ "\n",
+ "ax.set_xlim(60.5, 50.0)\n",
+ "ax.set_ylim(-33.5, -25.5)\n",
+ "plt.legend(loc=0);"
+ ]
+ }
+ ],
+ "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.6.2"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/Basics/dm_butler_skymap.py b/Basics/dm_butler_skymap.py
new file mode 100644
index 00000000..0827f88e
--- /dev/null
+++ b/Basics/dm_butler_skymap.py
@@ -0,0 +1,341 @@
+#!/usr/bin/env python
+# coding: utf-8
+
+# # Plotting the DC2 Run1.1p skyMap
+#
Owner: **Jim Chiang** ([@jchiang87](https://github.com/LSSTDESC/DC2-analysis/issues/new?body=@jchiang87))
+#
Last Verified to Run: **2018-10-26** (by @yymao)
+#
+# In this notebook, we show how to use the data butler to obtain information on the skyMap used in the coadd analyses performed by the DRP pipeline. These skyMaps are composed of tracts and patches on the sky. Each tract is a rectangular region of the sky with a common map projection; each tract is further divided into rectangular patches, which use the same tract coordinate system and which are a convenient size for processing the coadd data. A more complete description of the skyMap geometry is given in the HSC Software Pipeline paper ([Bosch et al. 2017](https://arxiv.org/abs/1705.06766)).
+#
+# Equipped with the info from the skyMap, we plot the tracts and patches that were used with the Run1.1p processing and overlay the WFD and uDDF simulation regions. We also use the butler to access the visit-level data and show how one can access the calexp (calibrated exposure) image data to obtain the PSF, zero-point, etc. as measured by the Stack for a given exposure. Finally, we show how to plot the sky region imaged on the focal plane for a given visit in two ways: the first using the CCD coordinates available from the calexps and the other using the lsst_sims code to compute those coordinates from the pointing information for the visit.
+#
+# ## Set Up
+
+# In[ ]:
+
+
+import os
+import glob
+import warnings
+import sqlite3
+import re
+
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.path import Path
+import matplotlib.patches as patches
+get_ipython().run_line_magic('matplotlib', 'inline')
+
+import lsst.afw.geom as afw_geom
+import lsst.afw.cameraGeom as cameraGeom
+import lsst.daf.persistence as dp
+# The lsst_sims code issues some ignorable warnings regarding ids used for querying the object
+# databases.
+
+#with warnings.catch_warnings():
+# warnings.simplefilter("ignore")
+# import lsst.sims.coordUtils
+# from lsst.sims.catUtils.utils import ObservationMetaDataGenerator
+# from lsst.sims.utils import getRotSkyPos
+
+#from desc_dc2_dm_data import REPOS
+
+
+# ## Plotting Functions
+# First, we define some functions to plot the tract, patch, and CCD regions on the sky. These are copied from [example code](https://github.com/yalsayyad/dm_notebooks/blob/master/desc-ssim/DESC-SSim%20Patch%20Geometry.ipynb) that Yusra AlSayyad presented at the [2017-06-29 SSim meeting](https://confluence.slac.stanford.edu/pages/viewpage.action?pageId=224461017).
+
+# In[ ]:
+
+
+def make_patch(vertexList, wcs=None):
+ """
+ Return a Path in sky coords from vertex list in pixel coords.
+
+ Parameters
+ ----------
+ vertexList: list of coordinates
+ These are the corners of the region to be plotted either in pixel coordinates or
+ sky coordinates.
+ wcs: lsst.afw.geom.skyWcs.skyWcs.SkyWcs [None]
+ The WCS object used to convert from pixel to sky coordinates.
+
+ Returns
+ -------
+ matplotlib.path.Path: The encapsulation of the vertex info that matplotlib uses to
+ plot a patch.
+ """
+ if wcs is not None:
+ skyPatchList = [wcs.pixelToSky(pos).getPosition(afw_geom.degrees)
+ for pos in vertexList]
+ else:
+ skyPatchList = vertexList
+ verts = [(coord[0], coord[1]) for coord in skyPatchList]
+ verts.append((0,0))
+ codes = [Path.MOVETO,
+ Path.LINETO,
+ Path.LINETO,
+ Path.LINETO,
+ Path.CLOSEPOLY,
+ ]
+ return Path(verts, codes)
+
+
+# In[ ]:
+
+
+def plot_skymap_tract(skyMap, tract=0, title=None, ax=None):
+ """
+ Plot a tract from a skyMap.
+
+ Parameters
+ ----------
+ skyMap: lsst.skyMap.SkyMap
+ The SkyMap object containing the tract and patch information.
+ tract: int [0]
+ The tract id of the desired tract to plot.
+ title: str [None]
+ Title of the tract plot. If None, the use `tract `.
+ ax: matplotlib.axes._subplots.AxesSubplot [None]
+ The subplot object to contain the tract plot. If None, then make a new one.
+
+ Returns
+ -------
+ matplotlib.axes._subplots.AxesSubplot: The subplot containing the tract plot.
+ """
+ if title is None:
+ title = 'tract {}'.format(tract)
+ tractInfo = skyMap[tract]
+ tractBox = afw_geom.Box2D(tractInfo.getBBox())
+ tractPosList = tractBox.getCorners()
+ wcs = tractInfo.getWcs()
+ xNum, yNum = tractInfo.getNumPatches()
+
+ if ax is None:
+ fig = plt.figure(figsize=(12,8))
+ ax = fig.add_subplot(111)
+
+ tract_center = wcs.pixelToSky(tractBox.getCenter()) .getPosition(afw_geom.degrees)
+ ax.text(tract_center[0], tract_center[1], '%d' % tract, size=16,
+ ha="center", va="center", color='blue')
+ for x in range(xNum):
+ for y in range(yNum):
+ patchInfo = tractInfo.getPatchInfo([x, y])
+ patchBox = afw_geom.Box2D(patchInfo.getOuterBBox())
+ pixelPatchList = patchBox.getCorners()
+ path = make_patch(pixelPatchList, wcs)
+ patch = patches.PathPatch(path, alpha=0.1, lw=1)
+ ax.add_patch(patch)
+ center = wcs.pixelToSky(patchBox.getCenter()) .getPosition(afw_geom.degrees)
+ ax.text(center[0], center[1], '%d,%d'%(x,y), size=6,
+ ha="center", va="center")
+
+ skyPosList = [wcs.pixelToSky(pos).getPosition(afw_geom.degrees)
+ for pos in tractPosList]
+ ax.set_xlim(max(coord[0] for coord in skyPosList) + 1,
+ min(coord[0] for coord in skyPosList) - 1)
+ ax.set_ylim(min(coord[1] for coord in skyPosList) - 1,
+ max(coord[1] for coord in skyPosList) + 1)
+ ax.grid(ls=':',color='gray')
+ ax.set_xlabel("RA (deg.)")
+ ax.set_ylabel("Dec (deg.)")
+ ax.set_title(title)
+ return ax
+
+
+# The following function loops over the available calexps as returned by the data butler to determine which CCDs to draw. Unfortunately, looping over those calexps using the butler is rather slow, but would be necessary if we wanted to access CCD-level information, like the PSF, from the calexps.
+# We include this function here for your edification, but then provide a faster function below.
+
+# In[ ]:
+
+
+def plot_focal_plane(butler, visit, ax, color='red'):
+ """
+ Plot the CCDs in the LSST focal plane using the coordinate information in the calexps.
+
+ Notes
+ -----
+ By looping over the available calexps, we only plot the CCDs for which image data
+ are available.
+
+ Parameters
+ ----------
+ butler: lsst.daf.persistence.Butler
+ The data butler serving up data from the desired repo.
+ visit: int
+ The visit or obsHistID number.
+ ax: matplotlib.axes._subplots.AxesSubplot
+ The matplotlib subplot object onto which to plot the focal plane.
+ color: str ['red']
+ Color to use for plotting the individual CCDs.
+
+ Returns
+ -------
+ matplotlib.axes._subplots.AxesSubplot: The subplot object used for plotting.
+ """
+ # We use the `subset` method to obtain all of the `datarefs` (i.e., references to calexp
+ # data in this case) that satisfy an "incomplete" dataId. For visit-level calexp data,
+ # a unique dataset would specify visit, raft, and sensor. If we just give the visit, then
+ # references to the available data for all of the CCDs would be returned.
+ dataId = dict(visit=visit)
+ datarefs = list(butler.subset('calexp', dataId=dataid))
+ for i, dataref in enumerate(datarefs):
+ calexp = dataref.get('calexp')
+ # We're not going to do anything with it here, but we can get the PSF from the calexp
+ # like this:
+ # psf = calexp.getPsf()
+ # and we can get the zero-point (in ADU) like this
+ # zero_point = calexp.getCalib().getFluxMag0()
+ ccd_box = afw_geom.Box2D(calexp.getBBox())
+ wcs = calexp.getWcs()
+ path = make_patch(ccd_box.getCorners(), wcs)
+ ccd = patches.PathPatch(path, alpha=0.2, lw=1, color=color)
+ ax.add_patch(ccd)
+ center = wcs.pixelToSky(ccd_box.getCenter()).getPosition(afw_geom.degrees)
+ return ax
+
+
+# The following fast version of the focal plane plotting code uses the lsst_sims package to obtain the location and orientation of the CCDs based on the pointing information for the desired visit. That pointing information is extracted from the dithered minion_1016 OpSim db that has been prepared for DC2. Since this code does not access the individual calexps for each CCD, it runs much faster. However, it assumes that the obs_lsstSim package was used in the analysis of the data, and it needs to use the inferred locations of the calexp files to determine if calexp data for a given CCD is available.
+
+# In[ ]:
+
+
+def plot_focal_plane_fast(butler, visit, ax, color='red', opsimdb=None):
+ """
+ Plot the CCDs in the LSST focal plane using CCD coordinates derived from the pointing
+ info using the lsst.sims code.
+
+ Notes
+ -----
+ This function assumes that the obs_lsstSims package was used to define the camera geometry
+ for the analysis of the simulated image data.
+
+ Parameters
+ ----------
+ butler: lsst.daf.persistence.Butler
+ The data butler serving up data from the desired repo.
+ visit: int
+ The visit or obsHistID number.
+ ax: matplotlib.axes._subplots.AxesSubplot
+ The matplotlib subplot object onto which to plot the focal plane.
+ color: str ['red']
+ Color to use for plotting the individual CCDs.
+ opsimDb: str [None]
+ Filename of the OpSim sqlite database. If None, then the dithered opsim db for Run1.1p
+ is used.
+
+ Returns
+ -------
+ matplotlib.axes._subplots.AxesSubplot: The subplot object used for plotting.
+ """
+ if opsimdb is None:
+ opsimdb = '/global/projecta/projectdirs/lsst/groups/SSim/DC2/minion_1016_desc_dithered_v4.db'
+ conn = sqlite3.connect(opsimdb)
+ obs_gen = ObservationMetaDataGenerator(database=opsimdb, driver='sqlite')
+
+ # The dithered pointing info was added to the baseline minion_1016 db. We query for the values
+ # used for the desired visit.
+ curs = conn.execute('''select descDitheredRA, descDitheredDec, descDitheredRotTelPos
+ from summary where obshistid={}'''.format(visit))
+ ra, dec, rottelpos = [np.degrees(x) for x in curs][0]
+
+ # An ObservationMetaData object used to pass the pointing info to the function in
+ # lsst.sims.coordUtils that provides the CCD coordinates.
+ obs_md = obs_gen.getObservationMetaData(obsHistID=visit, boundType='circle', boundLength=0.1)[0]
+ obs_md.pointingRA = ra
+ obs_md.pointingDec = dec
+ obs_md.OpsimMetaData['rotTelPos'] = rottelpos
+
+ # Convert the rotation angle of the sky relative to the telescope to the sky angle relative to
+ # the camera.
+ obs_md.rotSkyPos = getRotSkyPos(ra, dec, obs_md, rottelpos)
+
+ # Use the butler to get the camera appropriate for this observation. If the data were from a
+ # different camera, e.g., DECam or HSC, the corresponding camera objects with the associated
+ # CCD geometries would be returned.
+ camera = butler.get('camera')
+
+ # Grab one of the calexps via its dataref so that we can ask for its filename and thereby infer
+ # the location on disk of all of the calexps for this visit.
+ dataref = list(butler.subset('calexp', visit=visit))[0]
+ calexp_path = os.path.dirname(os.path.dirname(dataref.get('calexp_filename')[0]))
+
+ # The following code is specific to the obs_lsstSim package and how it names CCDs
+ # (e.g., "R:2,2 S:1,1") and formulates the path components for writing to disk. This
+ # code would not work for a different obs_ package/camera implementation.
+
+ # Re-order the CCD vertex list returned by the lsst_sims code so that a rectangle is plotted.
+ corner_index = (np.array([0, 1, 3, 2]),)
+ for det in camera:
+ # Skip the WAVEFRONT and GUIDER CCDs
+ if det.getType() != cameraGeom.SCIENCE:
+ continue
+ detname = det.getName()
+ raft, sensor = re.match(r'R:?(\d,?\d)[_ ]S:?(\d,?\d)', detname).groups()
+ raft = 'R' + raft.replace(',', '')
+ sensor = 'S{}.fits'.format(sensor.replace(',', ''))
+ if os.path.isfile(os.path.join(calexp_path, raft, sensor)):
+ corners = np.array(lsst.sims.coordUtils.getCornerRaDec(detname, camera, obs_md))
+ path = make_patch(corners[corner_index])
+ ccd = patches.PathPatch(path, alpha=0.2, lw=1, color=color)
+ ax.add_patch(ccd)
+
+ return ax
+
+
+# The following function just plots the boundaries of the Run1.1p regions as described in the [Run 1.1p Specifications document](https://docs.google.com/document/d/1aQOPL9smeDlhtlwDrp39Zuu2q8DKivDaHLQX3_omwOI/edit).
+
+# In[ ]:
+
+
+def plot_Run1_1p_region(ax):
+ """Plot the WFD and uDDF regions for Run1.1p."""
+ uddf_ra = [53.764, 52.486, 52.479, 53.771, 53.764]
+ uddf_dec = [-27.533, -27.533, -28.667, -28.667, -27.533]
+
+ wfd_ra = [52.25, 52.11, 58.02, 57.87, 52.25]
+ wfd_dec = [-27.25, -32.25, -32.25, -27.25, -27.25]
+
+ ax.errorbar(wfd_ra, wfd_dec, fmt='-', color='green',
+ label='protoDC2 WFD boundary')
+ ax.errorbar(uddf_ra, uddf_dec, fmt='-', color='red',
+ label='uDDF boundary')
+ return ax
+
+
+# ## Making the Plot
+# We now use the above functions to plot the tracts and patches on the sky.
+
+# In[ ]:
+
+"""
+repo = REPOS['1.2p']
+butler = dp.Butler(repo)
+
+# Single visit to overlay:
+visit = 219976
+
+# WARNING: this is a hack to search the coadd folder for the tracts that have data.
+# Unfortunately, this information is not directly accessible from the data butler.
+ref_path = butler.getUri('deepCoadd_forced_src', tract=4851, patch='0,0', filter='i')
+ref_path, success, _ = ref_path.partition('/4851/0,0')
+if not success:
+ ref_path = ref_path.partition('/i_t4851_p0,0')[0]
+tracts = sorted(set((int(re.search(r'(\d{4})', t).groups()[0]) for t in os.listdir(ref_path) if os.path.isdir(os.path.join(ref_path, t)))))
+
+# First, loop over all the tracts, plotting them as gray, numbered, rectangles:
+ax = None
+skyMap = butler.get('deepCoadd_skyMap')
+for tract in tracts:
+ ax = plot_skymap_tract(skyMap, tract=tract, title='', ax=ax)
+
+# Now overlay a single focal plane, for the chosen visit, in violet:
+plot_focal_plane_fast(butler, visit, ax, color="violet")
+
+# Finally, overlay the Run 1.1p main survey and uDDF regions, and then adjust the axes:
+plot_Run1_1p_region(ax)
+
+ax.set_xlim(60.5, 50.0)
+ax.set_ylim(-33.5, -25.5)
+plt.legend(loc=0);
+
+"""
\ No newline at end of file
diff --git a/GettingStarted/templates/template_Notebook.ipynb b/GettingStarted/templates/template_Notebook.ipynb
index db8831b9..70eb1e62 100644
--- a/GettingStarted/templates/template_Notebook.ipynb
+++ b/GettingStarted/templates/template_Notebook.ipynb
@@ -47,17 +47,9 @@
},
{
"cell_type": "code",
- "execution_count": 1,
+ "execution_count": null,
"metadata": {},
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "python: can't open file 'setup.py': [Errno 2] No such file or directory\n"
- ]
- }
- ],
+ "outputs": [],
"source": [
"! cd .. && python setup.py -q develop --user && cd -"
]
@@ -71,7 +63,7 @@
},
{
"cell_type": "code",
- "execution_count": 2,
+ "execution_count": null,
"metadata": {},
"outputs": [],
"source": [
@@ -125,7 +117,12 @@
"import numpy as np\n",
"import matplotlib as mpl\n",
"import matplotlib.pyplot as plt\n",
- "from IPython.display import IFrame, display, Markdown"
+ "from IPython.display import IFrame, display, Markdown\n",
+ "import warnings\n",
+ "\n",
+ "# Filter some warnings printed by v16.0 of the stack\n",
+ "warnings.simplefilter(\"ignore\", category=FutureWarning)\n",
+ "warnings.simplefilter(\"ignore\", category=UserWarning)"
]
},
{
diff --git a/ImageProcessing/Re-RunHSC.ipynb b/ImageProcessing/Re-RunHSC.ipynb
index c6e77690..8ea9b216 100644
--- a/ImageProcessing/Re-RunHSC.ipynb
+++ b/ImageProcessing/Re-RunHSC.ipynb
@@ -6,13 +6,13 @@
"source": [
"# HSC Re-Run: Making Forced Photometry Light Curves from Scratch\n",
"\n",
- "
Owner: **Justin Myles** ([@jtmyles](https://github.com/LSSTScienceCollaborations/StackClub/issues/new?body=@jtmyles))\n",
+ "
Owners: **Justin Myles** ([@jtmyles](https://github.com/LSSTScienceCollaborations/StackClub/issues/new?body=@jtmyles)), **Phil Marshall** ([@drphilmarshall](https://github.com/LSSTScienceCollaborations/StackClub/issues/new?body=@drphilmarshall))\n",
"
Last Verified to Run: **2018-09-21**\n",
"
Verified Stack Release: **16.0**\n",
"\n",
"This project addresses issue [#63: HSC Re-run](https://github.com/LSSTScienceCollaborations/StackClub/issues/63)\n",
"\n",
- "This notebook demonstrates the pipeline described in the [LSST Science Piplines data processing tutorial](https://pipelines.lsst.io/), from ingesting images (using the [obs_subaru](https://github.com/lsst/obs_subaru) package) through image processing, coaddition, source detection and object measurement all the way through to measuring forced photometry light curves in a small patch of the HSC sky (in the [ci_hsc](https://github.com/lsst/ci_hsc/) repository). \n",
+ "This notebook demonstrates the pipeline described in the [LSST Science Piplines data processing tutorial](https://pipelines.lsst.io/), from ingesting images (using the [obs_subaru](https://github.com/lsst/obs_subaru) package) through image processing, coaddition, source detection and object measurement all the way through to measuring forced photometry light curves in a small patch of the HSC sky (in the [ci_hsc](https://github.com/lsst/ci_hsc/) repository). It does this by calling a `bash` script, having first identified a minimal data set for demonstration purposes. \n",
"\n",
"### Learning Objectives:\n",
"After working through and studying this notebook you should be able to understand how to use the DRP pipeline from image visualization through to a forced photometry light curve. Specific learning objectives include: \n",
@@ -63,7 +63,7 @@
"source": [
"## Pipeline Preview\n",
"\n",
- "Before we unpack the pipeline described in the [LSST Science Piplines data processing tutorial](https://pipelines.lsst.io/), let's look at the complete set of command line tasks assembled into an end-to-end data reduction script."
+ "The pipeline described in the [LSST Science Piplines data processing tutorial](https://pipelines.lsst.io/) contains a complete set of command line tasks that can be assembled into an end-to-end data reduction pipeline script. Let's see what this script looks like."
]
},
{
@@ -175,7 +175,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "# Part 2: Calibrating single frames\n",
+ "## Part 2: Calibrating single frames\n",
"https://pipelines.lsst.io/getting-started/processccd.html\n",
"\n",
"Part II runs the following command-line task:\n",
diff --git a/stackclub/__init__.py b/stackclub/__init__.py
index 6c9efb13..78552e31 100644
--- a/stackclub/__init__.py
+++ b/stackclub/__init__.py
@@ -1,3 +1,4 @@
from .where_is import *
from .wimport import *
-from .nbimport import *
\ No newline at end of file
+from .nbimport import *
+from .taster import *
\ No newline at end of file
diff --git a/stackclub/taster.py b/stackclub/taster.py
new file mode 100644
index 00000000..c3692c3d
--- /dev/null
+++ b/stackclub/taster.py
@@ -0,0 +1,258 @@
+import numpy as np
+from IPython.display import display, Markdown
+
+class Taster(object):
+ """
+ Worker for tasting the datasets in a Butler's repo (based mostly off of querying metadata).
+ Instantiate with a repo.
+ """
+ def __init__(self, repo, vb=False, path_to_tracts=''):
+ self.repo = repo
+ # Instantiate a butler, or report failure:
+ from lsst.daf.persistence import Butler
+ try:
+ self.butler = Butler(repo)
+ except:
+ self.butler = None
+ print("Warning: failed to instantiate a butler to get data from repo '"+repo+"'")
+ return None
+ # Set up some internal variables:
+ self.vb = vb
+ self.exists = {}
+ self.existence = False
+ self.counts = {}
+ self.tracts = []
+ self.path_to_tracts = path_to_tracts
+ if path_to_tracts != '':
+ try:
+ self.skymap_butler = Butler(repo + path_to_tracts)
+ except:
+ self.skymap_butler = None
+ print("Warning: failed to find a skyMap for the path " + repo + path_to_tracts)
+ return
+
+ def what_exists(self, all=False):
+ """
+ Check for the existence of various useful things.
+
+ Parameters
+ ==========
+ all: boolean
+ If true, the method will check all possible dataset types
+
+ Returns
+ =======
+ exists: dict
+ Checklist of what exists (True) and what does not (False)
+ """
+ # Get mappers for all tested repos
+ from lsst.obs.hsc import HscMapper
+ from lsst.obs.comCam import ComCamMapper
+ #from lsst.obs.lsst import LsstCamMapper
+ from lsst.obs.ctio0m9 import Ctio0m9Mapper
+
+ #select proper mapper
+ if self.repo.find('hsc') != -1: mapper = HscMapper(root=self.repo)
+ elif self.repo.find('comCam') != -1: mapper = ComCamMapper(root=self.repo)
+ #elif self.repo.find('DC2') != -1: mapper = LsstCamMapper(root=self.repo)
+ elif self.repo.find('ctio0m9') != -1: mapper = Ctio0m9Mapper(root=self.repo)
+ else: print("Unable to locate Mapper file in specified repo. Check that you selected a valid repo.")
+
+
+ if all:
+ #collect a list of all possible dataset types
+ mapper = HscMapper(root=self.repo)
+ all_dataset_types = mapper.getDatasetTypes()
+
+ remove = ['_config', '_filename', '_md', '_sub', '_len', '_schema', '_metadata']
+
+ interesting = []
+ for dataset_type in all_dataset_types:
+ keep = True
+ for word in remove:
+ if word in dataset_type:
+ keep = False
+ if keep:
+ interesting.append(dataset_type)
+
+ else:
+ interesting = ['raw', 'calexp', 'src', 'deepCoadd_calexp', 'deepCoadd_meas']
+
+ self.look_for_datasets_of_type(interesting)
+ self.look_for_skymap()
+ self.existence = True
+ return
+
+ def look_for_datasets_of_type(self, datasettypes):
+ """
+ Check whether dataset of given type is in the metadata.
+
+ Parameters
+ ==========
+ datasettype: list of strings
+ Types of dataset to check for, eg 'calexp', 'raw', 'wcs' etc.
+ """
+ datasets_that_exist = []
+ datasets_that_do_not_exist = []
+
+ for datasettype in datasettypes:
+ try:
+ datasetkeys = self.butler.getKeys(datasettype)
+ onekey = list(datasetkeys.keys())[0]
+ metadata = self.butler.queryMetadata(datasettype, [onekey])
+ #if self.vb: print("{} dataset exists.".format(datasettype))
+ datasets_that_exist.append(datasettype)
+ self.exists[datasettype] = True
+ except:
+ #if self.vb: print("{} dataset doesn't exist.".format(datasettype))
+ datasets_that_do_not_exist.append(datasettype)
+ self.exists[datasettype] = False
+
+ #Organize output
+ if self.vb:
+ print("Datasets that exist\n-------------------")
+ print(datasets_that_exist)
+ print("\nDatasets that do not exist\n--------------------------")
+ print(datasets_that_do_not_exist)
+
+ return
+
+ def look_for_skymap(self):
+ """
+ Check for the existence of a skymap.
+ """
+ try:
+ self.skyMap = self.skymap_butler.get('deepCoadd_skyMap')
+ self.exists['deepCoadd_skyMap'] = True
+ if self.vb: print("\nSkymap\n-------------------\ndeepCoadd_skyMap exists.")
+ except:
+ self.skyMap = None
+ self.exists['deepCoadd_skyMap'] = False
+ if self.vb: print("\nSkymap\n-------------------\ndeepCoadd_skyMap doesn't exist.")
+ return
+
+
+
+ def estimate_sky_area(self):
+ """
+ Use available skymap to estimate sky area covered by tracts and patches.
+
+ Returns
+ =======
+ area: float
+ Sky area in square degrees
+ """
+ if self.skyMap is None: return None
+
+ area_label = 'Total Sky Area (deg$^2$)'
+ if area_label in self.counts.keys():
+ return self.counts[area_label]
+
+ # Collect tracts from files
+ import os, glob
+ tracts = sorted([int(os.path.basename(x)) for x in
+ glob.glob(os.path.join(self.repo + self.path_to_tracts, 'deepCoadd-results', 'merged', '*'))])
+
+ self.tracts = tracts
+ self.counts['Number of Tracts'] = len(tracts)
+
+ # Note: We'd like to do this with the butler, but it appears 'tracts' have to be
+ # specified in the dataId to be queried, so the queryMetadata method fails
+
+ # Calculate area from all tracts
+ total_area = 0.0 #deg^2
+ plotting_vertices = []
+ for test_tract in tracts:
+ # Get inner vertices for tract
+ tractInfo = self.skyMap[test_tract]
+ vertices = tractInfo._vertexCoordList
+ plotting_vertices.append(vertices)
+
+ # Calculate area of box
+ av_dec = 0.5 * (vertices[2][1] + vertices[0][1])
+ av_dec = av_dec.asRadians()
+ delta_ra_raw = vertices[0][0] - vertices[1][0]
+ delta_ra = delta_ra_raw.asDegrees() * np.cos(av_dec)
+ delta_dec= vertices[2][1] - vertices[0][1]
+ area = delta_ra * delta_dec.asDegrees()
+
+ # Combine areas
+ total_area += area
+
+ if self.vb: print(area_label, ": ", total_area)
+
+ # Round of the total area for table purposes
+ self.counts[area_label] = round(total_area, 2)
+ return self.counts[area_label]
+
+ def count_things(self):
+ """
+ Count the available number of calexp visits, sensors, fields etc.
+ """
+ # Collect numbers of images of various kinds:
+ if self.exists['calexp']:
+ self.counts['Number of Visits'] = \
+ len(self.butler.queryMetadata('calexp', ['visit']))
+ self.counts['Number of Pointings'] = \
+ len(self.butler.queryMetadata('calexp', ['pointing']))
+ self.counts['Number of Sensor Visits'] = \
+ len(self.butler.queryMetadata('calexp', ['ccd']))
+ self.counts['Number of Fields'] = \
+ len(self.butler.queryMetadata('calexp', ['field']))
+ self.counts['Number of Filters'] = \
+ len(self.butler.queryMetadata('calexp', ['filter']))
+ # Collect number of objects from Source Catalog
+ if self.exists['src']:
+ self.counts['Number of Sources'] = \
+ len(self.butler.queryMetadata('src', ['id']))
+ return
+
+ def plot_sky_coverage(self):
+ import matplotlib.pyplot as plt
+ fig = plt.figure()
+
+ for tract in self.tracts:
+ tractInfo = self.skyMap[tract]
+
+ corners = [(x[0].asDegrees(), x[1].asDegrees()) for x in tractInfo.getVertexList()]
+ x = [k[0] for k in corners] + [corners[0][0]]
+ y = [k[1] for k in corners] + [corners[0][1]]
+
+ plt.plot(x,y, color='b')
+
+ plt.xlabel('RA (deg)')
+ plt.ylabel('Dec (deg)')
+ plt.title('2D Projection of Sky Coverage')
+
+ plt.show()
+ return
+
+
+ def report(self):
+ """
+ Print a nice report of the data available in this repo.
+ """
+ # First check what's there:
+ if not self.existence: self.what_exists()
+
+ # Then, get the numbers:
+ self.count_things()
+ self.estimate_sky_area()
+
+ # A nice bold section heading:
+ display(Markdown('### Main Repo: %s' % self.repo))
+ if self.path_to_tracts != '':
+ display(Markdown('### Specified Tract Directory: %s' %self.path_to_tracts))
+
+ # Make a table of the collected metadata
+ output_table = "| Metadata Characteristics | | \n | :---: | --- | \n "
+ for key in self.counts.keys():
+ output_table += "| %s | %s | \n" %(key, self.counts[key])
+
+ # Display it:
+ display(Markdown(output_table))
+
+ # Plot sky coverage
+ self.plot_sky_coverage()
+
+ return
diff --git a/stackclub/where_is.py b/stackclub/where_is.py
index 9ff5ae45..ab8f7ab0 100644
--- a/stackclub/where_is.py
+++ b/stackclub/where_is.py
@@ -1,13 +1,15 @@
-def where_is(object, in_the=None):
+def where_is(object, in_the='source', assuming_its_a=None):
"""
Print a markdown hyperlink to the source code of `object`.
Parameters
----------
- object: python object
- The class or function you are looking for.
+ object: python object or string
+ The class or function you are looking for, or the name of a python object or file.
in_the: string, optional
The kind of place you want to look in: `['source', 'repo', 'technotes']`
+ assuming_its_a: string, optional
+ The kind of object you think you have: `['cmdlinetask'], default=None
Examples
--------
@@ -17,32 +19,42 @@ def where_is(object, in_the=None):
>>> where_is(Butler.get, in_the='source')
>>> where_is(Butler, in_the='repo')
>>> where_is(Butler, in_the='technotes')
+ >>> where_is("makeDiscreteSkyMap.py", in_the="source", assuming_its_a="cmdlinetask")
Notes
-----
See also the `FindingDocs tutorial notebook `_ for a working demo.
"""
+ # Deal with string object names - useful for locating command line tasks:
+ if isinstance(object, str):
+ objectname = object
+ if in_the == 'source' and assuming_its_a == None:
+ raise ValueError('Cannot locate task/object `'+object+'` in the source by name. Either pass in an object, or use the "assuming_its_a" kwarg to guess what kind of object it is.')
+ if assuming_its_a == "cmdlinetask":
+ modulename = 'lsst.pipe.tasks.'+objectname
+
+ elif isinstance(object, module):
+ # Locate the module that contains the desired object, and break its name into pieces:
+ modulename = object.__module__
+ objectname = object.__name__
- # Locate the module that contains the desired object, and break its name into pieces:
- modulename = object.__module__
- pieces = str.split(modulename,'.')
- objectname = object.__name__
+ else:
+ raise TypeError('Expecting objects of type "string" or "module"')
- # Form the URL, and a useful markdown representation of it:
- if in_the is None: in_the = 'source'
-
+ # Form the URL, and a useful markdown representation of it:
if in_the == 'source':
+ pieces = str.split(modulename,'.')
URL = 'https://github.com/'+pieces[0]+'/'+pieces[1]+'_'+pieces[2] \
+ '/blob/master/python/'+pieces[0]+'/'+pieces[1]+'/'+pieces[2]+'/'+pieces[3]+'.py'
- link = '['+modulename+']('+URL+')'
+ link = '[`'+modulename+'`]('+URL+')'
elif in_the == 'repo':
- URL = 'https://github.com/search?l=Python&q=user%3Alsst+'+objectname+'&type=Code'
+ URL = 'https://github.com/search?l=Python&q=org%3Alsst+'+objectname+'&type=Code'
link = '[searching for `'+objectname+'` in the `lsst` repo]('+URL+')'
elif in_the == 'technotes':
- URL = 'https://github.com/search?l=reStructuredText&q=user%3Alsst-dm+'+objectname+'&type=Code'
+ URL = 'https://github.com/search?l=reStructuredText&q=org%3Alsst-dm+'+objectname+'&type=Code'
link = '[searching for `'+objectname+'` in the `lsst-dm` technotes]('+URL+')'
else: