diff --git a/Basics/Calexp_guided_tour.ipynb b/Basics/Calexp_guided_tour.ipynb index dc7f37d5..98a95e41 100644 --- a/Basics/Calexp_guided_tour.ipynb +++ b/Basics/Calexp_guided_tour.ipynb @@ -460,7 +460,7 @@ "bbox = afwGeom.Box2I()\n", "bbox.include(afwGeom.Point2I(2200,3200))\n", "bbox.include(afwGeom.Point2I(2800,3800))\n", - "cutout = calexp.Factory(calexp, bbox, afwImage.LOCAL)" + "cutout = calexp[bbox]" ] }, { @@ -510,7 +510,7 @@ "metadata": {}, "outputs": [], "source": [ - "cutout_calexp = butler.get('calexp_sub', bbox=bbox, immediate=True, dataId=dataId)\n", + "cutout_calexp = butler.get('calexp_sub', bbox=bbox, dataId=dataId)\n", "cutout_calexp.getDimensions()" ] }, @@ -654,7 +654,16 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "The result of the `set` command above shows that a calexp and a coadd have the same methods." + "The result of the `set` command above shows that a calexp and a coadd have the same methods. This is expected, because they are the same class." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(calexp.__class__, coadd.__class__)" ] }, { diff --git a/Deblending/README.rst b/Deblending/README.rst new file mode 100644 index 00000000..6fd8286c --- /dev/null +++ b/Deblending/README.rst @@ -0,0 +1,36 @@ +Deblending +========== + +This folder contains a set of tutorial notebooks exploring the deblending of LSST objects. See the index table below for links to the notebook code, and an auto-rendered view of the notebook with outputs. + + +.. list-table:: + :widths: 10 20 10 10 + :header-rows: 1 + + * - Notebook + - Short description + - Links + - Owner + + + * - **SCARLET Tutorial** + - Introduction to the SCARLET deblender, how to configure and run it. + - `ipynb `_, + `rendered `_ + + .. image:: https://github.com/LSSTScienceCollaborations/StackClub/blob/rendered/Deblending/log/scarlet_tutorial.svg + :target: https://github.com/LSSTScienceCollaborations/StackClub/blob/rendered/Deblending/log/scarlet_tutorial.log + + - `Fred Moolekamp `_ + + + * - **Deblending in DRP** + - Where and how the deblending happens, in the DRP pipeline. + - `ipynb `_, + `rendered `_ + + .. image:: https://github.com/LSSTScienceCollaborations/StackClub/blob/rendered/Deblending/log/lsst_stack_deblender.svg + :target: https://github.com/LSSTScienceCollaborations/StackClub/blob/rendered/Deblending/log/lsst_stack_deblender.log + + - `Fred Moolekamp `_ diff --git a/Deblending/lsst_stack_deblender.ipynb b/Deblending/lsst_stack_deblender.ipynb new file mode 100755 index 00000000..bcfa264b --- /dev/null +++ b/Deblending/lsst_stack_deblender.ipynb @@ -0,0 +1,635 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Using the LSST Stack Multiband Deblender \n", + "
Owner(s): **Fred Moolekamp** ([@fred3m](https://github.com/LSSTScienceCollaborations/StackClub/issues/new?body=@fred3m))\n", + "
Last Verified to Run: **2018-08-17**\n", + "
Verified Stack Release: **w_2018_32**\n", + "\n", + "This tutorial is designed to illustrate how to execute the multiband deblender (*scarlet*) in the LSST stack. This includes a brief introduction to LSST stack objects including:\n", + "\n", + " - Geometry classes from `lsst.geom`, such as points and boxes.\n", + " - Higher-level astronomical primitives from `lsst.afw`, such as the `Image`, `Exposure`, and `Psf` classes.\n", + " - Our core algorithmic `Task` classes, including those for source detection, deblending, and measurement.\n", + " \n", + "We'll be working with coadded images made from Subaru Hyper Suprime-Cam (HSC) data in the COSMOS field. We've taken a recent LSST reprocessing of the HSC-SSP UltraDeep COSMOS field (see [this page](https://confluence.lsstcorp.org/display/DM/S18+HSC+PDR1+reprocessing) for information on that reprocessing, and [this page](https://hsc-release.mtk.nao.ac.jp/doc/) for the data), and added simulated stars from a scaled [SDSS catalog](http://www.sdss.org/dr14/data_access/value-added-catalogs/?vac_id=photometry-of-crowded-fields-in-sdss-for-galactic-globular-and-open-clusters). The result is a very deep image (deeper than the 10-year LSST Deep-Wide-Fast survey, though not as deep as LSST Deep Drilling fields will be) with both a large number of galaxies and region full of stars.\n", + "\n", + "This tutorial is based on Jim Bosch's globular cluster tutorial, however in it's present state *scarlet* is unable to process the crowded field (most likely) due to poor initial conditions for the sources in the field. So instead we use a region of the image outside of the cluster.\n", + "\n", + "### Learning Objectives:\n", + "\n", + "After working through this tutorial you should be able to: \n", + "1. Configure and run the LSST multiband deblender on a test list of objects;\n", + "2. Understand its task context in the DRP pipeline.\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": [ + "## Imports\n", + "\n", + "We'll start with some standard imports of both LSST and third-party packages." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "from lsst.daf.persistence import Butler\n", + "from lsst.geom import Box2I, Box2D, Point2I, Point2D, Extent2I, Extent2D\n", + "from lsst.afw.image import Exposure, Image, PARENT, MultibandExposure, MultibandImage\n", + "from lsst.afw.detection import MultibandFootprint" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Reading Data\n", + "\n", + "We'll be retrieving data using the `Butler` tool, which manages where various datasets are stored on the filesystem (and can in principle manage datasets that aren't even stored as files, though all of these are).\n", + "\n", + "We start by creating a `Butler` instance, pointing it at a *Data Repository* (which here is just a root directory)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "butler = Butler(\"/project/jbosch/tutorials/lsst2018/data\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Datasets managed by a butler are identified by a dictionary *Data ID* (specifying things like the visit number or sky patch) and a string *DatasetType* (such as a particular image or catalog). Different DatasetTypes have different keys, while different instances of the same Dataset Type have different values. All of the datasets we use in this tutorial will correspond to the same patch of sky, so they'll have at least the keys in the dictionary in the next cell (they will also have `filter`, but with different values):" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dataId = {\"tract\": 9813, \"patch\": \"4,4\"}" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can now use those to load a set of *grizy* coadds, which we'll put directly in a dictionary. The result of each `Butler.get` call is in this case an `lsst.afw.image.Exposure` object, an image that actually contains three \"planes\" (the main image, a bit mask, and a variance image) as well as many other objects that describe the image, such as its PSF and WCS. Note that we (confusingly) use `Exposures` to hold coadd images as well as true single-exposure images, but combine them into a `MultibandExposure`, which contains an exposure in each band.\n", + "\n", + "The DatasetType here is `deepCoadd_calexp` (a coadd on which we've already done some additional processing, such as subtracting the background and setting some mask values), and the extra `filter` argument gets appended to the Data ID." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "filters = \"grizy\"\n", + "coadds = [butler.get(\"deepCoadd_calexp\", dataId, filter=\"HSC-{}\".format(f.upper())) for f in filters]\n", + "coadds = MultibandExposure.fromExposures(filters, coadds)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Making and displaying color composite images\n", + "\n", + "We'll start by just looking at the images, as 3-color composites. We'll use astropy to build those as a nice way to demonstrate how to get NumPy arrays from the `MultibandImage` objects (the images in `coadds`). (LSST also has code to make 3-color composites using the same algorithm, and in fact the Astropy implementation is based on ours, but now that it's in Astropy we'll probably retire ours.)\n", + "\n", + "We'll just use matplotlib to display the images themselves. We'll use Firefly for other image display tasks later, but while Firefly itself supports color-composites, work on our preferred composition algorithm is still in progress, and we haven't quite finished connecting that functionality to the Python client we'll demonstrate here." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from astropy.visualization import make_lupton_rgb\n", + "import matplotlib.pyplot as plt\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We'll use the following function a few times to display color images. It's worth reading through the implementation carefully to see what's going on." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def showRGB(image, bgr=\"gri\", ax=None, fp=None, figsize=(8,8), stretch=1, Q=10):\n", + " \"\"\"Display an RGB color composite image with matplotlib.\n", + " \n", + " Parameters\n", + " ----------\n", + " image : `MultibandImage`\n", + " `MultibandImage` to display.\n", + " bgr : sequence\n", + " A 3-element sequence of filter names (i.e. keys of the exps dict) indicating what band\n", + " to use for each channel. If `image` only has three filters then this parameter is ignored\n", + " and the filters in the image are used.\n", + " ax : `matplotlib.axes.Axes`\n", + " Axis in a `matplotlib.Figure` to display the image.\n", + " If `axis` is `None` then a new figure is created.\n", + " fp: `lsst.afw.detection.Footprint`\n", + " Footprint that contains the peak catalog for peaks in the image.\n", + " If `fp` is `None` then no peak positions are plotted.\n", + " figsize: tuple\n", + " Size of the `matplotlib.Figure` created.\n", + " If `ax` is not `None` then this parameter is ignored.\n", + " stretch: int\n", + " The linear stretch of the image.\n", + " Q: int\n", + " The Asinh softening parameter.\n", + " \"\"\"\n", + " # If the image only has 3 bands, reverse the order of the bands to produce the RGB image\n", + " if len(image) == 3:\n", + " bgr = image.filters\n", + " # Extract the primary image component of each Exposure with the .image property, and use .array to get a NumPy array view.\n", + " rgb = make_lupton_rgb(image_r=image[bgr[2]].array, # numpy array for the r channel\n", + " image_g=image[bgr[1]].array, # numpy array for the g channel\n", + " image_b=image[bgr[0]].array, # numpy array for the b channel\n", + " stretch=stretch, Q=Q) # parameters used to stretch and scale the pixel values\n", + " if ax is None:\n", + " fig = plt.figure(figsize=figsize)\n", + " ax = fig.add_subplot(1,1,1)\n", + " \n", + " # Exposure.getBBox() returns a Box2I, a box with integer pixel coordinates that correspond to the centers of pixels.\n", + " # Matplotlib's `extent` argument expects to receive the coordinates of the edges of pixels, which is what\n", + " # this Box2D (a box with floating-point coordinates) represents.\n", + " integerPixelBBox = image[bgr[0]].getBBox()\n", + " bbox = Box2D(integerPixelBBox)\n", + " ax.imshow(rgb, interpolation='nearest', origin='lower', extent=(bbox.getMinX(), bbox.getMaxX(), bbox.getMinY(), bbox.getMaxY()))\n", + " if fp is not None:\n", + " for peak in fp.getPeaks():\n", + " ax.plot(peak.getIx(), peak.getIy(), \"bx\", mew=2)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Notice that we can slice `MultibandImage` objects (as well a `MultibandExposure` objects) along the filter dimension using the filter names as indices. Like `Exposure` objects, `MultibandExposure` objects have `image`, `mask`, and `variance` properties that contain the image, mask plane, and variance of the `Exposure` respectively. For now we will only worry about the `image` property, although internal deblending and measurement algorithms make use of all three objects (when available)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "showRGB(coadds[:\"z\"].image, figsize=(10, 10))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "showRGB(coadds[\"i\":].image, figsize=(10, 10))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Those images are a full \"patch\", which is our usual unit of processing for coadds - it's about the same size as a single LSST sensor (exactly the same in pixels, smaller in terms of area because these use HSC's smaller pixel scale). That's a bit unweildy (just because waiting for processing to happen isn't fun in a tutorial setting), so we'll reload our dict with sub-images centered on the region of interest.\n", + "\n", + "Note that we can load the sub-images directly with the `butler`, by appending `_sub` to the DatasetType and passing a `bbox` argument. If you want to see the region of the image with the cluster, use `clusterBBox` below, however as mentioned above, that region is too memory intensive for the current version of *scarlet*. Instead use `sampleBBox` to select a sub-region of the image (note that we add a small frame around each blend to include more background regions, which are important for the detection algorithm)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "frame = 50\n", + "clusterBBox = Box2I(Point2I(18325, 17725), Extent2I(400, 350))\n", + "\n", + "#sampleBBox = Box2I(Point2I(18699-frame, 17138-frame), Extent2I(93+2*frame, 104+2*frame))\n", + "#sampleBBox = Box2I(Point2I(16424-frame, 17806-frame), Extent2I(55+2*frame, 62+2*frame))\n", + "#sampleBBox = Box2I(Point2I(17838-frame, 18945-frame), Extent2I(111+2*frame, 103+2*frame))\n", + "sampleBBox = Box2I(Point2I(19141-frame, 18228-frame), Extent2I(63+2*frame, 87+2*frame))\n", + "\n", + "subset = coadds[:, sampleBBox]\n", + "# Due to a bug in the code the PSF isn't copied properly.\n", + "# The code below copies the PSF into the `MultibandExposure`,\n", + "# but will be unecessary in the future\n", + "for f in subset.filters:\n", + " subset[f].setPsf(coadds[f].getPsf())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "showRGB(subset.image)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Basic Processing\n", + "\n", + "Now we'll try the regular LSST processing tasks, with a simpler configuration than we usually use to process coadds, just to avoid being distracted by complexity. This includes\n", + "\n", + " - Detection (`SourceDetectionTask`): given an `Exposure`, find above-threshold regions and peaks within them (`Footprints`), and create a *parent* source for each `Footprint`.\n", + " - Deblending (`MultibandDeblendTask`): given a `MultibandExposure` and a catalog of parent sources, create a *child* source for each peak in every `Footprint` that contains more than one peak. Each child source is given a `HeavyFootprint`, which contains both the pixel region that source covers and the fractional pixel values associated with that source. A `SourceDeblendTask` is also available using the single band SDSS-HSC deblender that takes a single band `Exposure`).\n", + " - Measurment (`SingleFrameMeasurementTask`): given an `Exposure` and a catalog of sources, run a set of \"measurement plugins\" on each source, using deblended pixel values if it is a child. Notice that measurement is still performed on single band catalogs, since none of the measurement algorithms work for multiband data.\n", + "\n", + "We'll start by importing these, along with the `SourceCatalog` class we'll use to hold the outputs." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from lsst.meas.algorithms import SourceDetectionTask\n", + "from lsst.meas.deblender import MultibandDeblendTask\n", + "from lsst.meas.base import SingleFrameMeasurementTask\n", + "from lsst.afw.table import SourceCatalog" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We'll now construct all of these `Tasks` before actually running any of them. That's because `SourceDeblendTask` and `SingleFrameMeasurementTask` are constructed with a `Schema` object that records what fields they'll produce, and they modify that schema when they're constructed by adding columns to it. When we run the tasks later, they'll need to be given a catalog that includes all of those columns, **but we can't add columns to a catalog that already exists**.\n", + "\n", + "To recap, the sequence looks like this:\n", + "\n", + " 1. Make a (mostly) empty schema.\n", + " 2. Construct all of the `Task`s (in the order you plan to run them), which adds columns to the schema.\n", + " 3. Make a `SourceCatalog` object from the *complete* schema.\n", + " 4. Pass the same `SourceCatalog` object to each `Task` when you run it." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "schema = SourceCatalog.Table.makeMinimalSchema()\n", + "\n", + "detectionTask = SourceDetectionTask(schema=schema)\n", + "\n", + "config = MultibandDeblendTask.ConfigClass()\n", + "config.usePsfConvolution = True\n", + "config.conserveFlux = True\n", + "config.maxIter = 100\n", + "deblendTask = MultibandDeblendTask(schema=schema, config=config)\n", + "\n", + "# We'll customize the configuration of measurement to just run a few plugins.\n", + "# The default list of plugins is much longer (and hence slower).\n", + "measureConfig = SingleFrameMeasurementTask.ConfigClass()\n", + "measureConfig.plugins.names = [\"base_SdssCentroid\", \"base_PsfFlux\", \"base_SkyCoord\"]\n", + "# \"Slots\" are aliases that provide easy access to certain plugins.\n", + "# Because we're not running the plugin these slots refer to by default,\n", + "# we need to disable them in the configuration.\n", + "measureConfig.slots.apFlux = None\n", + "measureConfig.slots.instFlux = None\n", + "measureConfig.slots.shape = None\n", + "measureConfig.slots.modelFlux = None\n", + "measureConfig.slots.calibFlux = None\n", + "measureTask = SingleFrameMeasurementTask(config=measureConfig, schema=schema)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The first step we'll run is detection, which actually returns a new `SourceCatalog` object rather than working on an existing one.\n", + "\n", + "Instead, it takes a `Table` object, which is sort of like a factory for records. We won't use it directly after this, and it isn't actually necessary to make a new `Table` every time you run `MultibandDetectionTask` (but you can only create one after you're done adding columns to the schema).\n", + "\n", + "`Task`s that return anything do so via a `lsst.pipe.base.Struct` object, which is just a simple collection of named attributes. The only return values we're interested is `sources`. That's our new `SourceCatalog`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "table = SourceCatalog.Table.make(schema)\n", + "detectionResult = detectionTask.run(table, subset[\"r\"])\n", + "catalog = detectionResult.sources" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's take a quick look at what's in that catalog. First off, we can look at its schema:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "catalog.schema" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note that this includes a lot of columns that were actually added by the deblend or measurement steps; those will all still be blank (`0` for integers or flags, `NaN` for floating-point columns).\n", + "\n", + "In fact, the only columns filled by `SourceDetectionTask` are the IDs. But it also attaches `Footprint` objects, which don't appear in the schema. You can retrieve the `Footprint` by calling `getFootprint()` on a row:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "footprint = catalog[0].getFootprint()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "`Footprints` have two components:\n", + " - a `SpanSet`, which represents an irregular region on an image via a list of (y, x0, x1) `Spans`;\n", + " - a `PeakCatalog`, a slightly different kind of catalog whose rows represent peaks within that `Footprint`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(footprint.getSpans())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(footprint.getPeaks())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "If we actually look at the footprints in the catalog we see that some have only a single peak, while others have multiple peaks that need to be deblended.\n", + "\n", + "To display only the pixels contained in the footprint (and not other pixels in the bounding box) we create a `MultibandFootprint`, which is a `HeavyFootprint` that contains a `SpanSet`, `PeakCatalog`, and `flux` values for all of the pixels in the `SpanSet`. In this case the `flux` is the total measured flux in the image, since no deblending has taken place yet." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "for src in catalog:\n", + " fp = src.getFootprint()\n", + " img = coadds[:,fp.getBBox()].image\n", + " mfp = MultibandFootprint.fromImages(coadds.filters, image=img, footprint=fp)\n", + " showRGB(mfp.getImage().image, fp=fp, figsize=(3,3))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "It's worth noting that while the peaks *can* have both an integer-valued position and a floating-point position, they're the same right now; `SourceDetectionTask` currently just finds the pixels that are local minima and doesn't try to find their sub-pixel locations. That's left to the centroider, which is part of the measurement stage.\n", + "\n", + "Before we can get to that point, we need to run the deblender:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fluxCatalog, templateCatalog = deblendTask.run(coadds, catalog)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "`MultibandDeblendTask` always returns two catalogs, a `templateCatalog` that contains the model outputs from *scarlet* and a `fluxCatalog`, which uses the *scarlet* models as weights to redistribute the flux from the input image (in other words they contain flux-conserved models). If `MultibandDeblendTask.config.saveTemplates` is `False`, then `templateCatalog` will be `None`. Similarly, if `MultibandDeblendTask.config.conserveFlux` is `False` then the `fluxCatalog` will be `None` (and the code will run slightly faster, since it doesn't have to reweight the flux, however this is a small fraction of the processing time).\n", + "\n", + "The deblender itself sets the `parent` column for each source, which is `0` for objects with no parent, and all of the columns that begin with `deblend_` and also adds new rows to the catalog for each child. It does *not* remove the parent rows it created those child rows from, and this is intentional, because we want to measure both \"interpretations\" of the blend family: one in which there is only one object (the parent version) and one in which there are several (the children). Before doing any science with the outputs of an LSST catalog, it's important to remove one of those interpretations (typically the parent one). That can be done by looking at the `deblend_nChild` and `parent` fields:\n", + "\n", + " - `parent` is the ID of the source from which this was deblended, or `0` if the source is itself a parent.\n", + " - `deblend_nChild` is the number of child sources this source has (so it's `0` for sources that are themselves children or were never blended).\n", + " \n", + "Together, these define two particularly useful filters:\n", + "\n", + " - `deblend_nChild == 0`: never-blended object or de-blended child\n", + " - `deblend_nChild == 0 and parent == 0`: never-blended object\n", + " \n", + "The first is what you'll usually want to use; the second is what to use if you're willing to throw away some objects (possibly many) because you don't trust the deblender.\n", + "\n", + "The last processing step for our purposes is running measurement, which must be done on each catalog, in each band (if we want measurements for all of them):" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "measureTask.run(templateCatalog[\"r\"], coadds['r'])\n", + "measureTask.run(templateCatalog[\"i\"], coadds['i'])\n", + "measureTask.run(fluxCatalog[\"r\"], coadds['r'])\n", + "measureTask.run(fluxCatalog[\"i\"], coadds['i'])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Due to an unfortunate bug in the deblender task, the resulting catalogs are not contiguous and we need to copy them into new objects to use them appropriately. This step can be avoided in the near future." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import lsst.afw.table as afwTable\n", + "\n", + "for f in filters:\n", + " _catalog = afwTable.SourceCatalog(templateCatalog[f].table.clone())\n", + " _catalog.extend(templateCatalog[f], deep=True)\n", + " templateCatalog[f] = _catalog\n", + " _catalog = afwTable.SourceCatalog(fluxCatalog[f].table.clone())\n", + " _catalog.extend(fluxCatalog[f], deep=True)\n", + " fluxCatalog[f] = _catalog" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Since we care about deblending (for the sake of this tutorial) lets look at the results from the 13th blend displayed above. We use the `HeavyFootprint`s from the catalog sources that have the same parent (parent 13 from above) to build a model for the entre scene, and to compare the results of the flux conserved and *scarlet* models. In the process we look at both the *scarlet* and flux conserved models." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "from lsst.afw.detection import MultibandFootprint\n", + "from lsst.afw.image import MultibandImage\n", + "\n", + "# Use the 13th parent in the blend\n", + "# Note: this is not the parent ID, but the 13th source in the catalog\n", + "parentIdx = 13\n", + "\n", + "# Create empty multiband images to model the entire scene\n", + "templateModel = MultibandImage.fromImages(coadds.filters,\n", + " [Image(fluxCatalog[\"r\"][parentIdx].getFootprint().getBBox(), dtype=np.float32)\n", + " for b in range(len(filters))])\n", + "fluxModel = MultibandImage.fromImages(coadds.filters,\n", + " [Image(fluxCatalog[\"r\"][parentIdx].getFootprint().getBBox(), dtype=np.float32)\n", + " for b in range(len(filters))])\n", + "\n", + "# Only use the subset catalogs with the same parent\n", + "parentId = fluxCatalog[\"r\"][parentIdx].get(\"id\")\n", + "fluxChildren = {b: fluxCatalog[b][fluxCatalog[b].get(\"parent\")==parentId] for b in filters}\n", + "templateChildren = {b: templateCatalog[b][templateCatalog[b].get(\"parent\")==parentId] for b in filters}\n", + "assert(len(fluxChildren)==len(templateChildren))\n", + "\n", + "for n in range(len(templateChildren[\"r\"])):\n", + " # Add the source model to the model of the entire scene\n", + " fp = MultibandFootprint(coadds.filters, [templateChildren[b][n].getFootprint() for b in filters])\n", + " _fp = MultibandFootprint(coadds.filters, [fluxChildren[b][n].getFootprint() for b in filters])\n", + " templateModel[:, fp.getBBox()].array += fp.getImage(fill=0).image.array\n", + " fluxModel[:, _fp.getBBox()].array += _fp.getImage(fill=0).image.array\n", + "\n", + " # Show the model\n", + " fig = plt.figure(figsize=(6, 3))\n", + " ax = [fig.add_subplot(1, 2, n+1) for n in range(2)]\n", + " ax[0].set_title(\"scarlet\")\n", + " ax[1].set_title(\"flux conserved\")\n", + " showRGB(fp.getImage().image, ax=ax[0])\n", + " showRGB(_fp.getImage().image, ax=ax[1])\n", + " plt.show()\n", + "\n", + "templateResidual = MultibandImage(coadds.filters,\n", + " coadds[:, templateModel.getBBox()].image.array - templateModel.array)\n", + "fluxResidual = MultibandImage(coadds.filters,\n", + " coadds[:, fluxModel.getBBox()].image.array - fluxModel.array)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Finally we look at the full models and the residuals. As expected, there are no residuals for the flux conserved model since all of the flux in the image (that is within the footprint) is added to one of the sources. In this particular case that works fine, but in instances where one or more sources were not detected this can cause one source to have its flux contaminated with its neighbor." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "for title, model, residual in [[\"scarlet\", templateModel, templateResidual], [\"flux conserved\", fluxModel, fluxResidual]]:\n", + " fig = plt.figure(figsize=(15,8))\n", + " ax = [fig.add_subplot(1, 2, n+1) for n in range(2)]\n", + " ax[0].set_title(\"{0} model\".format(title))\n", + " ax[1].set_title(\"{0} residual\".format(title))\n", + " showRGB(model, ax=ax[0])\n", + " showRGB(residual ,ax=ax[1], Q=0)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Exercises\n", + "\n", + "- Use some of the other`sampleBBox` regions and run through the code again, from source detection through measurment and blending displays. Don't foget to change the parent index to view only the children of the correct blend.\n", + "- Play around with other *scarlet* constraints, such as adding an L0 penalty. This should help the code execute faster, as one of the main reasons for the slow down is unecessarily large boxes surrounding the smaller sources. See https://github.com/lsst/meas_deblender/blob/master/python/lsst/meas/deblender/deblend.py#L462-L545 for a description of the other configuration options that can be passed to the deblender\n", + "\n", + "Note: in order to try out different constraints the `meas_deblender` package requires an upgrade that has not been pushed to master yet. To use the latest changes, from your terminal session you must execute the following steps:\n", + "\n", + "```bash\n", + "~$ cp /project/fred3m/tutorials/lsst2018/.user_setups ~/notebooks\n", + "~$ source ~/notebooks/.user_setups\n", + "```\n", + "\n", + "You will have to restart the kernel for this notebook session in order for the changes to take place." + ] + } + ], + "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/Deblending/scarlet_tutorial.ipynb b/Deblending/scarlet_tutorial.ipynb new file mode 100755 index 00000000..b98af66c --- /dev/null +++ b/Deblending/scarlet_tutorial.ipynb @@ -0,0 +1,484 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Deblending with *Scarlet*\n", + "
Owner(s): **Fred Moolekamp** ([@fred3m](https://github.com/LSSTScienceCollaborations/StackClub/issues/new?body=@fred3m))\n", + "
Last Verified to Run: **2018-08-17**\n", + "
Verified Stack Release: **w_2018_32**\n", + "\n", + "The purpose of this tutorial is to familiarize you with the basics of using *scarlet* to model blended scenes, and how tweaking various objects and parameters affects the resulting model. A tutorial that is more specific to using scarlet in the context of the LSST DM Science Pipelines is also available.\n", + "\n", + "### Learning Objectives:\n", + "\n", + "After working through this tutorial you should be able to: \n", + "1. Configure and run _scarlet_ on a test list of objects;\n", + "2. Understand its various model assumptions and applied constraints.\n", + "\n", + "Before attempting this tutorial it will be useful to read the [introduction](http://scarlet.readthedocs.io/en/latest/user_docs.html) to the *scarlet* User Guide, and many of the exercises below may require referencing the *scarlet* [docs](http://scarlet.readthedocs.io/en/latest/index.html).\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": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "# Import the necessary libraries\n", + "import os\n", + "\n", + "%matplotlib inline\n", + "import matplotlib\n", + "import matplotlib.pyplot as plt\n", + "# don't interpolate the pixels\n", + "matplotlib.rc('image', interpolation='none')\n", + "\n", + "import numpy as np\n", + "import scarlet\n", + "import scarlet.display" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Display functions\n", + "\n", + "Below are several usful functions used throughout this tutorial to visualize the data and models." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "# Display the sources\n", + "def display_sources(sources, image, norm=None, subset=None, combine=False, show_sed=True, filter_indices=None):\n", + " \"\"\"Display the data and model for all sources in a blend\n", + "\n", + " This convenience function is used to display all (or a subset) of\n", + " the sources and (optionally) their SED's.\n", + " \"\"\"\n", + " if subset is None:\n", + " # Show all sources in the blend\n", + " subset = range(len(sources))\n", + " if filter_indices is None:\n", + " filter_indices = [3,2,1]\n", + " for m in subset:\n", + " # Load the model for the source\n", + " src = sources[m]\n", + " model = [comp.get_model() for comp in src]\n", + "\n", + " # Select the image patch the overlaps with the source and convert it to an RGB image\n", + " img_rgb = scarlet.display.img_to_rgb(image[src[0].bb], filter_indices=filter_indices, norm=norm)\n", + "\n", + " # Build a model for each component in the model\n", + " rgb = []\n", + " for _model in model:\n", + " # Convert the model to an RGB image\n", + " _rgb = scarlet.display.img_to_rgb(_model, filter_indices=filter_indices, norm=norm)\n", + " rgb.append(_rgb)\n", + "\n", + " # Display the image and model\n", + " figsize = [6,3]\n", + " columns = 2\n", + " # Calculate the number of columns needed and shape of the figure\n", + " if show_sed:\n", + " figsize[0] += 3\n", + " columns += 1\n", + " if not combine:\n", + " figsize[0] += 3*(len(model)-1)\n", + " columns += len(model)-1\n", + " # Build the figure\n", + " fig = plt.figure(figsize=figsize)\n", + " ax = [fig.add_subplot(1,columns,n+1) for n in range(columns)]\n", + " ax[0].imshow(img_rgb)\n", + " ax[0].set_title(\"Data: Source {0}\".format(m))\n", + " for n, _rgb in enumerate(rgb):\n", + " ax[n+1].imshow(_rgb)\n", + " if combine:\n", + " ax[n+1].set_title(\"Initial Model\")\n", + " else:\n", + " ax[n+1].set_title(\"Component {0}\".format(n))\n", + " if show_sed:\n", + " for comp in src:\n", + " ax[-1].plot(comp.sed)\n", + " ax[-1].set_title(\"SED\")\n", + " ax[-1].set_xlabel(\"Band\")\n", + " ax[-1].set_ylabel(\"Intensity\")\n", + " # Mark the current source in the image\n", + " y,x = src[0].center\n", + " ax[0].plot(x-src[0].bb[2].start, y-src[0].bb[1].start, 'x', color=\"#5af916\", mew=2)\n", + " plt.tight_layout()\n", + " plt.show()\n", + "\n", + "def display_model_residual(images, blend, peaks, norm, filter_indices=None):\n", + " \"\"\"Display the data, model, and residual for a given result\n", + " \"\"\"\n", + " if filter_indices is None:\n", + " filter_indices = [3,2,1]\n", + " model = blend.get_model()\n", + " residual = images-model\n", + " print(\"Data range: {0:.3f} to {1:.3f}\\nresidual range: {2:.3f} to {3:.3f}\\nrms: {4:.3f}\".format(\n", + " np.min(images),\n", + " np.max(images),\n", + " np.min(residual),\n", + " np.max(residual),\n", + " np.sqrt(np.std(residual)**2+np.mean(residual)**2)\n", + " ))\n", + " # Create RGB images\n", + " img_rgb = scarlet.display.img_to_rgb(images, filter_indices=filter_indices, norm=norm)\n", + " model_rgb = scarlet.display.img_to_rgb(model, filter_indices=filter_indices, norm=norm)\n", + " residual_norm = scarlet.display.Linear(img=residual)\n", + " residual_rgb = scarlet.display.img_to_rgb(residual, filter_indices=filter_indices, norm=residual_norm)\n", + "\n", + " # Show the data, model, and residual\n", + " fig = plt.figure(figsize=(15,5))\n", + " ax = [fig.add_subplot(1,3,n+1) for n in range(3)]\n", + " ax[0].imshow(img_rgb)\n", + " ax[0].set_title(\"Data\")\n", + " ax[1].imshow(model_rgb)\n", + " ax[1].set_title(\"Model\")\n", + " ax[2].imshow(residual_rgb)\n", + " ax[2].set_title(\"Residual\")\n", + " for k,component in enumerate(blend.components):\n", + " y,x = component.center\n", + " #px, py = peaks[k]\n", + " ax[0].plot(x, y, \"gx\")\n", + " #ax[0].plot(px, py, \"rx\")\n", + " ax[1].text(x, y, k, color=\"r\")\n", + " plt.show()\n", + "\n", + "def show_psfs(psfs, filters, norm=None):\n", + " rows = int(np.ceil(len(psfs)/3))\n", + " columns = min(len(psfs), 3)\n", + " figsize = (45/columns, rows*5)\n", + " fig = plt.figure(figsize=figsize)\n", + " ax = [fig.add_subplot(rows, columns, n+1) for n in range(len(psfs))]\n", + " for n, psf in enumerate(psfs):\n", + " im = ax[n].imshow(psf, norm=norm)\n", + " ax[n].set_title(\"{0}-band PSF\".format(filters[n]))\n", + " plt.colorbar(im, ax=ax[n])\n", + " plt.show()\n", + "\n", + "def display_diff_kernels(psf_blend, diff_kernels):\n", + " model = psf_blend.get_model()\n", + " for b, component in enumerate(psf_blend.components):\n", + " fig = plt.figure(figsize=(15,2.5))\n", + " ax = [fig.add_subplot(1,4,n+1) for n in range(4)]\n", + " # Display the psf\n", + " ax[0].set_title(\"psf\")\n", + " _img = ax[0].imshow(psfs[b])\n", + " fig.colorbar(_img, ax=ax[0])\n", + " # Display the model\n", + " ax[1].set_title(\"modeled psf\")\n", + " _model = np.ma.array(model[b], mask=model[b]==0)\n", + " _img = ax[1].imshow(_model)\n", + " fig.colorbar(_img, ax=ax[1])\n", + " # Display the difference kernel\n", + " ax[2].set_title(\"difference kernel\")\n", + " _img = ax[2].imshow(np.ma.array(diff_kernels[b], mask=diff_kernels[b]==0))\n", + " fig.colorbar(_img, ax=ax[2])\n", + " # Display the residual\n", + " ax[3].set_title(\"residual\")\n", + " residual = psfs[b]-model[b]\n", + " vabs = np.max(np.abs(residual))\n", + " _img = ax[3].imshow(residual, vmin=-vabs, vmax=vabs, cmap='seismic')\n", + " fig.colorbar(_img, ax=ax[3])\n", + " plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Load and Display the data\n", + "\n", + "The `file_path` points to a directory with 147 HSC blends from the COSMOS field detected by the LSST pipeline. Changing `idx` below will select a different blend." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Load the sample images\n", + "idx = 53\n", + "file_path = \"/project/shared/data/testdata_deblender/real_data/hsc_cosmos/not_matched\"\n", + "files = os.listdir(file_path)\n", + "data = np.load(os.path.join(file_path, files[idx]))\n", + "image = data[\"images\"]\n", + "wmap = data[\"weights\"]\n", + "peaks = data[\"peaks\"]\n", + "psfs = data[\"psfs\"]\n", + "filters = [\"G\", \"R\", \"I\", \"Z\", \"Y\"]\n", + "# Only a rough estimate of the background is needed\n", + "# to initialize and resize the sources\n", + "bg_rms = np.std(image, axis=(1,2))\n", + "print(\"Background RMS: {0}\".format(bg_rms))\n", + "\n", + "# Use Asinh scaling for the images\n", + "norm = scarlet.display.Asinh(img=image, Q=10)\n", + "# Map i,r,g -> RGB\n", + "filter_indices = [3,2,1]\n", + "# Convert the image to an RGB image\n", + "img_rgb = scarlet.display.img_to_rgb(image, filter_indices=filter_indices, norm=norm)\n", + "plt.imshow(img_rgb)\n", + "plt.title(\"Image: {0}\".format(idx))\n", + "for src in peaks:\n", + " plt.plot(src[0], src[1], \"rx\", mew=2)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Initializing Sources\n", + "\n", + "Astrophysical objects are modeled in scarlet as a collection of components, where each component has a single SED that is constant over it's morphology (band independent intensity). So a single source might have multiple components, like a bulge and disk, or a single component.\n", + "\n", + "The different classes that inherit from `Source` mainly differ in how they are initialized, and otherwise behave similarly during the optimization routine. This section illustrates the differences between different source initialization classes.\n", + "\n", + "The simplest source is a single component intialized with only a single pixel (at the center of the object) turned on.\n", + "\n", + "### *WARNING* \n", + "Scarlet accepts source positions using the numpy/C++ convention of (y,x), which is different than the astropy and LSST stack convention of (x,y)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sources = [scarlet.PointSource((peak[1], peak[0]), image) for peak in peaks]\n", + "\n", + "# Display the initial guess for each source\n", + "display_sources(sources, image, norm=norm)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Exercise:\n", + "\n", + "* Experiment with the above code by using `ExtendedSource`, which initializes each object as a single component with maximum flux at the peak that falls off monotonically and has 180 degree symmetry; and using `MultiComponentSource`, which models a source as two components (a bulge and a disk) that are each symmetric and montonically decreasing from the peak.\n", + "\n", + "# Deblending a scene\n", + "\n", + "The `Blend` class contains the list of sources, the image, and any other configuration parameters necessary to fit the data, including routines to fit the center positions and resize the bounding box containing the sources (if necessary). Once a blend has been initialized with a list of sources, the image and background RMS values must be set (the background RMS is used to determine when to truncate the bounding box around a source)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "blend = scarlet.Blend(sources)\n", + "blend.set_data(image, bg_rms=bg_rms)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next we can fit a model, given a maximum number of iterations and the relative error required for convergence." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "blend.fit(100, 1e-2)\n", + "print(\"Deblending completed in {0} iterations\".format(blend.it))\n", + "display_model_residual(image, blend, peaks, norm)\n", + "display_sources(sources, image, norm)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Exercises\n", + "\n", + "* Experiment by running the above code using different source models (for example `ExtendedSource`) to see how initializtion affects the belnding results.\n", + "\n", + "* The code above initialized the sources at their exact centers. Try offsetting the initial positions by `0.5` pixels in `x` and/or `y` and passing a `shift_center=0` argument when initializing the source. This prevents the source from updating its position, so notice how that affects the resulting model." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Constraints\n", + "\n", + "The above models used the default constraints: perfect symmetry and a weighted monotonicity that decreases from the peak. So each source is defined (internally during initialization) with the constraints" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "import scarlet.constraint as sc\n", + "constraints = (sc.SimpleConstraint(),\n", + " sc.DirectMonotonicityConstraint(use_nearest=False),\n", + " sc.DirectSymmetryConstraint())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "where `SimpleConstraint` forces the SED and morphology to be non-negative, the SED to be normalized to unity, and the peak to have some (minimal) flux at the center." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Exercises\n", + "\n", + "* Go back to the source initialization cell and pass a custom set of constraints. For example, pass `DirectSymmetryConstraint` a number between 0 and 1 to set the level of symmetry required, or eliminate the symmetry constraint altogether and see how that affects deblending.\n", + "\n", + "* Set `use_nearest=True` in the `DirectMonotonicityConstraint`.\n", + "\n", + "* Add `L0Constraint` or `L1Constraint` to the list of constraints and observe the results." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Configuration\n", + "\n", + "There are additional configuration paramters that can be used to initialize a source, as described in http://scarlet.readthedocs.io/en/latest/config.html#Configuration-(scarlet.config).\n", + "\n", + "## Exercises\n", + "\n", + "* Initialize the sources with a custom configuration where `refine_skip=2`, which updates the positions and box sizes on every other step, and see how the results compare" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# PSF Deconvolution\n", + "\n", + "When analyzing real images the PSF will be different in each band unless they have been PSF matched. In general deblending should not be performed on PSF matched coadds, as matching will increase the blending in bands with better seeing. Instead scarlet can be used to build a deconvolved model which is a more sparse (and less blended) representation of the data, and convolve the model in each band to compare to the input data.\n", + "\n", + "To initialize a source with a PSF, pass the PSF as an input to the new source:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "scarlet.ExtendedSource(peaks[0], image, bg_rms, psf=psfs)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Partial PSF Deconvolution\n", + "\n", + "As discussed in the tutorial http://scarlet.readthedocs.io/en/latest/psf_matching.html, the data is noisy and the fully deconvolved scene is undersampled, making the application of the constraints and and full convolution kernel unstable and prone to biases. Instead we can create a target PSF and model the sources in the partially deconvolved target PSF scene.\n", + "\n", + "First we need to specify the target PSF. *scarlet* includes a `fit_target_psf` function to fit the PSF in each band to either a `moffat`, `gaussian`, or `double_gaussian` function. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import scarlet.psf_match\n", + "\n", + "show_psfs(psfs, filters)\n", + "\n", + "# Find the target PSF\n", + "target_psf = scarlet.psf_match.fit_target_psf(psfs, scarlet.psf_match.moffat)\n", + "plt.imshow(target_psf)\n", + "plt.title(\"target PSF\")\n", + "plt.colorbar()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Once we have the target PSF we can find the difference kernel in each band using *scarlet*. The `build_diff_kernels` function basically treats the PSF image as a blend, where the PSF in each band is a monochromatic source, and fits the difference kernels using the minimum number of pixels necessary." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "diff_kernels, psf_blend = scarlet.psf_match.build_diff_kernels(psfs, target_psf)\n", + "display_diff_kernels(psf_blend, diff_kernels)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Exercises\n", + "\n", + "* Try building the difference kernels while varying the parameters in `build_diff_kernels`, for example using larger and smaller values for `l0_thresh`.\n", + "\n", + "* Go back up to source initialization and use `psf=psfs` to fully deconvolve the scene and fit the blend\n", + "\n", + "* Try the same thing but set `psf=diff_kernels` for each source to partially deconvolve the scene." + ] + } + ], + "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/GettingStarted/FindingDocs.ipynb b/GettingStarted/FindingDocs.ipynb index 1fd994da..ca7c8947 100644 --- a/GettingStarted/FindingDocs.ipynb +++ b/GettingStarted/FindingDocs.ipynb @@ -27,11 +27,11 @@ "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", + "We'll need the `stackclub` package to be installed. If you are not developing this package, and you have permission to write to your base python site-packages, 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:" + "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, like this:" ] }, { diff --git a/GettingStarted/GettingStarted.md b/GettingStarted/GettingStarted.md index 9dd0878c..46c9fe6f 100644 --- a/GettingStarted/GettingStarted.md +++ b/GettingStarted/GettingStarted.md @@ -4,20 +4,20 @@ _Greg Madejski and [Phil Marshall](https://github.com/LSSTScienceCollaborations/ We are developing tutorial notebooks on remote JupyterLab instances, to short-circuit the DM stack installation process and get used to working in the notebook aspect of the LSST science platform. In these notes we provide: -* [Notes on how to get set up on the LSST Science Platform (LSP) JupyterLab Notebook Aspect at NCSA](https://github.com/LSSTScienceCollaborations/StackClub/blob/master/GettingStarted/GettingStarted.md#accessing-the-lsst-science-platform) +* [Notes on how to get set up on the LSST Science Platform (LSP) JupyterLab Notebook Aspect at the LSST Data Facility at NCSA](https://github.com/LSSTScienceCollaborations/StackClub/blob/master/GettingStarted/GettingStarted.md#accessing-the-lsst-science-platform) * [Help with getting set up to run and edit the Stack Club tutorial notebooks](https://github.com/LSSTScienceCollaborations/StackClub/blob/master/GettingStarted/GettingStarted.md#running-and-contributing-to-the-stack-club-notebooks) ## Accessing the LSST Science Platform -The [LSST Science Platform (LSP) Notebook Aspect Documentation](https://nb.lsst.io/) provides an introduction to the NCSA system, including how to gain access and then how to use JupyterLab once you are in. -Getting in to NCSA takes involves getting an NCSA account, and then figuring out VPN access. +The [LSST Science Platform (LSP) Notebook Aspect Documentation](https://nb.lsst.io/) provides an introduction to the system, including how to gain access and then how to use JupyterLab once you are in. +Getting on to the LSP involves getting an NCSA account, and then figuring out VPN access. -#### Getting an NCSA Account -Contact Phil (DM @drphilmarshall on LSSTC Slack) to get an NCSA Stack Club account. You'll need to provide your full name (first and last) and your email address. You'll (eventually) get an email invitation to [create an account at NCSA](https://identity.ncsa.illinois.edu/) (including choosing a username of 8 characters or fewer). After you have submitted your form, it typically takes 24 hours for your account to be set up: set an alarm to come back the next day! +#### Getting an LSST Science Platform Account +The Stack Club has a limited number of active LSST Science Platform accounts it can support. To join the Stack Club and request one of these accounts, please fill out the [Stack Club Membership Application Form](https://goo.gl/forms/588KlPTFfkEEFFUu2). You'll need to agree to abide by the [Rules](../Rules.md), and then provide your full name (first and last) and your email address. If your application is successful, you'll get an email with instructions on how to set up your LSP account. -#### Accessing NCSA via its VPN +#### Accessing the LSP via its VPN At present, unless you are on an approved network, you must use the [NCSA virtual private network (VPN)](https://wiki.ncsa.illinois.edu/display/cybersec/Virtual+Private+Network+%28VPN%29+Service). -The recommended method is to use Cisco's AnyConnect with DUO two-factor authentication. Detailed instructions are available on the [NCSA VPN site](https://wiki.ncsa.illinois.edu/display/cybersec/Virtual+Private+Network+%28VPN%29+Service#VirtualPrivateNetwork(VPN)Service-UsingtheCiscoAnyConnectVPNClient(Required)). +The recommended method is to use Cisco's AnyConnect with DUO two-factor authentication (verified on Mac and Linux). Detailed instructions are available on the [NCSA VPN site](https://wiki.ncsa.illinois.edu/display/cybersec/Virtual+Private+Network+%28VPN%29+Service#VirtualPrivateNetwork(VPN)Service-UsingtheCiscoAnyConnectVPNClient(Required)). > You can get AnyConnect by pointing your browser at https://sslvpn.ncsa.illinois.edu/ and selecting the `ncsa-vpn-default` option (this will only work if you have a java-compatible browser, like firefox esr version<=52). If you already have the AnyConnect client installed, open it up and enter `sslvpn.ncsa.illinois.edu/` in its connection window. @@ -25,6 +25,8 @@ The recommended method is to use Cisco's AnyConnect with DUO two-factor authenti If you forget your password it can be reset following the instructions [here](https://developer.lsst.io/services/lsst-dev.html?highlight=reset#lsst-dev-password). If you have problems connecting to the NCSA services you can check their status and submit a help ticket [here](https://confluence.lsstcorp.org/display/DM/LSST+Service+Status+page). +For a Linux install, you may need to pre-install [`openconnect`](http://www.infradead.org/openconnect/) from your favorite package manager. + #### Starting up the LSST Science Platform JupyterLab Notebook Aspect Once the VPN connection is established, you should be able to navigate to the the JupyterLab instance at **https://lsst-lspdev.ncsa.illinois.edu/nb**. Select the `Release` and `medium` options on the Spawner Options landing page, and then hit the "Spawn" button. You'll (eventually) end up on the JupyterLab launcher, where you can use the file manager in the left hand side bar to open your Jupyter notebooks, or start terminal or notebook editor tabs from the buttons provided. You should see the pre-installed `notebook-demo` notebooks in the file manager, for example. @@ -61,18 +63,20 @@ Broadly useful, small datasets are available in `/project/shared/data` - this i Larger datasets are available in `/datasets`. This is a read-only folder. #### The Stack Club Library -The [`stackclub` folder in this repo](../stackclub) is a python package containing a number of utility functions and classes for use in tutorial notebooks. You can browse its documentation at https://stackclub.readthedocs.io/. If you are not developing this package, you can install it using pip, like this: -``` -pip install git+git://github.com/LSSTScienceCollaborations/StackClub.git#egg=stackclub -``` -However, if you are contributing notebooks it is likely that you'll need to develop the `stackclub` package as well -(eg by adding modules to it), and so you'll need to make a local, editable installation. In the top level folder of your local clone of the StackClub repo, do: +The [`stackclub` folder in this repo](../stackclub) is a python package containing a number of utility functions and classes for use in tutorial notebooks. You can browse its documentation at https://stackclub.readthedocs.io/. +If you are contributing notebooks, you may want or need to develop the `stackclub` package as well +(eg by adding modules to it), and so its best to setup the package installation to be local and editable. In the top level folder of your local clone of the StackClub repo, do: ``` python setup.py -q develop --user ``` -This will put the `stackclub` folder on your path. You may find the following lines useful to add to your notebook as you develop the library: +This will put the repo's `stackclub` folder on your path. When developing the package, you may find it useful to add the following lines to your notebook: ```python %load_ext autoreload %autoreload 2 ``` -This enables you to repeatedly `import stackclub` as you update the library code. +This enables you to repeatedly `import stackclub` as you update the library code. The above lines are in the [template notebook](templates/template_Notebook.ipynb), for your convenience. + +If you are not developing this package, and you have permission to write to your base python site-packages, you can install it using pip, like this: +``` +pip install git+git://github.com/LSSTScienceCollaborations/StackClub.git#egg=stackclub +``` diff --git a/GettingStarted/HelloWorld.ipynb b/GettingStarted/HelloWorld.ipynb index ee48b344..8fccd6fd 100644 --- a/GettingStarted/HelloWorld.ipynb +++ b/GettingStarted/HelloWorld.ipynb @@ -277,6 +277,22 @@ "source": [ "take_that_first_baby_step(and_follow_up=\"it's me, Phil, using python!\")" ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(\"Rob was here, thanks all!\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { diff --git a/ImageProcessing/BrighterFatterCorrection.ipynb b/ImageProcessing/BrighterFatterCorrection.ipynb new file mode 100644 index 00000000..abcb766e --- /dev/null +++ b/ImageProcessing/BrighterFatterCorrection.ipynb @@ -0,0 +1,754 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Analysis of Beam Simulator Images and Brighter-fatter Correction\n", + "
Owner(s): **Andrew Bradshaw** ([@andrewkbradshaw](https://github.com/LSSTScienceCollaborations/StackClub/issues/new?body=@andrewkbradshaw))\n", + "
Last Verified to Run: **2018-09-14**\n", + "
Verified Stack Release: **16.0 and 16.0+22 (w_2018_31)**\n", + "\n", + "This notebook demonstrates the [brighter-fatter systematic error](https://arxiv.org/abs/1402.0725) on images of stars and galaxies illuminated on an ITL-3800C-002 CCD at the [UC Davis LSST beam simulator laboratory](https://arxiv.org/abs/1411.5667). Using a series of images at increasing exposure times, we demonstrate the broadening of image profiles on DM stack shape measurements, and a [possible correction method](https://arxiv.org/abs/1711.06273) which iteratively applies a kernel to restore electrons to the pixels from which they were deflected. To keep things simple, for now we skip most DM stack instrument signature removal (ISR) and work on a subset of images which are already processed arrays (500x500) of electrons.\n", + "\n", + "### Learning Objectives:\n", + "\n", + "After working through this tutorial you should be able to: \n", + "1. Characterize and measure objects (stars/galaxies) in LSST beam simulator images\n", + "2. Test the Brighter-Fatter kernel correction method on those images\n", + "3. Build your own tests of stack ISR algorithms\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": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import pandas as pd\n", + "import matplotlib.pyplot as plt\n", + "from matplotlib.colors import LogNorm\n", + "from itertools import cycle\n", + "from astropy.io import fits\n", + "import time,glob,os\n", + "\n", + "\n", + "# if running stack v16.0, silence a long matplotlib Agg warning with:\n", + "import warnings\n", + "warnings.filterwarnings(\"ignore\", category=UserWarning)\n", + "\n", + "%matplotlib inline\n", + "\n", + "# What version of the Stack am I using?\n", + "! echo $HOSTNAME\n", + "! eups list -s | grep lsst_distrib\n", + "\n", + "# make a directory to write the catalogs\n", + "username=os.environ.get('USERNAME')\n", + "cat_dir='/home/'+username+'/DATA/beamsim/'\n", + "if not os.path.exists(cat_dir):\n", + " ! mkdir /home/$USER/DATA/beamsim/" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Step 1: Read in an image\n", + "Cut-outs of beam simulator star/galaxy images have been placed in the shared data directory at `/project/shared/data/beamsim/bfcorr/`. We skip (for now) most of the instrument signature removal (ISR) steps because these are preprocessed images (bias subtracted, gain corrected). We instead start by reading in one of those `.fits` files and making an image plane `afwImage.ExposureF` as well as a variance plane (based upon the image), which is then ready for characterization and calibration in the following cells." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import lsst.afw.image as afwImage\n", + "from lsst.ip.isr.isrFunctions import updateVariance\n", + "\n", + "# where the data lives, choosing one image to start\n", + "imnum=19 # for this dataset, choose 0-19 as an example\n", + "fitsglob='/project/shared/data/beamsim/bfcorr/*part.fits'\n", + "fitsfilename = np.sort(glob.glob(fitsglob))[imnum] \n", + "\n", + "# Read in a single image to an afwImage.ImageF object\n", + "image_array=afwImage.ImageF.readFits(fitsfilename)\n", + "image = afwImage.ImageF(image_array)\n", + "exposure = afwImage.ExposureF(image.getBBox())\n", + "exposure.setImage(image)\n", + "hdr=fits.getheader(fitsfilename) # the header has some useful info in it\n", + "print(\"Read in \",fitsfilename.split('/')[-1])\n", + "\n", + "# Set the variance plane using the image plane via updateVariance function\n", + "gain = 1.0 # because these images are already gain corrected\n", + "readNoise = 10.0 # in electrons\n", + "updateVariance(exposure.maskedImage, gain, readNoise)\n", + "\n", + "# Another way of setting variance and/or masks?\n", + "#mask = afwImage.makeMaskFromArray(np.zeros((4000,4072)).astype('int32'))\n", + "#variance = afwImage.makeImageFromArray((readNoise**2 + image_array.array())\n", + "#masked_image = afwImage.MaskedImageF(image, mask, variance)\n", + "#exposure = afwImage.ExposureF(masked_image)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now visualize the image and its electron distribution using matplotlib. Things to note: 1) the array is (purposefully) tilted with respect to the pixel grid, 2) most pixel values are at the background/sky level (a function of the mask opacity and illumination), but there is a pileup of counts around ~200k electrons indicating full well and saturation in some of the brightest pixels of the image" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure(figsize=(12,5)),plt.subplots_adjust(wspace=.3)\n", + "plt.suptitle('Star/galaxy beam sim image and histogram \\n'+fitsfilename.split('/')[-1])\n", + "\n", + "plt.subplot(121)\n", + "plt.imshow(exposure.image.array,vmax=1e3,origin='lower')\n", + "plt.colorbar(label='electrons')\n", + "\n", + "plt.subplot(122)\n", + "plt.hist(exposure.image.array.flatten(),bins=1000,histtype='step')\n", + "plt.yscale('log')#,plt.xscale('log')\n", + "plt.xlabel('Number of electrons in pixel'),plt.ylabel('Number of pixels')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Step 2: Perform image characterization and initial measurement\n", + "We now perform a base-level characterization of the image using the stack. We set some configuration settings which are specific to our sestup which has a very small optical PSF, setting a PSF size and turning off some other aspects such as cosmic ray rejection because of this." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from lsst.pipe.tasks.characterizeImage import CharacterizeImageTask, CharacterizeImageConfig\n", + "import lsst.meas.extensions.shapeHSM\n", + "\n", + "# first set a few configs that are specific to our beam simulator data\n", + "charConfig = CharacterizeImageConfig()\n", + "#this set the fwhm of the simple PSF to that of optics\n", + "charConfig.installSimplePsf.fwhm = .2\n", + "charConfig.doMeasurePsf = False\n", + "charConfig.doApCorr = False # necessary\n", + "charConfig.repair.doCosmicRay = False \n", + "# we do have some cosmic rays, but we also have subpixel mask features and an undersampled PSF\n", + "charConfig.detection.background.binSize = 10 # worth playing around with\n", + "#charConfig.background.binSize = 50\n", + "charConfig.detection.minPixels = 2 # also worth playing around with\n", + "\n", + "# Add the HSM (Hirata/Seljak/Mandelbaum) adaptive moments shape measurement plugin\n", + "charConfig.measurement.plugins.names |= [\"ext_shapeHSM_HsmSourceMoments\"]\n", + "# to configure hsm you would do something like\n", + "# charConfig.measurement.plugins[\"ext_shapeHSM_hsmSourceMoments\"].addFlux = True\n", + "# (see sfm.py in meas_base for all the configuration options for the measurement task)\n", + "\n", + "charTask = CharacterizeImageTask(config=charConfig)\n", + "\n", + "charTask.characterize?\n", + "# use charTask.run instead of characterize for v16.0+22\n", + "# could also perform similar functions with processCcdTask.run()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Display which plugins are being used for measurement\n", + "charConfig.measurement.plugins.active " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "tstart=time.time()\n", + "charResult = charTask.characterize(exposure) # charTask.run(exposure) stack v16.0+22\n", + "print(\"Characterization took \",str(time.time()-tstart)[:4],\" seconds\")\n", + "print(\"Detected \",len(charResult.sourceCat),\" objects \")\n", + "\n", + "plt.title('X/Y locations of detections')\n", + "plt.plot(charResult.sourceCat['base_SdssCentroid_x'],charResult.sourceCat['base_SdssCentroid_y'],'r.')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This figure illustrates the centroids of detections made during characterization. Note that most objects have been detected (except for one), and that there are several spurious detections which are not on our grid. Further visualization of these will be done in the Firefly window a few cells below." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# display some of the source catalog measurements filtered by searchword\n", + "searchword='shape'\n", + "for name in charResult.sourceCat.schema.getOrderedNames():\n", + " if searchword in name.lower():\n", + " print(name)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Looking at the mask plane, which started off as all zeros\n", + "# and now has some values of 2^5\n", + "maskfoo=exposure.mask\n", + "print(\"Unique mask plane values: \",np.unique(maskfoo.array))\n", + "print(\"Mask dictionary entries: \",maskfoo.getMaskPlaneDict())\n", + "\n", + "plt.figure(figsize=(12,5)),plt.subplots_adjust(wspace=.3)\n", + "plt.subplot(121)\n", + "plt.imshow(maskfoo.array,origin='lower'),plt.colorbar()\n", + "plt.subplot(122)\n", + "plt.hist(maskfoo.array.flatten()),plt.xlabel('Mask plane values')\n", + "plt.yscale('log')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The above figures illustrate the new mask plane of the exposure object which was added and modified during characterization. Values of 0 and 5 are now seen, which correspond to unassociated pixels and those which are \"detected\". Further visualization of the mask plane can be seen in the Firefly cell down below." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Step 3: Further image calibration and measurement\n", + "This builds on the exposure output from characterization, using the new mask plane as well as the source catalog. Similar to the characterization, we turn off some processing which is suited to our particular setup. For this dataset a calibrate task is almost unncessary (as it is not on-sky data and we don't have a reference catalog), however, it does provide a background-subtracted image and for completeness it is included here. The steps in calibration that are turned on/off can be seen by printing the calibration config object." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from lsst.pipe.tasks.calibrate import CalibrateTask, CalibrateConfig\n", + "\n", + "calConfig = CalibrateConfig()\n", + "calConfig.doAstrometry = False\n", + "calConfig.doPhotoCal = False\n", + "calConfig.doApCorr = False\n", + "calConfig.doDeblend = False # these are well-separated objects, deblending adds time & trouble\n", + "# these images should have a uniform background, so measure it\n", + "# on scales which are larger than the objects\n", + "calConfig.detection.background.binSize = 50\n", + "calConfig.detection.minPixels = 15\n", + "calConfig.measurement.plugins.names |= [\"ext_shapeHSM_HsmSourceMoments\"]\n", + "# to configure hsm you would do something like\n", + "#charConfig.measurement.plugins[\"ext_shapeHSM_hsmSourceMoments\"].addFlux = True\n", + "\n", + "calTask = CalibrateTask(config= calConfig, icSourceSchema=charResult.sourceCat.schema)\n", + "\n", + "#calTask.run? # for stack v16.0+22 \n", + "calTask.calibrate?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "tstart=time.time()\n", + "# for stack v16.0+22, change to calTask.run(charResult.exposure)\n", + "calResult = calTask.calibrate(charResult.exposure, background=charResult.background,\n", + " icSourceCat = charResult.sourceCat)\n", + "\n", + "print(\"Calibration took \",str(time.time()-tstart)[:4],\" seconds\")\n", + "print(\"Detected \",len(calResult.sourceCat),\" objects \")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Below we look at some of the measurements in the source catalog which has been attached to the calibration result. We also save the source catalog to `$fitsfilename.cat` in `/home/$USER/beamsim/`, which was created in the first cell. This will allow the results from each image to be read in after these measurements are performed on each image." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(\"Catalogs will be saved to: \"+cat_dir)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "src=calResult.sourceCat #.copy(deep=True) ?\n", + "#print(src.asAstropy)\n", + "\n", + "# catalog directory\n", + "src.writeFits(cat_dir+fitsfilename.split('/')[-1].replace('.fits','.cat'))\n", + "# read back in and access via:\n", + "#catalog=fits.open(fitsfilename+'.cat')\n", + "#catalog[1].data['base_SdssShape_xx'] etc.\n", + "\n", + "par_names=['base_SdssShape_xx','base_SdssShape_yy','base_SdssShape_flux']\n", + "par_mins=[0,0,0]\n", + "par_maxs=[5,5,1e6]\n", + "n_par=len(par_names)\n", + "\n", + "\n", + "plt.figure(figsize=(5*n_par,6)),plt.subplots_adjust(wspace=.25)\n", + "for par_name,par_min,par_max,i in zip(par_names,par_mins,par_maxs,range(n_par)):\n", + " plt.subplot(2,n_par,i+1)\n", + " plt.scatter(src['base_SdssCentroid_x'],src['base_SdssCentroid_y'],c=src[par_name],marker='o',vmin=par_min,vmax=par_max)\n", + " plt.xlabel('X'),plt.ylabel('Y'),plt.colorbar(label=par_name)\n", + "\n", + "\n", + " plt.subplot(2,n_par,n_par+i+1)\n", + " plt.hist(src[par_name],range=[par_min,par_max],bins=20,histtype='step')\n", + " plt.xlabel(par_name)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The above figures show the 2-dimensional distribution of detected objects measured parameter values and their histogram. By default, two shape parameters (in pixels) and a flux measurement (in electrons) are shown, but this can be modified through the `par_names` variable in the cell above. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Optional step: Display the image with Firefly and overlay detected objects\n", + "This is a nice interface for looking at measurements and images together, and it is much more powerful than demonstrated below (see other stack club notebooks for demonstration). From this display, it is clear that some objects are detected as two or more, complicating downstream measurements." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import lsst.afw.display as afwDisplay\n", + "\n", + "# Firefly client imports\n", + "from firefly_client import FireflyClient\n", + "\n", + "# Standard libraries in support of Firefly display\n", + "from urllib.parse import urlparse, urlunparse, ParseResult\n", + "from IPython.display import IFrame, display, Markdown\n", + "import os\n", + "\n", + "# Info for Firefly server connection\n", + "my_channel = '{}_test_channel'.format(os.environ['USER'])\n", + "server = 'https://lsst-lspdev.ncsa.illinois.edu'\n", + "ff='{}/firefly/slate.html?__wsch={}'.format(server, my_channel)\n", + "IFrame(ff,1000,600) # initiate the window\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# set the backend and attach to the waiting display channel\n", + "afwDisplay.setDefaultBackend('firefly')\n", + "afw_display = afwDisplay.getDisplay(frame=1, \n", + " name=my_channel)\n", + "\n", + "# Open the exposure (Firefly knows about mask planes)\n", + "afw_display.mtv(exposure)\n", + "\n", + "# Now overplot sources from the src table onto the image display using the Display’s dot method \n", + "# It is more efficient to send a batch of updates to the display, \n", + "# so we enclose the loop in a display.Buffering context, like this:\n", + "\n", + "afw_display.erase() #\n", + "\n", + "with afw_display.Buffering():\n", + " for record in src[:]:\n", + " afw_display.dot('o', record.getX(), record.getY(), size=20, ctype='orange')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Step 4: Apply the brighter-fatter kernel correction to an image\n", + "This brighter fatter correction method takes in a \"kernel\" (derived from theory or flat fields) which models the broadening of incident image profiles assuming the pixel boundary displacement can be represented as the gradient of a scalar field. Given a kernel and this assumption, the incident image profile can in theory be reconstructed using an iterative process, which we test here using our beam simulator images. See [this paper](https://arxiv.org/abs/1711.06273) and the IsrTask docstring below for more details about the theory and its assumptions. The kernel used here is not generated by the stack but rather through similar code which was written at UC Davis by Craig Lage. Future additions to the notebook will use a stack-generated kernel when available." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from lsst.ip.isr.isrTask import IsrTask # brighterFatterCorrection lives here\n", + "isr=IsrTask()\n", + "\n", + "pre_bfcorr_exposure=exposure.clone() #save a copy of the pre-bf corrected image\n", + "\n", + "isr.brighterFatterCorrection?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Read in the kernel (determined from e.g. simulations or flat fields)\n", + "kernel=fits.getdata('/project/shared/data/beamsim/bfcorr/BF_kernel-ITL_3800C_002.fits')\n", + "exposure=pre_bfcorr_exposure.clone() # save the pre-bf correction image\n", + "\n", + "# define the maximum number of iterations and threshold for differencing convergence (e-)\n", + "bf_maxiter,bf_threshold=20,10\n", + "\n", + "# Perform the correction\n", + "tstart=time.time()\n", + "isr.brighterFatterCorrection(exposure,kernel,bf_maxiter,bf_threshold,False)\n", + "print(\"Brighter-fatter correction took\",time.time()-tstart,\" seconds\")\n", + "#takes 99 seconds for 4kx4k exposure, 21x21 kernel, 20 iterations, 10 thresh\n", + "\n", + "# Plot kernel and image differences\n", + "plt.figure(),plt.title('BF kernel')\n", + "plt.imshow(kernel),plt.colorbar()\n", + "\n", + "imagediff=(pre_bfcorr_exposure.image.array-exposure.image.array)\n", + "imagediffpct=np.sum(imagediff)/np.sum(pre_bfcorr_exposure.image.array)*100.\n", + "print(str(imagediffpct)[:5],' percent change in flux')\n", + "\n", + "plt.figure(figsize=(16,10))\n", + "plt.subplot(231),plt.title('Before')\n", + "plt.imshow(pre_bfcorr_exposure.image.array,vmin=0,vmax=1e3,origin='lower'),plt.colorbar()\n", + "plt.subplot(232),plt.title('After')\n", + "plt.imshow(exposure.image.array,vmin=0,vmax=1e3,origin='lower'),plt.colorbar()\n", + "plt.subplot(233),plt.title('Before - After')\n", + "vmin,vmax=-50,50\n", + "plt.imshow(imagediff,vmin=vmin,vmax=vmax,origin='lower'),plt.colorbar()\n", + "\n", + "nbins=1000\n", + "plt.subplot(234)\n", + "plt.hist(pre_bfcorr_exposure.image.array.flatten(),bins=nbins,histtype='step',label='before')\n", + "plt.yscale('log')\n", + "plt.subplot(235)\n", + "plt.hist(exposure.image.array.flatten(),bins=nbins,histtype='step',label='after')\n", + "plt.yscale('log')\n", + "plt.subplot(236)\n", + "plt.hist(imagediff.flatten(),bins=nbins,histtype='step',label='difference')\n", + "plt.yscale('log')\n", + "plt.legend()\n", + "plt.xlabel('Pixel values [e-]')\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The above figures illustrate the way that the brighter-fatter correction works: by iteratively convolving a physically-motivated kernel with the electron image to redistribute charge from the periphery to the center of objects." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Step 5: Run the above steps (with and without brighter-fatter correction) on 20 exposures of increasing exposure time\n", + "Here we re-do all of the previous work, which was done with one image, on a series of images with increasing exposure times. We will generate this series of catalogs both with and without applying the brighter-fatter correction, allowing us to test the fidelity of the brighter-fatter correction with our beam simulator images. To do this in a simple way, we create a function to perform all of the above tasks, called `make_bf_catalogs`, which only takes in a list of filenames but uses some of the same global configuration values (`charTask.config` and `calTask.config`) which we set above." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fitsglob='/project/shared/data/beamsim/bfcorr/*part.fits' # fits filenames to read in\n", + "fitsfilelist=np.sort(glob.glob(fitsglob))\n", + "\n", + "def make_bf_catalogs(fitsfilelist,do_bf_corr=False,do_verbose_print=True):\n", + " for fitsfilename in fitsfilelist:\n", + " tstart=time.time()\n", + " image_array=afwImage.ImageF.readFits(fitsfilename)\n", + " image = afwImage.ImageF(image_array)\n", + "\n", + " exposure = afwImage.ExposureF(image.getBBox())\n", + " exposure.setImage(image)\n", + "\n", + " updateVariance(exposure.maskedImage, gain, readNoise)\n", + " \n", + " # start the characterization and measurement, \n", + " # optionally beginning with the brighter-fatter correction\n", + " if do_bf_corr:\n", + " isr.brighterFatterCorrection(exposure,kernel,bf_maxiter,bf_threshold,False)\n", + " # print(\"Brighter-fatter correction took\",str(time.time()-tstart)[:4],\" seconds\")\n", + " # for stack v16.0+22 use charTask.run() and calTask.run()\n", + " charResult = charTask.characterize(exposure) \n", + " calResult = calTask.calibrate(charResult.exposure, background=charResult.background,\n", + " icSourceCat = charResult.sourceCat)\n", + " src=calResult.sourceCat\n", + "\n", + " # write out the source catalog, appending -bfcorr for the corrected catalogs\n", + " catfilename=cat_dir+fitsfilename.replace('.fits','.cat').split('/')[-1]#\n", + " if do_bf_corr: catfilename=catfilename.replace('.cat','-bfcorr.cat')\n", + " src.writeFits(catfilename)\n", + "\n", + " if do_verbose_print: \n", + " print(fitsfilename.split('/')[-1],\" char. & calib. took \",\n", + " str(time.time()-tstart)[:4],\" seconds to measure \",\n", + " len(calResult.sourceCat),\" objects \")\n", + "\n", + " \n", + "# Run the catalog maker on the series of uncorrected and corrected images\n", + "make_bf_catalogs(fitsfilelist,do_bf_corr=True,do_verbose_print=True)\n", + "make_bf_catalogs(fitsfilelist,do_bf_corr=False,do_verbose_print=False)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now read in those catalogs, both corrected and uncorrected (this could be improved with e.g. pandas)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "cat_arr = []\n", + "catglob=cat_dir+'ITL*part.cat' # uncorrected catalogs\n", + "for catfilename in np.sort(glob.glob(catglob)): cat_arr.append(fits.getdata(catfilename))\n", + "\n", + "bf_cat_arr = []\n", + "catglob=cat_dir+'ITL*part-bfcorr.cat' # corrected catalogs\n", + "for catfilename in np.sort(glob.glob(catglob)): bf_cat_arr.append(fits.getdata(catfilename))\n", + "ncats=len(cat_arr)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Show issues with multiply detected sources which which we remedy with matching rejection\n", + "for i in range(ncats):\n", + " xfoo,yfoo=cat_arr[i]['base_SdssCentroid_x'],cat_arr[i]['base_SdssCentroid_y']\n", + " plt.plot(xfoo,yfoo,'o',alpha=.4)\n", + "plt.title('Centroids of sequential exposures')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The above image illustrates a problem with comparing images with different exposure times. Namely, that different sets of objects may be detected. To remedy this, we use a fiducial frame as reference and simply match the catalogs by looking for *single* object matches within a specified distance of those fiducial objects. We then collect a shape measurement (e.g. `base_SdssShape_xx/yy`) for that object as well as a brightness measurement (e.g. `base_SdssShape_flux`) to test for a trend in size vs. brightness." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fidframe=10 # frame number to compare to\n", + "maxdist=.5 # max distance to match objects between frames\n", + "\n", + "# choose which stack measurements to use for centroids and shape\n", + "# +TODO use 'ext_shapeHSM_HsmSourceMoments_xx','ext_shapeHSM_HsmSourceMoments_yy'\n", + "cen_param1,cen_param2='base_SdssCentroid_x','base_SdssCentroid_y'\n", + "bf_param1,bf_param2='base_SdssShape_xx','base_SdssShape_yy'\n", + "flux_param='base_GaussianFlux_flux' # could try 'base_CircularApertureFlux_17_0_flux' or 'base_SdssShape_flux'\n", + "\n", + "# get the centroids (used for matching) from the fiducial frame \n", + "x0s,y0s=cat_arr[fidframe][cen_param1],cat_arr[fidframe][cen_param2]\n", + "nspots=len(x0s)\n", + "\n", + "# make an array to hold that number of objects and their centroid/shape/flux measurements\n", + "# the 8 rows collect x/y centroid, x/y shape, x/y corrected shape, flux, and corrected flux\n", + "bf_dat=np.empty((ncats,nspots),\n", + " dtype=np.dtype([('x', float), ('y', float),('shapex', float), ('shapey', float),\n", + " ('corrshapex', float), ('corrshapey', float),\n", + " ('flux', float), ('corrflux', float)]))\n", + "bf_dat[:]=np.nan # so that un-matched objects aren't plotted/used by default\n", + "\n", + "\n", + "# loop over catalogs\n", + "for i in range(ncats):\n", + " # get the centroids of objects in the bf-corrected and uncorrected images\n", + " x1,y1=cat_arr[i][cen_param1],cat_arr[i][cen_param2]\n", + " x1_bf,y1_bf=bf_cat_arr[i][cen_param1],bf_cat_arr[i][cen_param2]\n", + " # loop over fiducial frame centroids to find matches\n", + " for j in range(nspots): \n", + " x0,y0=x0s[j],y0s[j] # fiducial centroid to match\n", + " # find objects in both catalogs which are within maxdist\n", + " bf_gd=np.where(np.sqrt((x1_bf-x0)**2+(y1_bf-y0)**2)250 - sz ) | (np.abs(yc-250)>250 - sz )): continue\n", + " stamp=exposure.getImage().array[yc-sz:yc+sz,xc-sz:xc+sz]\n", + " \n", + " # show the stamp with log scale (1,max)\n", + " plt.subplot(131),plt.title('stamp '+str(index).zfill(3)+' (log scale)')\n", + " plt.imshow(stamp,origin='lower',norm=LogNorm(1,stamp.max())),plt.colorbar()\n", + " \n", + " # x size vs flux\n", + " plt.subplot(132),plt.title('x (row) size vs. flux')\n", + " plt.plot(bf_dat['flux'][:,index],bf_dat['shapex'][:,index],'r.',label='Uncorrected')\n", + " plt.plot(bf_dat['corrflux'][:,index],bf_dat['corrshapex'][:,index],'g.',label='Corrected')\n", + " plt.xlabel(flux_param),plt.ylabel(bf_param1),plt.xscale('log')\n", + " plt.legend(loc='upper left')\n", + "\n", + " # y size vs flux\n", + " plt.subplot(133),plt.title('y (column) size vs. flux')\n", + " plt.plot(bf_dat['flux'][:,index],bf_dat['shapey'][:,index],'r.',label='Uncorrected')\n", + " plt.plot(bf_dat['corrflux'][:,index],bf_dat['corrshapey'][:,index],'g.',label='Corrected')\n", + " plt.xlabel(flux_param),plt.ylabel(bf_param2)\n", + " plt.xscale('log')\n", + " plt.savefig(cat_dir+str(index).zfill(5)+'bfcorr.png')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The above figures illustrate the brighter-fatter effect (slight increasing size with flux) in the red dots, and the corrected image analysis in green. Curiously, some of the objects indicate that the default correction method is properly correcting star-like objects, but over- or under-correcting the effect in galaxy images. This could be due to a violation of some of the underlying assumptions in the method, including the small-pixel approximation or the linearity of kernel correction. Some of the remaining trends could be related to an increase in signal-to-noise in the images, however this is a universally applicable issue with shape measurement and is beyond the scope of this notebook. In some of the figures, a rapid increase in size can be seen at the highest fluxes, indicating saturation of pixel wells which is unrelated to the brighter-fatter effect." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now plot the flux lost/gained in the process of brighter-fatter correction, by subtracting the flux of the corrected images from the uncorrected ones. The flux measurement is the same as the one used in the above figures and is measured in electrons." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "colorpalette=cycle(plt.cm.viridis(np.linspace(0,1,len(indexfoo))))\n", + "stylepalette=cycle(['s','*','o'])\n", + "plt.figure(figsize=(8,5))\n", + "for nfoo in indexfoo:\n", + " flux_foo=bf_dat['flux'][:,nfoo]\n", + " fluxdiffpct_foo=(bf_dat['flux'][:,nfoo]-bf_dat['corrflux'][:,nfoo])/bf_dat['flux'][:,nfoo]*100.\n", + " plt.plot(flux_foo,fluxdiffpct_foo,label=str(nfoo).zfill(3),c=next(colorpalette),marker=next(stylepalette))\n", + "plt.xscale('log')\n", + "plt.legend()\n", + "plt.xlabel(flux_param,fontsize=20)\n", + "plt.ylabel('Flux change of correction \\n (before - after) [%]',fontsize=20)\n", + "#plt.savefig(cat_dir+'BF_corr_flux_change.png',dpi=150)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# +TODO other ways of doing matching, catalog stacking" + ] + } + ], + "metadata": { + "anaconda-cloud": {}, + "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/ImageProcessing/README.rst b/ImageProcessing/README.rst index 36bdd516..9cb2e0e0 100644 --- a/ImageProcessing/README.rst +++ b/ImageProcessing/README.rst @@ -24,3 +24,14 @@ This folder contains a set of tutorial notebooks exploring the image processing - `Alex Drlica-Wagner `_ + + * - **BrighterFatterCorrection.ipynb** + - Analysis of Beam Simulator Images and Brighter-fatter Correction. + - `ipynb `_, + `rendered `_ + + .. image:: https://github.com/LSSTScienceCollaborations/StackClub/blob/rendered/ImageProcessing/log/BrighterFatterCorrection.svg + :target: https://github.com/LSSTScienceCollaborations/StackClub/blob/rendered/ImageProcessing/log/BrighterFatterCorrection.log + + - `Andrew Bradshaw `_ + diff --git a/Meetings.md b/Meetings.md index effd7653..6bcb98b6 100644 --- a/Meetings.md +++ b/Meetings.md @@ -1,16 +1,31 @@ # Stack Club Meetings -Inidividual session recordings are linked below. +Individual session links and recordings are given below, most recent meeting at the top. | Session | Date | Topic | Links | |---|---|---|---| -| Phase 1, Session 1 | Friday May 25, 2018 [(video)](https://stanford.zoom.us/recording/share/xA33Pv0oq_g5l6a0CaJ0az01mbROy_gyGLDEqIR92FOwIumekTziMw) | Visualization with Firefly | [SQuaRE notebook](https://github.com/lsst-sqre/notebook-demo/blob/master/Firefly.ipynb) | -| Phase 1, Session 2 | Wednesday June 6, 2018 [(video)](https://stanford.zoom.us/recording/share/YZad6BLPZFCjhgLSckrpis7w6Ekyr61VhhIvtFnjR_-wIumekTziMw) | Commissioning Team Bootcamp Report | [LSST Commissioning Team Notebooks](https://github.com/lsst-com/notebooks) | -| Phase 1, Session 3 | Friday June 22, 2018 | | | -| Phase 1, Session 4 | Friday July 6, 2018 [(video)](https://stanford.zoom.us/recording/share/1ZHCNdwRZnhwq8sb1TPvznug-AusUCCIV55N0DUF-LawIumekTziMw) | Project Discussion | [Topic List](https://docs.google.com/document/d/1PSA1uWwTfs9CweatpxF8CEPGBYRY5ZaXB39JzXYE7_U/edit#heading=h.txq6h6bpxzkd) | -| Phase 1, Session 5 | Friday July 13, 2018 [(video)](https://stanford.zoom.us/recording/share/QT8r75yuXR1sjZVkHh4MstfBLJ80wKubuqSvW4s3gfGwIumekTziMw) | Hack Session | [Project List](https://github.com/LSSTScienceCollaborations/StackClub/issues?q=is%3Aopen+label%3Aproject+sort%3Aupdated-desc) | -| Phase 1, Session 6 | Friday July 20, 2018 [(video)](https://stanford.zoom.us/recording/share/XqFx95GJ7zlSZTOVBPZz4l8WGYUmj7EyNsuF6vofMtewIumekTziMw) | Hack Session | [CalExp Tour](https://github.com/LSSTScienceCollaborations/StackClub/blob/project/calexp-tour/stargaser/Basics/Calexp_guided_tour.ipynb) | -| Phase 1, Session 7 | Friday July 27, 2018 [(video)](https://stanford.zoom.us/recording/share/AOFd8Q8yH4lHI6aTLylqRBgcusMUERz-ksiULX4rRL2wIumekTziMw) | Hack Session | [VPN set-up change](https://github.com/LSSTScienceCollaborations/StackClub/blob/master/GettingStarted/GettingStarted.md#accessing-ncsa-via-its-vpn) | +| Phase 2, Session 6 | Friday September 21, 2018 [(video)](https://stanford.zoom.us/recording/share/8_fQYpnZFh2jDLE4LHKjEfdiQ28kjFxGu5jSdTuzdE2wIumekTziMw) | Hacking, Syllabus discussion | ["Course" topic list](https://docs.google.com/document/d/1PSA1uWwTfs9CweatpxF8CEPGBYRY5ZaXB39JzXYE7_U/edit?ts=5ba52b5e#heading=h.txq6h6bpxzkd) | +| Phase 2, Session 5 | Friday September 14, 2018 [(video)](https://stanford.zoom.us/recording/share/-IiuluXvCcOdD-L8FNQSmnB29-f8lU2pTfPyyahcJ1uwIumekTziMw) | Tutorial walkthrough, Hacking | [HSC Re-Run Script and Notebook](https://github.com/LSSTScienceCollaborations/StackClub/blob/project/hsc-re-run/ImageProcessing/Re-RunHSC.ipynb) | +| Phase 2, Session 4 | Friday September 7, 2018 [(video)](https://stanford.zoom.us/recording/share/ZlkFudy5hMTeR-GZVOgo_oGd0R9Q4dkrN6-aJMfelGawIumekTziMw) | Notebook walkthrough, Hacking | [Guided Tour of an AFW Table](https://github.com/LSSTScienceCollaborations/StackClub/blob/project/afw_table/ishasan/Basics/afw_table_guided_tour.ipynb) | +| Phase 2, Session 3 | Friday August 31, 2018 [(video)](https://stanford.zoom.us/recording/share/U7_XJvwjNlUh4N7g3ytBbKtTQHl-fLS0tqiBhAxZrEmwIumekTziMw) | Live code review, hacking | [Brighter-Fatter Correction with Beamsim Data](https://github.com/LSSTScienceCollaborations/StackClub/blob/project/beamsim/andrewkbradshaw/ImageProcessing/BrighterFatterCorrection.ipynb) | +| Phase 2, Session 2 | Friday August 24, 2018 [(video)](https://stanford.zoom.us/recording/share/share/Xii8Utw9RX5rqGUn8a_barg6NDBcRuzkmDIjDrUds82wIumekTziMw) | Interactive visualization, new member start-ups, hacking | [Bokeh/HoloViews Demo](https://github.com/LSSTScienceCollaborations/StackClub/blob/project/bokeh_holoviews_datashader/bechtol/Visualization/bokeh_holoviews_datashader.ipynb) | +| LSST2018 Launch | Monday August 13, 2018 [(video)](https://stanford.zoom.us/recording/share/qyunKljpUWaFQBneuiL4PnbxTB-tf1BvttELFVPJHnuwIumekTziMw) | PCW welcome, discussion, hacking | [Introduction to the LSST Stack Club](https://docs.google.com/presentation/d/1LWShGi-YLqWoxPvewkg-JOpb67WKZyI0YQcSFrmAl14/edit#slide=id.p1) | + + +### Phase 1 Sessions + +Before the August 2018 PCW "launch", we met as a small group, putting together our first notebooks to get the Stack Club started. + +| Session | Date | Topic | Links | +|---|---|---|---| +| Phase 1, Session 9 | Friday August 10, 2018 [(video)](https://stanford.zoom.us/recording/share/d5skJMVG1L-XhtyV6xZA6gI0pN552NQjfNuFUMymocCwIumekTziMw) | Hack Session | [Rules, `stackclub` library package, CIT with `beavis-ci`](https://github.com/LSSTScienceCollaborations/StackClub/issues/85) | | Phase 1, Session 8 | Friday August 3, 2018 [(video)](https://stanford.zoom.us/recording/share/Pnin7IjBNCyGrOCKgyXTvRFFbcuE_eG6tN6QWkQtsvmwIumekTziMw) | PCW Planning | [Source Detection: Low Surface Brightness Galaxies](https://github.com/LSSTScienceCollaborations/StackClub/blob/master/SourceDetection/LowSurfaceBrightness.ipynb) | +| Phase 1, Session 7 | Friday July 27, 2018 [(video)](https://stanford.zoom.us/recording/share/AOFd8Q8yH4lHI6aTLylqRBgcusMUERz-ksiULX4rRL2wIumekTziMw) | Hack Session | [VPN set-up change](https://github.com/LSSTScienceCollaborations/StackClub/blob/master/GettingStarted/GettingStarted.md#accessing-ncsa-via-its-vpn) | +| Phase 1, Session 6 | Friday July 20, 2018 [(video)](https://stanford.zoom.us/recording/share/XqFx95GJ7zlSZTOVBPZz4l8WGYUmj7EyNsuF6vofMtewIumekTziMw) | Hack Session | [CalExp Tour](https://github.com/LSSTScienceCollaborations/StackClub/blob/project/calexp-tour/stargaser/Basics/Calexp_guided_tour.ipynb) | +| Phase 1, Session 5 | Friday July 13, 2018 [(video)](https://stanford.zoom.us/recording/share/QT8r75yuXR1sjZVkHh4MstfBLJ80wKubuqSvW4s3gfGwIumekTziMw) | Hack Session | [Project List](https://github.com/LSSTScienceCollaborations/StackClub/issues?q=is%3Aopen+label%3Aproject+sort%3Aupdated-desc) | +| Phase 1, Session 4 | Friday July 6, 2018 [(video)](https://stanford.zoom.us/recording/share/1ZHCNdwRZnhwq8sb1TPvznug-AusUCCIV55N0DUF-LawIumekTziMw) | Project Discussion | [Topic List](https://docs.google.com/document/d/1PSA1uWwTfs9CweatpxF8CEPGBYRY5ZaXB39JzXYE7_U/edit#heading=h.txq6h6bpxzkd) | +| Phase 1, Session 3 | Friday June 22, 2018 | | | +| Phase 1, Session 2 | Wednesday June 6, 2018 [(video)](https://stanford.zoom.us/recording/share/YZad6BLPZFCjhgLSckrpis7w6Ekyr61VhhIvtFnjR_-wIumekTziMw) | Commissioning Team Bootcamp Report | [LSST Commissioning Team Notebooks](https://github.com/lsst-com/notebooks) | +| Phase 1, Session 1 | Friday May 25, 2018 [(video)](https://stanford.zoom.us/recording/share/xA33Pv0oq_g5l6a0CaJ0az01mbROy_gyGLDEqIR92FOwIumekTziMw) | Visualization with Firefly | [SQuaRE Firefly demo](https://github.com/lsst-sqre/notebook-demo/blob/master/Firefly.ipynb) | diff --git a/README.md b/README.md index 9f4ec115..f02b30c6 100644 --- a/README.md +++ b/README.md @@ -17,16 +17,20 @@ the LSST Science Collaborations. | Visualization | Displaying images and catalogs. | [StackClub/Visualization](Visualization) | | Image Processing | From raw images to `calexp`s and `coadd`s. | [StackClub/ImageProcessing](ImageProcessing) | | SourceDetection | Detection of sources in images - including low surface brightness galaxies. | [StackClub/SourceDetection](SourceDetection) | +| Deblending | Deblending the objects | [StackClub/Deblending](Deblending) | | Validation | Tools for validating Stack outputs, example validation analyses | [StackClub/Validation](Validation) | * [Stack Club projects](https://github.com/LSSTScienceCollaborations/StackClub/labels/project), as defined by Stack Club members - follow [this link](https://github.com/LSSTScienceCollaborations/StackClub/labels/project) to see what people are working on. [Unassigned projects](https://github.com/LSSTScienceCollaborations/StackClub/issues?utf8=%E2%9C%93&q=is%3Aopen+label%3Aproject+no%3Aassignee) are available for new members to take on! * [Working list of target topics, with links to tutorial seeds](https://docs.google.com/document/d/1PSA1uWwTfs9CweatpxF8CEPGBYRY5ZaXB39JzXYE7_U/edit#), for help in defining a new Stack Club project. This list is a fairly comprehensive collection of existing project and community tutorial web pages and demo notebooks, from which seeds can be drawn. +## Joining the Stack Club +If you would like to join the Stack Club, please fill out this short **[application form](https://goo.gl/forms/588KlPTFfkEEFFUu2)**. (Basically you'll be asked to agree to abide by the [Stack Club Rules](Rules.md), and then give enough contact information to request an account on the LSST Science Platform.) If you are not ready to commit time to working on a Stack Club project, you can still follow along by [watching](https://github.com/LSSTScienceCollaborations/StackClub/subscription) this repo and joining the [#stack-club LSSTC Slack channel](https://lsstc.slack.com/messages/C9YRAS4HM/). + ## Contributing -New Stack Club members: please see the [notes on getting started](GettingStarted/GettingStarted.md) - they'll walk you onto you new LSST Science Platform account, and then show you how to work on your tutorial notebooks. Also, please note the [Stack Club Rules](Rules.md) that we all agree to abide by. +New Stack Club members: please see the [notes on getting started](GettingStarted/GettingStarted.md) - they'll walk you onto you new LSST Science Platform account, and then show you how to work on your tutorial notebooks. -Everyone else: we welcome pull requests! Feel free to fork this repo and send us a pull request. And if you are interested in joining the Stack Club, please drop one of us a line, or come and find us in the [#stack-club](https://lsstc.slack.com/messages/C9YRAS4HM) LSSTC Slack channel. +Everyone else: we welcome pull requests! Feel free to fork this repo and send us a pull request. And if you are interested in joining the Stack Club, please drop one of us a line, or come and find us in the [#stack-club](https://lsstc.slack.com/messages/C9YRAS4HM) LSSTC Slack channel. > When preparing a pull request, please note the [standards](https://github.com/LSSTScienceCollaborations/StackClub/blob/master/GettingStarted/GettingStarted.md#standards) that we are trying to uphold. @@ -38,7 +42,7 @@ We welcome your input! Please post questions and suggestions in the * Alex Drlica-Wagner (Fermilab, [@kadrlica](https://github.com/LSSTScienceCollaborations/StackClub/issues/new?body=@kadrlica)) * Phil Marshall (SLAC, [@drphilmarshall](https://github.com/LSSTScienceCollaborations/StackClub/issues/new?body=@drphilmarshall)) -The Club meets periodically via Zoom, but you can find us on LSSTC Slack at [#stack-club](https://lsstc.slack.com/messages/C9YRAS4HM). You can also watch the tutorial walkthroughs in the Club sessions in the videos linked from our [Meetings page](Meetings.md). +The Club meets periodically via Zoom, but you can find us on LSSTC Slack at [#stack-club](https://lsstc.slack.com/messages/C9YRAS4HM). You can also watch the tutorial walkthroughs in the Club sessions in the videos linked from our [Meetings page](Meetings.md). If you are just looking for the application form, it's [here](https://goo.gl/forms/588KlPTFfkEEFFUu2). ## License @@ -52,4 +56,4 @@ but you can't blame us if it doesn't do what you want. ## More About This Project -Following a successful LSSTC "Enabling Science" proposal, we put together a 3-phase plan, which you can read about in more detail [here](https://docs.google.com/document/d/103kzjOklSUWo5MJP9B-EsnAdO7V6bstTC_mzBvd0NIk/edit#). Phase 0 involved collecting existing tutorials and identifying potential club members from around the LSST Science Collaborations. Then, in Phase 1 (late May 2018 to mid August 2018) we worked together in a small group to turn a subset of those existing "seed" tutorials into community-maintained Jupyter notebooks, for display at the August LSST 2018 Project and Community Workshop (PCW) in Tucson. At that meeting, we opened up to a larger group of LSST science collaboration members, extending and spinning off the initial set of notebooks. +Following a successful LSSTC "Enabling Science" proposal, we put together a 3-phase plan, which you can read about in more detail [here](https://docs.google.com/document/d/103kzjOklSUWo5MJP9B-EsnAdO7V6bstTC_mzBvd0NIk/edit#). Phase 0 involved collecting existing tutorials and identifying potential club members from around the LSST Science Collaborations. Then, in Phase 1 (late May 2018 to mid August 2018) we worked together in a small group to turn a subset of those existing "seed" tutorials into community-maintained Jupyter notebooks, for display at the August LSST 2018 Project and Community Workshop (PCW) in Tucson. At that meeting, we opened up to a larger group of LSST science collaboration members, extending and spinning off the initial set of notebooks. diff --git a/SourceDetection/LowSurfaceBrightness.ipynb b/SourceDetection/LowSurfaceBrightness.ipynb index 50797916..babc6873 100644 --- a/SourceDetection/LowSurfaceBrightness.ipynb +++ b/SourceDetection/LowSurfaceBrightness.ipynb @@ -7,11 +7,11 @@ "# Low-Surface Brightness Source Detection\n", "\n", "
Owner: **Alex Drlica-Wagner** ([@kadrlica](https://github.com/LSSTScienceCollaborations/StackClub/issues/new?body=@kadrlica))\n", - "
Last Verified to Run: **2018-07-22**\n", - "
Verified Stack Release: **w201829**\n", + "
Last Verified to Run: **2018-08-13**\n", + "
Verified Stack Release: **v16.0**\n", "\n", - "This Notebook demonstrates how to run the source detection, measurment, and deblending algorithms with a focus on optimizing for low-surface brightness object detection. It attempts to split out the source detection and measurement algorithms from `processCCD` and apply them to the search for low surface brightness galaxies. The content of this notebook builds off of Robert Lupton's [Greco LSB.ipynb](https://github.com/RobertLuptonTheGood/notebooks/blob/master/Demos/Greco%20LSB.ipynb) with some source detection and measurement details from [Tune Detection.ipynb](https://github.com/RobertLuptonTheGood/notebooks/blob/master/Demos/Tune%20Detection.ipynb) and [Kron.ipynb](https://github.com/RobertLuptonTheGood/notebooks/blob/master/Demos/Kron.ipynb).\n", - "Interaction with `lsst.afw.display` was also improved by studying Michael Wood-Vasey's [DC2_Postage Stamps.ipynb](https://github.com/LSSTDESC/DC2_Repo/blob/master/Notebooks/DC2_Postage_Stamps.ipynb).\n", + "This Notebook demonstrates how to run the source detection, measurment, and deblending algorithms with a focus on optimizing for low-surface brightness object detection. It attempts to split out the source detection and measurement algorithms from `processCCD` and apply them to the search for low surface brightness galaxies. The content of this notebook builds off of an analysis from Johnny Greco, adapted into notebook form in Robert Lupton's [Greco LSB.ipynb](https://github.com/RobertLuptonTheGood/notebooks/blob/master/Demos/Greco%20LSB.ipynb). Some source detection and measurement details come from [Tune Detection.ipynb](https://github.com/RobertLuptonTheGood/notebooks/blob/master/Demos/Tune%20Detection.ipynb) and [Kron.ipynb](https://github.com/RobertLuptonTheGood/notebooks/blob/master/Demos/Kron.ipynb).\n", + "Interaction with `lsst.afw.display` was also improved by studying Michael Wood-Vasey's [DC2_Postage Stamps.ipynb](https://github.com/LSSTDESC/DC2-analysis/blob/master/tutorials/dm_butler_postage_stamps.ipynb).\n", "\n", "### Learning Objectives:\n", "After working through and studying this notebook you should be able to\n", @@ -111,7 +111,7 @@ "source": [ "### Data access\n", "\n", - "Here we use the `butler` to access a `calexp` from the Twinkles data subset. More information on the `butler` can be found in [Butler_Tutorial.ipynb](), while a deeper examination of the `calexp` object can be found in [Calexp_Tutorial.ipynb](). We expect the user to have a working knowledge of these objects." + "Here we use the `butler` to access a `calexp` from the Twinkles data subset. More information on the `butler` will be available in `Butler_Tutorial.ipynb` (TBD), while a deeper examination of the `calexp` object can be found in [Calexp_guided_tour.ipynb](https://github.com/LSSTScienceCollaborations/StackClub/blob/master/Basics/Calexp_guided_tour.ipynb). We expect the user to have a working knowledge of these objects." ] }, { @@ -449,7 +449,7 @@ "source": [ "if False:\n", " sources.writeFits(\"outputTable.fits\")\n", - " exposure.writeFits(\"example1-out.fits\")" + " calexp.writeFits(\"example1-out.fits\")" ] }, { @@ -629,6 +629,13 @@ "plt.gca().axis('off')" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The procedure above masks all `DETECTED` pixels with a simplistic selection. However, the Stack provides an alternative mechanism for more directed masking through the [meas.base.noiseReplacer](http://doxygen.lsst.codes/stack/doxygen/x_masterDoxyDoc/classlsst_1_1meas_1_1base_1_1noise_replacer_1_1_noise_replacer.html). The `noiseReplacer` takes as input the calexp object and the footprint set (`fpset`) returned `sourceDetectionTask`. With the `noiseReplacer`, it is possible to selectively replace a subset of the sources." + ] + }, { "cell_type": "markdown", "metadata": {}, diff --git a/Visualization/AFW_Display_Demo.ipynb b/Visualization/AFW_Display_Demo.ipynb new file mode 100644 index 00000000..8a7b7ae1 --- /dev/null +++ b/Visualization/AFW_Display_Demo.ipynb @@ -0,0 +1,376 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "# **Demo of lsst.afw.display -- displaying images using the LSST DM Astronomical Framework library**\n", + "\n", + "**Owner:** Brant Robertson ([@brantr](https://github.com/LSSTScienceCollaborations/StackClub/issues/new?body=@brantr)) \n", + "**Level:** Introductory \n", + "**Last Verified to Run:** 2018-08-24 \n", + "**Verified Stack Release:** v16.0 \n", + "\n", + "## **Learning Objectives:**\n", + "\n", + "In this tutorial we will \n", + "\n", + "* Show how to access the `lsst.afw.display` routines.\n", + "\n", + "* Use the LSST data Butler to access processed data and inspect it visually.\n", + "\n", + "This tutorial is designed to help users get a brief feel for the `lsst.afw.display` library that enables the visual inspection of data. The [`lsst.afw` library](https://github.com/lsst/afw) provides an \"Astronomical Framework\" (afw) while the `lsst.daf.*` libraries (see, e.g., [daf_base](https://github.com/lsst/daf_base)) provides a Data Access Framework (daf). Both libraries are used in this tutorial, with the `lsst.daf.persistence` library used to access a calibrated exposure (calexp) and the `lsst.afw.display` library used to show the exposure image on the screen.\n", + "\n", + "This tutorial made use of the [`LowSurfaceBrightness.ipynb` StackClub notebook](https://nbviewer.jupyter.org/github/LSSTScienceCollaborations/StackClub/blob/rendered/SourceDetection/LowSurfaceBrightness.nbconvert.ipynb) by [Alex Drlica-Wagner](https://github.com/LSSTScienceCollaborations/StackClub/issues/new?body=@kadrlica). More examples of the use of `lsst.afw.display` can be found in the [Stack ](https://pipelines.lsst.io/getting-started/display.html)." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "## **Step 0) Import Common Python Libraries**\n", + "\n", + "The [`matplotlib`](https://matplotlib.org/), [`numpy`](http://www.numpy.org/), and [`astropy`](http://www.astropy.org/) libraries are widely used Python libraries for plotting, scientific computing, and astronomical data analysis. We will use these packages in common ways below, including the `matplotlib.pyplot` plotting sublibrary. We also import the [`warnings` library](https://docs.python.org/2/library/warnings.html) to prevent some routine warning messages from printing to the screen." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#allow for matplotlib to create inline plots in our notebook\n", + "%matplotlib inline \n", + "import numpy as np #imports numpy with the alias np\n", + "import matplotlib.pyplot as plt #imports matplotlib.pyplot as plt\n", + "import warnings #imports the warnings library" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's go ahead and import from `astropy` the image stretch limits from the familiar `zscale()` function." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from astropy.visualization import ZScaleInterval #This function allows use to use the `zscale()` rescaling limits function familiar from, e.g., DS9, to adjust the image stretch.\n", + "zscale = ZScaleInterval() #create an alias to the `ZScaleInterval()` function" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "And let the kernel know that we're happy not to have some useful warnings printed during this tutorial." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "warnings.simplefilter(\"ignore\", category=FutureWarning) #prevent some helpful but ancillary warning messages from printing during some LSST DM Release calls\n", + "warnings.simplefilter(\"ignore\", category=UserWarning) #prevent some helpful but ancillary warning messages from printing during some LSST DM Release calls" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As a last preparatory task, we set the parameters of `matplotlib.pyplot` to give us a large default size for an image." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.rcParams['figure.figsize'] = (8.0, 8.0) #set a large default size for our images" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "## **Step 1) Loading the LSST DM Stack**\n", + "\n", + "To manipulate data, the LSST DM Stack provides a `Butler` that enables generic access routines to DM-generated data. For more information, see [the Data Butler entry in the LSST Software User Guide](https://confluence.lsstcorp.org/display/LSWUG/Data+Butler). In order to access a calibrated exposure from data stored in the format required by the LSST Data Butler, we must load the `lsst.daf.persistence` library to produce a Butler instance from the data." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import lsst.daf.persistence as dafPersist #load lsst.daf.persistence to gain access to a Butler instance" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next, we need to load the `lsst.afw.display` library to gain access to the image visualization routines we'd like to use." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import lsst.afw.display as afwDisplay #load lsst.afw.display to gain access to image visualization routines." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## **Step 2) Importing Data to Visualize**\n", + "\n", + "To plot an image to the screen, we must first load some data. In this tutorial, we will use the `Twinkles` simulated images available in the StackClub data repository. These data sit in the data directory `/project/shared/data/Twinkles_subset/output_data_v2` and contain a set of data produced in generating a calibrated exposure by the DM Stack. These data are organized in a structure that enables a DM Stack `Butler` instance to be generated and provide access to a single filter image (in this case `r` band), a specific detector raft (2,2), a specific sensor in the raft (`1,1`) and a specific visit (in this case, 235 -- note only one band is available per visit in this example).\n", + "\n", + "Once we define a string that contains the data directory, we start the `Butler` instance using the `lsst.daf.persistence` library alias `dafPersist` and its `Butler` class. The `Butler` object is initialized with a string containing the data directory we wish to access. Running the cell may take a few moments." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "datadir = \"/project/shared/data/Twinkles_subset/output_data_v2\" #our data directory containing the Twinkles data organized as Butler expects\n", + "butler = dafPersist.Butler(datadir) #create an instance of the Butler, which we call `butler`, with access to our data" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "With the `Butler` instance now generated using our data directory, we can retrieve the desired calibrated exposure by telling the butler which filter, raft, sensor, and visit we wish to view. To do this, we define dictionary with the required information." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Grab a calexp of interest\n", + "dataId = {'filter': 'r', 'raft': '2,2', 'sensor': '1,1', 'visit': 235} #Define a dictionary with the filter, raft, sensor, and visit we wish to view\n", + "calexp = butler.get('calexp', **dataId) #retrieve the data using the `butler` instance and its function `get()`" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## **Step 3.1) Use AFWDisplay to Visualize the Image**\n", + "\n", + "Now, with a `Butler` instance defined and a calibrated exposure retrieved, we can use [`lsst.afw.display`](https://github.com/lsst/afw/tree/master/python/lsst/afw/display) to visualize the data. The next task is to let AFWDisplay know that we want it to enroll `matplotlib` as our default display backend. To do this, we use the `setDefaultBackend()` function. Remember that we made an alias to `lsst.afw.display` called `afwDisplay`, so we'll use that to call `setDefaultBackend()`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "afwDisplay.setDefaultBackend('matplotlib') # Use lsst.afw.display with the matplotlib backend" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We are now set to display the image. To do this, we:\n", + "\n", + "* First create a `matplotlib.pyplot` figure using `plt.figure()` -- this will be familiar to anyone with experience using `matplotlib`.\n", + "* Then create an alias to the `lsst.afw.display.Display` method that will allow us to display the data to the screen. This alias will be called `afw_display`.\n", + "* Before showing the data on the screen, we have to decide how to apply an image stretch given the data. The algorithm we'll use is `asinh` familiar from SDSS images, with a range of values set by `zscale`. To do this, we use the `scale()` function provided by `lsst.afw.display`. See the `scale()` function definition in the [`interface.py` file of the lsst.afw.display library](https://github.com/lsst/afw/blob/master/python/lsst/afw/display/interface.py).\n", + "* Finally, we can display the image. Do do this, we provide the `mtv()` method the `image` member of our calibrated image retrieved by the `butler`. We can then use `plt.show()` to display our figure.\n", + "\n", + "All these tasks are best done within the same notebook cell." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure() #create a matplotlib.pyplot figure\n", + "afw_display = afwDisplay.Display() #get an alias to the lsst.afw.display.Display() method\n", + "afw_display.scale('asinh', 'zscale') #set the image stretch algorithm and range\n", + "afw_display.mtv(calexp.image) #load the image into the display\n", + "plt.show() #show the corresponding pyplot figure" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**Congrats!** We've plotted an image using `lsst.afw.display`!" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## **Step 3.2) Use AFWDisplay to Visualize the Image and Mask Plane**\n", + "\n", + "The `calexp` returned by the butler contains more than just the image pixel values (see the [calexp tutorial](https://github.com/LSSTScienceCollaborations/StackClub/blob/master/Basics/Calexp_guided_tour.ipynb) for more details). One other component is the mask plane associated with the image. `AFWDisplay` provides a nice pre-packaged interface for overplotting the mask associated with an image. A mask is composed of a set of \"mask planes\", 2D binary bit maps corresponding to pixels that are masked for various reasons (see [here](https://pipelines.lsst.io/v/DM-11392/getting-started/display.html#interpreting-displayed-mask-colors) for more details)." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We'll follow the same steps as above to display the image, but we'll add a few modifications\n", + "\n", + "* We explicitly set the transparency of the overplotted mask (0 = transparent, 1 = opaque)\n", + "* We explicitly set the color of the 'DETECTED' mask plane to 'blue' (i.e. all pixels associated with detected objects).\n", + "* We pass the full `calexp` object to `mtv` instead of just the image plane." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure() #create a matplotlib.pyplot figure\n", + "afw_display = afwDisplay.Display() #get an alias to the lsst.afw.display.Display() method\n", + "afw_display.scale('asinh', 'zscale') #set the image stretch algorithm and range\n", + "afw_display.setMaskTransparency(0.4) #set the transparency of the mask plane (1 = opaque)\n", + "afw_display.setMaskPlaneColor('DETECTED','blue') #set the color for a single plane in the mask\n", + "afw_display.mtv(calexp) #load the image and mask plane into the display\n", + "plt.show() #show the corresponding pyplot figure" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The `afw_display` object contains more information about the mask planes that can be accessed" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(\"Mask plane bit definitions:\\n\", afw_display.getMaskPlaneColor()) # Print the colors associated to each plane in the mask\n", + "print(\"\\nMask plane methods:\\n\")\n", + "help(afw_display.setMaskPlaneColor)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## **Step 4) More Information about lsst.afw.display**\n", + "\n", + "To get some more information about `lsst.afw.display`, we can print the method list to see what's available. The next cell will print `lsst.afw.display` methods to the screen." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "method_list = [func for func in dir(afw_display) if callable(getattr(afw_display, func))]\n", + "print(method_list)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "If you'd like to learn more about any given function, please see the [`lsst.afw.display` source code](https://github.com/lsst/afw/tree/master/python/lsst/afw/display).\n", + "\n", + "You can also read the API documentation about the above functions using the Jupyter notebook `help()` function:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "help(afw_display.scale)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "help(afw_display.mtv)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### **Further Documentation**\n", + "\n", + "If you'd like some more information on `lsst.afw.display`, please have a look at the following websites:\n", + "\n", + "* [Info on image indexing conventions.](https://github.com/lsst/afw/blob/master/doc/lsst.afw.image/indexing-conventions.rst) \n", + "* [afw.display Doxygen website](http://doxygen.lsst.codes/stack/doxygen/x_masterDoxyDoc/namespacelsst_1_1afw_1_1display.html) \n", + "* [afw.display GitHub website](https://github.com/RobertLuptonTheGood/afw/tree/master/python/lsst/afw/display) \n", + "* [The `pipelines.lsst.io` Getting Started on Image Display website.](https://pipelines.lsst.io/getting-started/display.html)" + ] + } + ], + "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 +} diff --git a/Visualization/README.rst b/Visualization/README.rst index 949fff7e..58928d29 100644 --- a/Visualization/README.rst +++ b/Visualization/README.rst @@ -14,12 +14,23 @@ See the index table below for links to the notebook code, and an auto-rendered v - Owner - * - **bokeh_holoviews_datashader.ipynb** - - Examples of interactive visualization with the Boken, HoloViews, and Datashader plotting packages available in PyViz suite of data analysis python modules; brushing and linking with large datasets + * - **Firefly Visualization Demo** + - Introduction to the Firefly interactive plotter and image viewer. + - `ipynb `_, `video `_ + - `Simon Krughoff `_ + + + * - **Interactive Visualization with Bokeh, HoloViews, and Datashader** + - Examples of interactive visualization with the Boken, HoloViews, and Datashader plotting packages available in PyViz suite of data analysis python modules; brushing and linking with large datasets - `ipynb `_, `rendered `_ .. image:: https://github.com/LSSTScienceCollaborations/StackClub/blob/rendered/Visualization/log/bokeh_holoviews_datashader.svg :target: https://github.com/LSSTScienceCollaborations/StackClub/blob/rendered/Visualization/log/bokeh_holoviews_datashader.log - - `Keith Bechtol `_ \ No newline at end of file + - `Keith Bechtol `_ + + * - **"With Globular" LSST 2018 Tutorial** + - General purpose tutorial including interactive Firefly visualization. + - `ipynb `_ + - `Jim Bosch `_ diff --git a/Visualization/bokeh_holoviews_datashader.ipynb b/Visualization/bokeh_holoviews_datashader.ipynb index 33993b00..53de322d 100644 --- a/Visualization/bokeh_holoviews_datashader.ipynb +++ b/Visualization/bokeh_holoviews_datashader.ipynb @@ -154,8 +154,10 @@ "metadata": {}, "outputs": [], "source": [ + "width = 300\n", + "\n", "# create a new plot and add a renderer\n", - "left = figure(tools=TOOLS_LEFT, plot_width=500, plot_height=500, output_backend=\"webgl\",\n", + "left = figure(tools=TOOLS_LEFT, plot_width=width, plot_height=width, output_backend=\"webgl\",\n", " title='Spatial: Centered on (RA, Dec) = (%.2f, %.2f)'%(ra_target, dec_target))\n", "left.circle('x0', 'y0', hover_color='firebrick', source=source,\n", " selection_fill_color='steelblue', selection_line_color='steelblue',\n", @@ -166,7 +168,7 @@ "left.yaxis.axis_label = 'Delta DEC'\n", "\n", "# create another new plot and add a renderer\n", - "right = figure(tools=TOOLS_RIGHT, plot_width=500, plot_height=500, output_backend=\"webgl\",\n", + "right = figure(tools=TOOLS_RIGHT, plot_width=width, plot_height=width, output_backend=\"webgl\",\n", " title='CMD')\n", "right.circle('x1', 'y1', hover_color='firebrick', source=source,\n", " selection_fill_color='steelblue', selection_line_color='steelblue',\n", @@ -237,6 +239,40 @@ "print(selection.index)" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "selected_points = points[selection.values]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "len(selected_points)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "selected_points" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "> For more help on selecting points in HoloViews, see the [user guide](http://build.holoviews.org/User_Guide/Indexing_and_Selecting_Data.html)." + ] + }, { "cell_type": "markdown", "metadata": {},