diff --git a/Basics/afw_table_guided_tour.ipynb b/Basics/afw_table_guided_tour.ipynb new file mode 100644 index 00000000..e10da507 --- /dev/null +++ b/Basics/afw_table_guided_tour.ipynb @@ -0,0 +1,1204 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "# afwTables: A Guided Tour\n", + "
Owner(s): **Imran Hasan** ([@ih64](https://github.com/LSSTScienceCollaborations/StackClub/issues/new?body=@ih64))\n", + "
Last Verified to Run: **2018-10-19**\n", + "
Verified Stack Release: **16.0**\n", + "\n", + "Catalogs of astronomical objects, and their many automated measurements, will be the primary data product that LSST provides, and queries of those catalogs will be the starting point for almost all LSST science analyses. On the way to filling the LSST database with these catalogs, the science pipelines will generate and manipulate a lot of internal tables; the python class that the Stack defines and uses for these tables is called an \"afwTable\". \n", + "\n", + "### Learning Objectives:\n", + "\n", + "After working through this tutorial you should be able to: \n", + "1. Make a bare bones afw schema and table;\n", + "2. Set and get values in a schema and table;\n", + "3. Read and write a source detection catalog table;\n", + "4. Learn to use source detection catalog methods, and to avoid common pitfalls;\n", + "5. Learn to use source match vectors.\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", + "\n", + "## Set-up" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The next few cells give you some options for your \"Set-up\" section - you may not need them all." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We'll need the `stackclub` package to be installed. If you are not developing this package, you can install it using `pip`, like this:\n", + "```\n", + "pip install git+git://github.com/LSSTScienceCollaborations/StackClub.git#egg=stackclub\n", + "```\n", + "If you are developing the `stackclub` package (eg by adding modules to it to support the Stack Club tutorial that you are writing, you'll need to make a local, editable installation. In the top level folder of the `StackClub` repo, do:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "! cd .. && python setup.py -q develop --user && cd -" + ] + }, + { + "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." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%load_ext autoreload\n", + "%autoreload 2" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "You can find the Stack version that this notebook is running by using eups list -s on the terminal command line:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# What version of the Stack am I using?\n", + "! echo $HOSTNAME\n", + "! eups list -s | grep lsst_distrib" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "For this tutorial we'll need the following modules:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "#%matplotlib ipympl\n", + "\n", + "import os\n", + "import warnings\n", + "import numpy as np\n", + "import matplotlib as mpl\n", + "import matplotlib.pyplot as plt\n", + "plt.style.use('seaborn-poster')\n", + "from IPython.display import IFrame, display, Markdown" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import lsst.daf.persistence as dafPersist\n", + "import lsst.daf.base as dafBase\n", + "import lsst.afw.geom as afwGeom\n", + "import lsst.afw.table as afwTable\n", + "from astropy.io import ascii" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Filter some warnings printed by v16.0 of the stack\n", + "# Courtesy of ADW\n", + "warnings.simplefilter(\"ignore\", category=FutureWarning)\n", + "warnings.simplefilter(\"ignore\", category=UserWarning)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Your first table\n", + "\n", + "To begin, we will make a bare-bones afw table so that we can clearly showcase some important concepts. First we will make the simplest possible table, by hand. While creating tables by hand will not likely be the standard use case, it is useful from a tutorial standpoint, as it will allow us to excercise some concepts one at a time" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "slideshow": { + "slide_type": "-" + } + }, + "outputs": [], + "source": [ + "# afw tables need a schema to tell the table how its data are organized\n", + "# Lets have a look at a simple schema:\n", + "min_schema = afwTable.SourceTable.makeMinimalSchema()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# But what is the schema exactly? Printing it out can be informative\n", + "print(min_schema)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Our schema contains 4 Fields: one for each celestial coordinate, an id that uniquely defines it, and a 'parent', which lists the id of the source this source was deblended from. We will deal with the parent column in more detail in a few cells, but for now you can ignore it.\n", + "\n", + "Each field has some accompanying information to go along with it. In addition to its name, we get a helpful docstring describing it. We also get the units that values for this field must have. For example, any value associated with the id key has to be a long integer, and all entries for celestial coordniates have to be instances of an Angle class. We will showcase the Angle class shortly.\n", + "\n", + "If printing out the schema gives you more information that you want, you can get the names. If the names are informative enough, this might be all you need." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "min_schema.getNames()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# We can also add another field to the schema, using a call pattern like this:\n", + "min_schema.addField(\"r_mag\", type=np.float32, doc=\"r band flux\", units=\"mag\")\n", + "# Lets make sure the field was added by printing out the schema once more:\n", + "print(min_schema)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "> We pause here to point out some caveats. \n", + "1. Schemas are append only. You can add new fields, but you cannot remove them. \n", + "2. The units you use have to be understood by astropy. You can find a list of acceptable units at the bottom of this page http://docs.astropy.org/en/stable/units/index.html#module-astropy.units\n", + "3. Specific types are allowed. The short and long of it is you may use floats, ints, longs, strings, Angle objects, and arrays. For more details you can go to the bottom of this page http://doxygen.lsst.codes/stack/doxygen/x_masterDoxyDoc/afw_table.html" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now that we have a schema, we can use it to make a table.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "min_table = afwTable.BaseCatalog(min_schema)\n", + "# our table is empty, and we can check this by looking at its length\n", + "print('our minimal table has {} rows'.format(len(min_table)))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we will add some data to our minimal catalog. Catalogs are collections of 'records', which actually contain the table's data. To add some data to our catalog, we will have to create some records to add on. Records must adhere to the schema that the table has, and so we must add data in field by field." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# make a new record.\n", + "rec = min_table.addNew()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# grab a hold of the keys for the record. We will use these to add data \n", + "keys = min_schema.extract('*') #this returns a dictionary of all the fields\n", + "\n", + "# access the dictionary one field at a time, and grab each field's key. \n", + "# note these are instances of a Key object, and not just simple strings.\n", + "id_key = keys['id'].key\n", + "ra_key = keys['coord_ra'].key\n", + "dec_key = keys['coord_dec'].key\n", + "parent_key = keys['parent'].key\n", + "r_mag_key = keys['r_mag'].key\n", + "\n", + "#use the keys to add data in our record\n", + "rec.set(id_key, 1)\n", + "rec.set(r_mag_key, 19.0)\n", + "rec.set(ra_key, afwGeom.Angle(33.89, units=afwGeom.degrees))\n", + "rec.set(dec_key, afwGeom.Angle(-42.1, units=afwGeom.degrees))\n", + "rec.set(parent_key, 0)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Notice to set the ra and dec, we needed to create `afwGeom.Angle` objects for them. _These are in units of radians by default._ To explicitly tell the object we want it to be in units of degrees, we passed the `afwGeom.dgrees` object in as a keyword argument. You can find more information on `afwGeom.Angle`s [here]( http://doxygen.lsst.codes/stack/doxygen/x_masterDoxyDoc/classlsst_1_1geom_1_1_angle.html)\n", + "\n", + "Additionally, we set the parent to zero. This means this record refers to the object before any deblending occoured. Lets look at our table now to see how it stands." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "min_table" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We will flesh out the parent column a bit more by adding our next record. Notice we can keep using the keys we defined above. Also notice our second record's parent is listed as 1. This means the object 2 was the result of being deblended from object 1, i.e. object 2 is a child object of object 1." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "rec = min_table.addNew()\n", + "rec.set(id_key, 2)\n", + "rec.set(r_mag_key, 18.5)\n", + "rec.set(ra_key, afwGeom.Angle(32.01, units=afwGeom.degrees))\n", + "rec.set(dec_key, afwGeom.Angle(42.5, units=afwGeom.degrees))\n", + "rec.set(parent_key, 1)\n", + "min_table" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "One more caveat to note: in the output in the cell above, the table prints coordinates in radians by default" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# your turn. add one more record to our table\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# now that we have multiple records in our table, we can select particular ones\n", + "min_table[1]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# or iterate over them\n", + "for rec in min_table:\n", + " print(rec.get(id_key))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# you can grab values from particular records by using our keys\n", + "min_table[1].get(id_key)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Using source catalogs produced by DM\n", + "\n", + "A more typical use case will be to read in a catalog that is produced by a DM process. We will show how to read in and work with a source catalog from the Twinkles data in the following section. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Data access\n", + "If you know the path to your source catalog, there is a quick way to read it in. However, it is often more powerful to the 'data butler' to fetch data for you. The butler knows about camera geometry, sensor characteristics, where data are located, and so forth. Having this anciliary information on hand is often very useful. For completeness we will demonstrate both ways of reading in a source catalog, with the note that it is largely considered better practice to use the data butler. \n", + "\n", + "The data butler deserves a tutorial in its own right, and so we will defer further details on it until later. For now, you may think of it as an abstraction that allows you to quickly fetch catalogs for you. The user just needs to point the butler to where to look and what to look for." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# here's the quick and dirty way:\n", + "file_path = '/project/shared/data/Twinkles_subset/output_data_v2/src/v235-fr/R22/S11.fits'\n", + "source_cat = afwTable.SourceCatalog.readFits(file_path)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# here's the way to get the same catalog with a butler:\n", + "\n", + "# first set up our butler by telling it to look at this twinkles directory\n", + "butler = dafPersist.Butler('/project/shared/data/Twinkles_subset/output_data_v2')\n", + "# now we put together a dataId that uniquely specifies a datum\n", + "dataId = {'filter': 'r', 'raft': '2,2', 'sensor': '1,1', 'visit': 235}\n", + "\n", + "# use the dataId and the 'src' to get the source catalog. \n", + "source_cat = butler.get('src', **dataId)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "A few comments are in order on questions you may be having about the butler, and the previous cell. First, there is no good way to know which `dataId`s exist. That means you have to know ahead of time which `dataId`s it makes sense to use. DM is working hard on fixing this. Second, the string `'src'` refers to a very specific data product in the DM philosophy, which is a catalog that contains _the results of different measurement algorithms on detected sources on an individual CCD image_. We will meet some other catalogs later in the tutorial. For now, lets get to know this `src` catalog." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### afw source catalog schemas" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#check its schema. Heads up, the schema is pretty big\n", + "source_cat.getSchema()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "These schemas tend to be pretty large, because every measurement algorithm will create several fields. There are handy ways of grabbing fields that are interesting to you. Suppose you are interested in HSM PSF shape measurements. We can use unix-like pattern matching with the extract method to search the schema. This returns a dictionary where the keys are the schema fields whose names match the pattern you specified, and the values are the fields themselves. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "source_cat.getSchema().extract('*HSM*Psf*')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "If we are just intested in the field names, we can do this:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "source_cat.getSchema().extract('*HSM*Psf*').keys()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "When we dumped the entire schema, the very bottom of the schema contained fields are named 'slot_'. These are called aliases in the schema, and can help you deal with any ambiguity in the table. For example, there are several algorithms used to measure the centroid, and many fileds with 'centroid' in their name as a result. If you want to have quick access to one algorithms measurement result, you can set up a slot alias for it. Lets do a working example on the first record in our table." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "slot_centroid = source_cat[0].getCentroid()\n", + "gauss_cent_x, gauss_cent_y = (source_cat['base_GaussianCentroid_x'][0], source_cat['base_GaussianCentroid_y'][0])\n", + "sdss_cent_x, sdss_cent_y = (source_cat['base_SdssCentroid_x'][0], source_cat['base_SdssCentroid_y'][0])\n", + "\n", + "print('gaussian centroid is {}, {}'.format(gauss_cent_x, gauss_cent_y))\n", + "print('sloan centroid is {}, {}'.format(sdss_cent_x, sdss_cent_y))\n", + "print('slot centroid is {}'.format(slot_centroid))\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# aliasing works with other methods\n", + "psf_flux_key = source_cat.getPsfFluxKey()\n", + "id_key = source_cat.getIdKey()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As advertised, the slot centroid and SDSS centroid are the same. We also used some syntactic sugar to access the `gaussian` centroids and `sdss` centroids, which will be familiar to you if you are an astropy tables user. \n", + "\n", + "> Speaking of astropy tables, you can make an astropy table version of a source catalog:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "source_cat.asAstropy()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we will set aside the schema for this table, and look at the table itself so we can examine its methods.\n", + "\n", + "### afw source catalogs\n", + "\n", + "Source catalogs support a lot of fast operations for common use cases which we will now discuss. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Sorting is supported. by default catalogs are sorted by id\n", + "source_cat.isSorted(id_key)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# You can cut on the catalog.\n", + "# e.g. Make a boolean array to only keep sources with positive psf flux:\n", + "psf_mask = source_cat.getPsfFlux() > 0\n", + "psf_mask &= np.isfinite(source_cat['slot_ApFlux_flux'])\n", + "psf_mask &= np.isfinite(source_cat['slot_ApFlux_fluxSigma'])\n", + "psf_mask &= np.isfinite(source_cat['base_ClassificationExtendedness_value'])\n", + "pos_flux = source_cat.subset(psf_mask)\n", + "\n", + "# You can sort on other keys too:\n", + "flux_key = pos_flux.getPsfFluxKey()\n", + "pos_flux.sort(flux_key)\n", + "pos_flux.isSorted(flux_key)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# You can get the children of particular objects.\n", + "# This is useful if you want to understand how one object was deblended, for example:\n", + "source_cat.getChildren(1010357503918) #the argument is the id of the parent object\n", + "\n", + "# Note that this will only work if the source catalog is sorted on id or parent" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "You can check if catalogs are contiguous in memory." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "source_cat.isContiguous()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "pos_flux.isContiguous()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Some operations are quicker if catalogs are contiguous in memory, like using numpy-like syntax to create masks. You can force the table to be contiguous if you make a deep copy of it. We will show how forcing the table to be contiguous makes the previous operation quicker. Although the speedup is marginal, and not statistically significant, it would be for a much larger catalog. Eli Rykoff performed some benchmark tests showing this is the case for a catalog with about half a million enteries. You can find the full details in a Slack thread [here](https://lsstc.slack.com/archives/C2JPL2DGD/p1525799998000344)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "pos_flux = pos_flux.copy(deep=True)\n", + "pos_flux.isContiguous()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Use the between method to get the indices of values within a range:\n", + "pos_flux.between(1e4,1e6,psf_flux_key)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# The slice object tells you the (start, stop, stride) for values that fit our query.\n", + "# You can check to see that the first record outside the slice is above the flux threshold\n", + "pos_flux[2390].getPsfFlux()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# and that the last element in the slice is inside the threshold\n", + "pos_flux[2389].getPsfFlux()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# your turn. confirm the lower limits of the between query\n", + "# pos_flux[...].getPsfFlux\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#the upper and lower bound methods work similarly\n", + "pos_flux.upper_bound(1e4, psf_flux_key)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "pos_flux.lower_bound(1e6, psf_flux_key)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Example: Star-Galaxy Separation\n", + "\n", + "Now that we have introduced the functionality of the source catalog and its schema, we will do a toy example of star-galaxy separation. This small demo will also flags and fields that users are use, and ultimately make a plot." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# let's select sources that were not deblended\n", + "select_mask = source_cat['deblend_nChild'] == 0\n", + "select_mask &= source_cat['parent'] == 0\n", + "\n", + "# use the extendedness column for a mock star/galaxy seperation\n", + "# we only want to use columns where this algorithm worked\n", + "# the flag is set true if there was a failure, so we invert the flag values here\n", + "select_mask &= ~source_cat['base_ClassificationExtendedness_flag']\n", + "\n", + "# we will also use the sloan shapes to measure size\n", + "select_mask &= ~source_cat['base_SdssShape_flag']\n", + "\n", + "# and a simple aperture flux for the brightness\n", + "select_mask &= ~source_cat['base_CircularApertureFlux_12_0_flag']\n", + "\n", + "# only consider sources with positive flux\n", + "select_mask &= source_cat['base_CircularApertureFlux_12_0_flux'] > 0\n", + "\n", + "size_bright_cat = source_cat.subset(select_mask)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we make a crude size magnitude diagram, color coding the data by their 'extendedness value'. The extendedness will be 1 for extended sources-like galaxies-and 0 for point sources-like stars. One hopes the stars will all live on the stellar locus..." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.scatter(np.log(size_bright_cat['base_CircularApertureFlux_12_0_flux']), \n", + " size_bright_cat['base_SdssShape_xx'] + size_bright_cat['base_SdssShape_yy'],\n", + " c=size_bright_cat['base_ClassificationExtendedness_value'],\n", + " s=4)\n", + "plt.xlabel('log flux')\n", + "plt.ylabel('size px^2')\n", + "plt.ylim([0,60]) #zoom in to make the stellar locus clearer" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Our plot shows some star galaxy separation, but also has other interesting features. Some detected sources appear to be smaller than the PSF, some of the point sources have a (crudely) calculated size that occupy the same parameter space as extended sources, and there are a few extremely faint detected point sources. We will leave it to you to delve into this mystery further as a homework assignment since we are primarily focused on understanding tables in this tutorial. By making this plot we exercised some of the methods of the catalog and its schema, to do a minimal analysis example." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Operations with multiple tables/catalogs\n", + "\n", + "In the next section we will show operations which involve two or more catalogs.\n", + "\n", + "#### Table concatenation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Grab a second catalog using the butler:\n", + "dataId2 = {'filter': 'r', 'raft': '2,2', 'sensor': '1,1', 'visit': 236}\n", + "source_cat2 = butler.get('src', dataId2)\n", + "\n", + "# Put our catalogs in a list:\n", + "catalogList = [source_cat, source_cat2]\n", + "\n", + "# The following concatenation function is courtesy of Jim Bosch:\n", + "def concatenate(catalogList):\n", + " from functools import reduce\n", + " \"\"\"Concatenate multiple catalogs (FITS tables from lsst.afw.table)\"\"\"\n", + " schema = catalogList[0].schema\n", + " for i, c in enumerate(catalogList[1:]):\n", + " #check that the schema in the tables are the same\n", + " #we can only cat them if this is true\n", + " if c.schema != schema:\n", + " raise RuntimeError(\"Schema for catalog %d not consistent\" % (i+1))\n", + "\n", + " # prepare the master catalog\n", + " out = afwTable.BaseCatalog(schema)\n", + " num = reduce(lambda n, c: n + len(c), catalogList, 0)\n", + " # set aside enough space for all the records and their pointers\n", + " out.reserve(num)\n", + "\n", + " # stick in all the records from all catalogs into the master catalog\n", + " for catalog in catalogList:\n", + " for record in catalog:\n", + " out.append(out.table.copyRecord(record))\n", + "\n", + " return out\n", + "\n", + "cat_source = concatenate(catalogList)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Catalog matching\n", + "\n", + "Quick positional matching is supported by the stack, and offers some useful functionality. In the next example, we will match one of the Twinkles truth catalogs against our stack-produced detection catalog. To be as DM like as possible, we will read in the truth catalog, make an `afwTable.BaseTable` version of the truth catalog to match against the src catalog produced by DM. As we do this, you'll see us re-use the code we showed above.\n", + "\n", + "First order of business is to grab and reformat the truth table." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "truth_table = ascii.read('/project/shared/data/Twinkles_subset/truth/sprinkled_lens_230_J2000.txt')\n", + "truth_table" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "min_schema = afwTable.SourceTable.makeMinimalSchema()\n", + "min_schema.addField(\"r_mag\", type=np.float32, doc=\"r band mag\", units=\"mag\")\n", + "truth_dm = afwTable.SourceCatalog(min_schema)\n", + "\n", + "# grab the keys to the schema\n", + "# grab a hold of the keys for the record. We will use these to add data \n", + "keys = min_schema.extract('*') #this returns a dictionary of all the fields\n", + "\n", + "# access the dictionary one field at a time, and grab each field's key\n", + "id_key = keys['id'].key\n", + "ra_key = keys['coord_ra'].key\n", + "dec_key = keys['coord_dec'].key\n", + "parent_key = keys['parent'].key\n", + "r_mag_key = keys['r_mag'].key\n", + "\n", + "# loop over the astropy table version of the truth catalog\n", + "# making new records in our DM catalog as we go along\n", + "for record in truth_table:\n", + " newRec = truth_dm.addNew()\n", + " newRec.set(id_key, record['galid'])\n", + " newRec.set(r_mag_key, record['mag'])\n", + " newRec.set(ra_key, afwGeom.Angle(record['ra'], afwGeom.degrees))\n", + " newRec.set(dec_key, afwGeom.Angle(record['dec'], afwGeom.degrees))\n", + "# newRec.set(ra_key, afwGeom.Angle(record['ra']*(np.pi/180.)))\n", + "# newRec.set(dec_key, afwGeom.Angle(record['dec']*(np.pi/180.)))\n", + "\n", + "# this makes the table contiguous \n", + "truth_dm = truth_dm.copy(deep=True)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we grab a src catalog that overlaps with the truth catalog. The README.txt file in the directory with the truth catalog tells us we need to use visit 230:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dataId = {'filter': 'r', 'raft': '2,2', 'sensor': '1,1', 'visit': 230}\n", + "src_twinkles = butler.get('src',dataId)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We will want to compare the magnitude of matched sources in both catalogs. As we saw from examining the `src` catalog schemas above, flux measurements from different photometry alogrithms are avaliable for every source, but the magnitudes are not explicitly given. We can get calibrated magnitudes from flux measurements by using the `calexp_calib` data product." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "calexp_calib_twinkles = butler.get('calexp_calib', dataId)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here's how to use the `calexp_calib` object to return magnitudes along with errors, given flux and flux errors:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "src_tw_mags, src_tw_magErr = calexp_calib_twinkles.getMagnitude(src_twinkles['base_GaussianFlux_flux'],\n", + " src_twinkles['base_GaussianFlux_fluxSigma'])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# remove sources with bag magnitudes and prune catalog\n", + "truth_mag_max = truth_dm['r_mag'].max()\n", + "truth_mag_min = truth_dm['r_mag'].min()\n", + "mask = np.isfinite(src_tw_mags)\n", + "mask &= (src_tw_mags < truth_mag_max)\n", + "mask &= (src_tw_mags > truth_mag_min)\n", + "src_twinkles = src_twinkles.subset(mask)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Before we do the matching, let's get a sense of the overlap of these catalogs by plotting them both in RA-DEC space:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.scatter(src_twinkles['coord_ra'], src_twinkles['coord_dec'], s=5, label='DM')\n", + "plt.scatter(truth_dm['coord_ra'],truth_dm['coord_dec'], s=20, label='Truth')\n", + "plt.legend();" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We may expect to find about 30 matches or so. In order to match catalogs, we must provide a `MatchControl` instance. The `MatchControl` provides configurations for catalog matching. It has three 'switches' in the form of class attributes. They are defined as follows:\n", + "\n", + "1. `findOnlyClosest`: True by default. If False, all other sources within a search radius are also matched \n", + "2. `includeMismatches`: False by default. If False, sources with no match are not reported in the match catalog. If True, sources with no match are included in the match catalog with Null as their match\n", + "3. `symmetricMatch`: False by default. If False, the match between source a from catalog a with source b from catalog b is reported alone. If True, the symmetric match between source b and a is also reported." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# get a match control, we will keep the default configuration\n", + "mc = afwTable.MatchControl()\n", + "\n", + "# match our two catalogs, setting the match threshold to be one arcsecond\n", + "matches = afwTable.matchRaDec(truth_dm, src_twinkles, afwGeom.Angle(1,afwGeom.arcseconds), mc)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "`afwTable.matchRaDec` returns a list, where each element is an instance of a `Match` class. The `Match` class has three attributes, which gives us information about the matched sources. Let us unpack this a bit before moving on to some analysis" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# how many sources were actually matched?\n", + "len(matches)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# lets examine the first element in the matches list\n", + "# we can grab the record corresponding to this source in the truth_dm catalog\n", + "# using the first attribute\n", + "matches[0].first" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# likewise the second attribute gives us the record from the src_twinkles catalog\n", + "\n", + "matches[0].second" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#finally the angular seperation is given in radians in the distance attribute\n", + "matches[0].distance" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now lets put this all together. We make a plot showing the angular separation between matched sources as a function of the truth catalog's magnitude:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "mag_truth = [m.first['r_mag'] for m in matches]\n", + "dist_resid = [m.distance for m in matches]\n", + "dist_resid = np.array(dist_resid)* 360*60*60/ (2*np.pi) #convert rad to arcseconds\n", + "\n", + "plt.scatter(mag_truth, dist_resid)\n", + "plt.ylabel('angular sep arcseconds');\n", + "plt.xlabel('r truth');" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In the previous example we only kept the nearest neighbor matches. Now we will show how you can collect *all* matches within the search radius, by overwriting the `findOnlyClosest` attribute of the match control." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#get a match control, we will keep the default configuration\n", + "mc = afwTable.MatchControl()\n", + "mc.findOnlyClosest = False\n", + "#match our two catalogs\n", + "matches_all = afwTable.matchRaDec(truth_dm, src_twinkles, afwGeom.Angle(1,afwGeom.arcseconds), mc)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's see if we get a few more matches! " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "len(matches_all)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To try to understand where these extra matches come from, we print out the id's from the truth catalog and their matched sources in the DM catalog, along with their angular separation and magnitude residual. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "for m in matches_all:\n", + " obj_id_1 = m.first.getId()\n", + " obj_id_2 = m.second.getId()\n", + " dist = m.distance* 360*60*60/ (2*np.pi)\n", + " mag_resid = m.first['r_mag'] - calexp_calib_twinkles.getMagnitude(m.second.getModelFlux())\n", + " print('id 1: {} id 2: {} distance {} mag {}'.format(obj_id_1, obj_id_2, dist, mag_resid))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can see that some of the matches have very similar magnitude residuals. Look at the 6th and 7th to last row, for example:\n", + "```\n", + "id 1: 279852058 id 2: 988882667634 distance 0.1571098610746062 mag -0.3084721056101962\n", + "\n", + "id 1: 279852058 id 2: 988882665917 distance 0.15712272456292592 mag -0.30847215544412876\n", + "```" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As it happens, the extra matches are due to deblending." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "src_twinkles.getChildren(988882665917)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "If you are only interested in the ids and the angular separation, you can pack the matches into a table." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "matches_table = afwTable.packMatches(matches)\n", + "matches_table" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "You can unpack the matches too:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "unpack_matches = afwTable.unpackMatches(matches_table, truth_dm, src_twinkles)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Hopefully this gives you some idea of the power of `afwTable`s in matching catalogs together, and understanding the deblending that has been carried out." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Summary\n", + "\n", + "In this tutorial we introduced afw tables. We introduced schemas, how to navigate them, add to them, and how to create tables based off of them. We covered data access to catalogs produced by DM using the data butler. We went over a some typical use cases for source catalogs, like catalog matching, and understanding bread and butter measurement algorithms." + ] + } + ], + "metadata": { + "celltoolbar": "Slideshow", + "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" + }, + "livereveal": { + "scroll": true, + "start_slideshow_at": "selected" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}