diff --git a/Basics/DataInventory.ipynb b/Basics/DataInventory.ipynb deleted file mode 100644 index 394a9a12..00000000 --- a/Basics/DataInventory.ipynb +++ /dev/null @@ -1,334 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": { - "slideshow": { - "slide_type": "slide" - } - }, - "source": [ - "# Exploring the Shared Datasets in the LSST Science Platform\n", - "
Owner(s): **Phil Marshall** ([@drphilmarshall](https://github.com/LSSTScienceCollaborations/StackClub/issues/new?body=@drphilmarshall)), \n", - "
Last Verified to Run: **2018-08-05**\n", - "
Verified Stack Release: **16.0**\n", - "\n", - "In this notebook we'll take a look at some of the datasets available on the LSST Science Platform. \n", - "\n", - "### Learning Objectives:\n", - "\n", - "After working through this tutorial you should be able to: \n", - "1. Start figuring out which of the available datasets is going to be of most use to you in any given project; \n", - "\n", - "When it is finished, you should be able to:\n", - "2. Plot the patches and tracts in a given dataset on the sky;\n", - "3. List the available catalogs in a given dataset.\n", - "\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": [ - "We'll need the `stackclub` package to be installed. If you are not developing this package, you can install it using `pip`, like this:\n", - "```\n", - "pip install git+git://github.com/LSSTScienceCollaborations/StackClub.git#egg=stackclub\n", - "```\n", - "If you are developing the `stackclub` package (eg by adding modules to it to support the Stack Club tutorial that you are writing, you'll need to make a local, editable installation. In the top level folder of the `StackClub` repo, do:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "! cd .. && python setup.py -q develop --user && cd -" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "When editing the `stackclub` package files, we want the latest version to be imported when we re-run the import command. To enable this, we need the %autoreload magic command." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "%load_ext autoreload\n", - "%autoreload 2" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "slideshow": { - "slide_type": "subslide" - } - }, - "source": [ - "For accessing the datasets using the Butler, and then visualizing the results, we'll need the following modules:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "%matplotlib inline\n", - "# %matplotlib ipympl\n", - "\n", - "import os, glob\n", - "import numpy as np\n", - "import matplotlib as mpl\n", - "import matplotlib.pyplot as plt\n", - "from IPython.display import IFrame, display, Markdown\n", - "\n", - "import stackclub" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import lsst.daf.persistence as dafPersist\n", - "import lsst.daf.base as dafBase\n", - "\n", - "import lsst.afw.math as afwMath\n", - "import lsst.afw.geom as afwGeom\n", - "\n", - "import lsst.afw.detection as afwDetect\n", - "import lsst.afw.image as afwImage\n", - "import lsst.afw.table as afwTable\n", - "\n", - "import lsst.afw.display as afwDisplay" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "You can find the Stack version that this notebook is running by using eups list -s on the terminal command line:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# What version of the Stack am I using?\n", - "! echo $HOSTNAME\n", - "! eups list -s | grep lsst_distrib" - ] - }, - { - "cell_type": "markdown", - "metadata": { - "slideshow": { - "slide_type": "slide" - } - }, - "source": [ - "## Listing the Available Datasets\n", - "First, let's look at what is currently available. There are two primary shared dataset folders in the LSP, the read-only `/datasets` folder, and the group-writeable folder `/projects/shared/datasets`. Let's see what's in there" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "**`/projects/shared/data`:** These datasets are designed to be small test sets, ideal for tutorials." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "shared_datasets = ! ls -d /project/shared/data/* | grep -v README\n", - "shared_datasets" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "%%bash\n", - "shared_datasets=$( ls -d /project/shared/data/* | grep -v README )\n", - "for dataset in $shared_datasets; do\n", - " du -sh $dataset\n", - "done" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "**`/datasets`:**\n", - "These are typically much bigger: to measure the size, uncomment the second cell below and edit it to target the dataset you are interested in. Running `du` on all folders takes several minutes." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "datasets = ! ls -d /datasets/* | grep -v USAGE | grep -v html\n", - "datasets" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# %%bash\n", - "# datasets=$( ls -d /datasets/* | grep -v USAGE | grep -v html )\n", - "# for dataset in $datasets; do\n", - "# du -h $dataset\n", - "# done" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Visualizing Sky Coverage\n", - "In this section, we'll plot the available patches and tracts in a given dataset on the sky, following the LSST DESC tutorial [dm_butler_skymap.ipynb](https://github.com/LSSTDESC/DC2-analysis/blob/master/tutorials/dm_butler_skymap.ipynb). In fact, we will _import_ this notebook, so that we can re-use its functions. This operation is handled by the `stackclub.wimport` function." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "dm_butler_skymap_notebook = \"https://github.com/LSSTDESC/DC2-analysis/raw/master/tutorials/dm_butler_skymap.ipynb\"\n", - "\n", - "skymapper = stackclub.wimport(dm_butler_skymap_notebook, vb=True)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "> BUG: remote notebooks are not yet `wimport`-able. A workaround could be to import the downloaded file explicitly. This is not yet working, hence the commented out failed attempt below." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# import sys, os\n", - "# import stackclub\n", - "# sys.path.append(os.getcwd() + '/.downloads')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# import dm_butler_skymap" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Now we can attempt to plot the available tracts, using the `plot_skymap_tract()` function." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# repo = \"/project/shared/data/Twinkles_subset/output_data_v2\"\n", - "repo = \"/datasets/hsc/repo/rerun/DM-13666/WIDE\"\n", - "butler = dafPersist.Butler(repo)\n", - "\n", - "# Glob the merged coadd folder for the tracts that have data. Unfortunately, this information is not\n", - "# directly accessible from the data butler.\n", - "tracts = sorted([int(os.path.basename(x)) for x in\n", - " glob.glob(os.path.join(repo, 'deepCoadd-results', 'merged', '*'))])\n", - "\n", - "# How many tracts do we have?\n", - "print(\"Found {} tracts\".format(len(tracts)))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "\"\"\"\n", - "Uncomment this cell when the `wimport` bug is fixed (or avoided).\n", - "\n", - "# Now, loop over all the tracts, plotting them as gray, numbered, rectangles:\n", - "ax = None\n", - "for tract in tracts:\n", - " skyMap = butler.get('deepCoadd_skyMap')\n", - " ax = skymapper.plot_skymap_tract(skyMap, tract=tract, title='', ax=ax)\n", - "\"\"\";" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Summary\n", - "\n", - "In this notebook we took a first look at the datasets available to us in two shared directories in the LSST science platform filesystem." - ] - } - ], - "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/ImageProcessing/Re-RunHSC.ipynb b/ImageProcessing/Re-RunHSC.ipynb index cd8cbacb..3d53533d 100644 --- a/ImageProcessing/Re-RunHSC.ipynb +++ b/ImageProcessing/Re-RunHSC.ipynb @@ -7,7 +7,8 @@ "# HSC Re-Run: Making Forced Photometry Light Curves from Scratch\n", "\n", "
Owner: **Justin Myles** ([@jtmyles](https://github.com/LSSTScienceCollaborations/StackClub/issues/new?body=@jtmyles))\n", - "
Last Verified to Run: **N/A -- in development**\n", + "
Last Verified to Run: **2018-08-10**\n", + "
Verified Stack Release: **16.0**\n", "\n", "This project addresses issue [#63: HSC Re-run](https://github.com/LSSTScienceCollaborations/StackClub/issues/63)\n", "\n", @@ -32,12 +33,15 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import os\n", + "os.system(\"eups list lsst_distrib\") #todo\n", "import sys\n", + "import matplotlib.pyplot as plt\n", + "%matplotlib inline\n", "import eups.setupcmd" ] }, @@ -63,14 +67,15 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 2, "metadata": {}, "outputs": [], "source": [ - "!eups list lsst_distrib\n", - "datarepo = \"/home/jmyles/repositories/ci_hsc/\"\n", - "datadir = \"/home/jmyles/DATA/\"\n", - "os.system(\"mkdir -p {}\".format(datadir))" + "USER = \"jmyles\"\n", + "# todo: $USER\n", + "datarepo = \"/home/{}/repositories/ci_hsc/\".format(USER)\n", + "datadir = \"/home/{}/DATA/\".format(USER)\n", + "os.system(\"mkdir -p {}\".format(datadir));" ] }, { @@ -110,6 +115,7 @@ "outputs": [], "source": [ "# ingest script\n", + "# todo: $USER ()\n", "!ingestImages.py /home/jmyles/DATA /home/jmyles/repositories/ci_hsc/raw/*.fits --mode=link" ] }, @@ -171,12 +177,10 @@ "outputs": [], "source": [ "\"\"\"\n", - "# review what data will be processed\n", - "!processCcd.py DATA --rerun processCcdOutputs --id --show data\n", - "# id allows you to select data by data ID\n", - "# unspecified id selects all raw data\n", - "# example IDs: raw, filter, visit, ccd, field\n", - "# show data turns on dry-run mode\n", + "!processCcd.py DATA --rerun processCcdOutputs --id\n", + "# all cl tasks write output datasets to a Butler repo\n", + "# --rerun configured to write to processCcdOutputs\n", + "# other option is --output\n", "\"\"\"" ] }, @@ -186,7 +190,52 @@ "metadata": {}, "outputs": [], "source": [ - "#!which processCcd.py" + "!which processCcd.py\n", + "\"\"\"\n", + "processCcd.py\n", + "from lsst.pipe.tasks.processCcd import ProcessCcdTask\n", + "\n", + "ProcessCcdTask.parseAndRun()\n", + "\"\"\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# show source of lsst.pipe.tasks.processCcd\n", + "# emacs /opt/lsst/software/stack/stack/miniconda3-4.3.21-10a4fa6/Linux64/pipe_tasks/16.0+1/python/lsst/pipe/tasks/processCcd.py\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from lsst.pipe.tasks.processCcd import ProcessCcdTask, ProcessCcdConfig\n", + "processCcdTaskInstance = ProcessCcdTask(butler=butler)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from stackclub import where_is\n", + "where_is(processCcdTaskInstance, in_the=\"source\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "processCcdConfig = ProcessCcdConfig()" ] }, { @@ -196,13 +245,24 @@ "outputs": [], "source": [ "\"\"\"\n", - "!processCcd.py DATA --rerun processCcdOutputs --id\n", - "# all cl tasks write output datasets to a Butler repo\n", - "# --rerun configured to write to processCcdOutputs\n", - "# other option is --output\n", + "# review what data will be processed\n", + "!processCcd.py DATA --rerun processCcdOutputs --id --show data\n", + "# id allows you to select data by data ID\n", + "# unspecified id selects all raw data\n", + "# example IDs: raw, filter, visit, ccd, field\n", + "# show data turns on dry-run mode\n", "\"\"\"" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#!which processCcd.py" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -245,6 +305,117 @@ "# - TAN: tangent-plane projection\n", "\"\"\"" ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Part 6: Multi-band catalog analysis\n", + "https://pipelines.lsst.io/getting-started/multiband-analysis.html" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import lsst.daf.persistence as dafPersist\n", + "datadir = \"/home/{}/DATA/\".format(os.environ['USER'])\n", + "butler = dafPersist.Butler(inputs=datadir + 'rerun/coaddForcedPhot/')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Access the sources identified from the coadd images" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "rSources = butler.get('deepCoadd_forced_src', {'filter': 'HSC-R', 'tract': 0, 'patch': '1,1'})\n", + "iSources = butler.get('deepCoadd_forced_src', {'filter': 'HSC-I', 'tract': 0, 'patch': '1,1'})" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Throw out negative fluxes, and convert fluxes to magnitudes." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "rCoaddCalib = butler.get('deepCoadd_calexp_calib', {'filter': 'HSC-R', 'tract': 0, 'patch': '1,1'})\n", + "iCoaddCalib = butler.get('deepCoadd_calexp_calib', {'filter': 'HSC-I', 'tract': 0, 'patch': '1,1'})\n", + "\n", + "rCoaddCalib.setThrowOnNegativeFlux(False)\n", + "iCoaddCalib.setThrowOnNegativeFlux(False)\n", + "\n", + "rMags = rCoaddCalib.getMagnitude(rSources['base_PsfFlux_flux'])\n", + "iMags = iCoaddCalib.getMagnitude(iSources['base_PsfFlux_flux'])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Select stars from catalog" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "deblended = rSources['deblend_nChild'] == 0\n", + "\n", + "refTable = butler.get('deepCoadd_ref', {'filter': 'HSC-R^HSC-I', 'tract': 0, 'patch': '1,1'})\n", + "inInnerRegions = refTable['detect_isPatchInner'] & refTable['detect_isTractInner'] # define inner regions\n", + "isSkyObject = refTable['merge_peak_sky'] # reject sky objects\n", + "isPrimary = refTable['detect_isPrimary']\n", + "\n", + "isStellar = iSources['base_ClassificationExtendedness_value'] < 1.\n", + "isGoodFlux = ~iSources['base_PsfFlux_flag']\n", + "selected = isPrimary & isStellar & isGoodFlux" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Make color-magnitude diagram." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.style.use('seaborn-notebook')\n", + "plt.figure(1, figsize=(4, 4), dpi=140)\n", + "plt.scatter(rMags[selected] - iMags[selected],\n", + " iMags[selected],\n", + " edgecolors='None', s=2, c='k')\n", + "plt.xlim(-0.5, 3)\n", + "plt.ylim(25, 14)\n", + "plt.xlabel('$r-i$')\n", + "plt.ylabel('$i$')\n", + "plt.subplots_adjust(left=0.125, bottom=0.1)\n", + "plt.show()" + ] } ], "metadata": { diff --git a/ImageProcessing/Re-RunHSC.sh b/ImageProcessing/Re-RunHSC.sh new file mode 100644 index 00000000..e5155312 --- /dev/null +++ b/ImageProcessing/Re-RunHSC.sh @@ -0,0 +1,131 @@ +: 'HSC Re-Run: Making Forced Photometry Light Curves from Scratch +Owner: **Justin Myles** (@jtmyles) +Last Verified to Run: **2018-09-05** +Verified Stack Release: **16.0** + +This project addresses issue #63: HSC Re-run + +This shell script runs the command-line tasks from the tutorial at pipelines.lsst.io for analysis +from raw images through source detection and forced photometry measurements. It is intended as an +intermediate step toward the end-goal of making a forced photometry lightcurve in the notebook at +StackClub/ImageProcessing/Re-RunHSC.ipynb + +Recommended to run with +$ bash Re-RunHSC.sh > output.txt +' + + +# Setup the LSST Stack +source /opt/lsst/software/stack/loadLSST.bash +eups list lsst_distrib +setup lsst_distrib + +# I. Setting up the Butler data repository + +date +echo "setup Butler" +setup -j -r /project/shared/data/ci_hsc +DATADIR="/home/$USER/DATA" +mkdir -p "$DATADIR" + +# A Butler needs a *mapper* file "to find and organize data in a format specific to each camera." +# We write this file to the data repository so that any instantiated Butler object knows which mapper to use. +echo lsst.obs.hsc.HscMapper > $DATADIR/_mapper + +# The injest script creates links in the instantiated butler repository to the original data files +ingestImages.py $DATADIR $CI_HSC_DIR/raw/*.fits --mode=link + +# Grab calibration files +installTransmissionCurves.py $DATADIR +ln -s $CI_HSC_DIR/CALIB/ $DATADIR/CALIB +mkdir -p $DATADIR/ref_cats +ln -s $CI_HSC_DIR/ps1_pv3_3pi_20170110 $DATADIR/ref_cats/ps1_pv3_3pi_20170110 + +date +echo "processCcd" +# II. Calibrate a single frame with processCcd.py +# Use calibration files to do CCD processing +processCcd.py $DATADIR --rerun processCcdOutputs --id + +# III. (omitted) Visualize images. + +# IV. Make coadds +date +echo "make coadds" +# IV. A. Make skymap +# A sky map is a tiling of the celestial sphere. It is composed of one or more tracts. +# A tract is composed of one or more overlapping patches. Each tract has a WCS. +# We define a skymap so that we can warp all of the exposure to fit on a single coordinate system +# This is a necessary step for making coadds + +date +echo "make skymap" +makeDiscreteSkyMap.py $DATADIR --id --rerun processCcdOutputs:coadd --config skyMap.projection="TAN" +# IV. B. Warp images onto skymap +date +echo "warp images" +makeCoaddTempExp.py $DATADIR --rerun coadd \ + --selectId filter=HSC-R \ + --id filter=HSC-R tract=0 patch=0,0^0,1^0,2^1,0^1,1^1,2^2,0^2,1^2,2 \ + --config doApplyUberCal=False doApplySkyCorr=False + + +makeCoaddTempExp.py $DATADIR --rerun coadd \ + --selectId filter=HSC-I \ + --id filter=HSC-I tract=0 patch=0,0^0,1^0,2^1,0^1,1^1,2^2,0^2,1^2,2 \ + --config doApplyUberCal=False doApplySkyCorr=False + + +# IV. C. Coadd warped images +# Now that we have warped images, we can perform coaddition to get deeper images +# The motivation for this is to have the deepest image possible for source detection +date +echo "coadd warped images" +assembleCoadd.py $DATADIR --rerun coadd \ + --selectId filter=HSC-R \ + --id filter=HSC-R tract=0 patch=0,0^0,1^0,2^1,0^1,1^1,2^2,0^2,1^2,2 + +assembleCoadd.py $DATADIR --rerun coadd \ + --selectId filter=HSC-I \ + --id filter=HSC-I tract=0 patch=0,0^0,1^0,2^1,0^1,1^1,2^2,0^2,1^2,2 + +# V. Measuring sources +date +echo "measure sources" +# V. A. Source detection +# As noted above, we do source detection on the deepest image possible. +echo "detect sources" +detectCoaddSources.py $DATADIR --rerun coadd:coaddPhot \ + --id filter=HSC-R tract=0 patch=0,0^0,1^0,2^1,0^1,1^1,2^2,0^2,1^2,2 + +detectCoaddSources.py $DATADIR --rerun coaddPhot \ + --id filter=HSC-I tract=0 patch=0,0^0,1^0,2^1,0^1,1^1,2^2,0^2,1^2,2 + +# V. B. Merge multi-band detection catalogs +# Ultimately, for photometry, we will need to deblend objects. +# In order to do this, we first merge the detected source catalogs. +date +echo "merge detection cats" +mergeCoaddDetections.py $DATADIR --rerun coaddPhot --id filter=HSC-R^HSC-I + +# V. C. Measure source catalogs on coadds +# Given a full source catalog, we can do regular photometry with implicit deblending. +date +echo "measure source cats on coadds" +measureCoaddSources.py $DATADIR --rerun coaddPhot --id filter=HSC-R +measureCoaddSources.py $DATADIR --rerun coaddPhot --id filter=HSC-I + +# V. D. Merge multi-band source catalogs from coadds +date +echo "merge source cats from coadds" +mergeCoaddMeasurements.py $DATADIR --rerun coaddPhot --id filter=HSC-R^HSC-I + +# V. E. Run forced photometry on coadds +# Given a full source catalog, we can do forced photometry with implicit deblending. +date +echo "run forcedphot on coadds" +forcedPhotCoadd.py $DATADIR --rerun coaddPhot:coaddForcedPhot --id filter=HSC-R +forcedPhotCoadd.py $DATADIR --rerun coaddForcedPhot --id filter=HSC-I + +# VI. Multi-band catalog analysis +# For analysis of the catalog, see part VI of StackClub/ImageProcessing/Re-RunHSC.ipynb diff --git a/Validation/image_quality_demo.ipynb b/Validation/image_quality_demo.ipynb deleted file mode 100644 index 13f4fcff..00000000 --- a/Validation/image_quality_demo.ipynb +++ /dev/null @@ -1,1135 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# LSST DM Stack Image Quality Walkthrough\n", - "\n", - "
Owner: **Keith Bechtol** ([@bechtol](https://github.com/LSSTScienceCollaborations/StackClub/issues/new?body=@bechtol))\n", - "
Last Verified to Run: **2018-08-03**\n", - "
Verified Stack Release: **v16.0**, **w201830**\n", - "\n", - "This notebook provides a walkthrough of several image quality measurements in the LSST DM Stack. The walkthrough is largely based on the following resources, which should be consulted for further details:\n", - "\n", - "* _The Hyper Suprime-Cam Software Pipeline_: https://arxiv.org/abs/1705.06766\n", - "* _The first-year shear catalog of the Subaru Hyper-Suprime-Cam SSP Survey_: https://arxiv.org/abs/1705.06745\n", - "* Systematics Tests In LEnsing (STILE): https://github.com/msimet/Stile\n", - "\n", - "This notebook uses a recent reprocessing of data from the first HSC Strategic Survey Program data release 1 (HSC-SSP DR1): \n", - "* _HSC Public Release Webpage_: https://hsc-release.mtk.nao.ac.jp/doc/\n", - "* _HSC-SSP DR1 Paper_: http://ads.nao.ac.jp/abs/2018PASJ...70S...8A\n", - "\n", - "### Learning Objectives\n", - "After working through and studying this notebook you should be able to\n", - " 1. Access PSF model sizes and ellipticities\n", - " 2. Access several different shape measurements with and without PSF correction for single-visit source catalogs and object catalogs derived from the full image stack\n", - " 3. Evaluate and visualize PSF model ellipticity residuals \n", - "\n", - "Other techniques that are demonstrated, but not empasized, in this notebook are\n", - " 1. Use `pandas` to compute aggregate statisitics and visualize the output \n", - " 2. Use `treecorr` to compute correlation functions\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." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Setup\n", - "You can find the Stack version by using `eups list -s` on the terminal command line." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# What version of the Stack am I using?\n", - "! echo $HOSTNAME\n", - "! eups list -s | grep lsst_distrib" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import numpy as np\n", - "from matplotlib import pyplot as plt\n", - "%matplotlib inline" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Matplotlib in Jupyter notebooks is a little fussy: let's ignore warnings from this and other modules." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import warnings\n", - "warnings.filterwarnings(\"ignore\", category=UserWarning) # eg warning from mpl.use(“Agg”)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Prelude: Data Sample\n", - "\n", - "Begin by identifying a data sample to use for illustration. In this example, we will use a recent reprocessing of the HSC-SSP DR1.\n", - "\n", - "Information available here: https://confluence.lsstcorp.org/display/DM/S18+HSC+PDR1+reprocessing\n", - "\n", - "The output repos are:\n", - "* /datasets/hsc/repo/rerun/DM-13666/UDEEP/\n", - "* /datasets/hsc/repo/rerun/DM-13666/DEEP/\n", - "* /datasets/hsc/repo/rerun/DM-13666/WIDE/\n", - "\n", - "Note that each of the data repositories contains all of the HSC visits, so one has to select by field to get the visits corresponding to a particular Strategic Survey Program (SSP) survey." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import lsst.daf.persistence as daf_persistence" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Instantiate the butler\n", - "band = 'HSC-I'\n", - "depth = 'WIDE' # WIDE, DEEP, UDEEP\n", - "field = 'SSP_WIDE' # SSP_WIDE, SSP_DEEP, SSP_UDEEP\n", - "butler = daf_persistence.Butler('/datasets/hsc/repo/rerun/DM-13666/%s/'%(depth))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Obtain a random dataid for a single-CCD image in the i-band\n", - "subset = butler.subset('calexp', dataId={'filter':band, 'field':field})\n", - "dataid = subset.cache[0]\n", - "print(dataid)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Use that dataid to get the image and associated source catalog\n", - "calexp = butler.get('calexp', **dataid)\n", - "src = butler.get('src', **dataid)\n", - "\n", - "# Note that the calexp contains a PSF object\n", - "psf = calexp.getPsf()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Part 1: PSF Model" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "shape = psf.computeShape()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### PSF Size\n", - "\n", - "This walkthrough focuses on the computation of PSF models and shape measurements using adaptive moments. Two common methods of computing the PSF size from adaptive moments are the trace radius and determinant radius.\n", - "\n", - "Adaptive second moments:\n", - "\n", - "$Q = \\begin{bmatrix}\n", - " I_{xx} & I_{xy} \\\\\n", - " I_{xy} & I_{yy}\n", - " \\end{bmatrix}$\n", - " \n", - "The trace radius is defined as $\\rm{Tr}(Q) = \\sqrt{\\frac{I_{xx} + I_{yy}}{2}}$ and determinant radius is defined as $\\rm{Det(Q)} = (I_{xx} I_{yy} - I_{xy}I_{xy})^\\frac{1}{4}$" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "i_xx, i_yy, i_xy = shape.getIxx(), shape.getIyy(), shape.getIxy()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Put in a few assert statements to show exactly what the Stack is doing" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Trace radius\n", - "assert np.isclose(shape.getTraceRadius(), np.sqrt((i_xx + i_yy) / 2.))\n", - "print('trace radius =', shape.getTraceRadius())" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Determinant radius\n", - "assert np.isclose(shape.getDeterminantRadius(), (i_xx * i_yy - i_xy**2)**(1. / 4.))\n", - "print('determinant radius =', shape.getDeterminantRadius())" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Comment Regarding Units\n", - "\n", - "What are the units of the trace radius and determinant radius? How do we convert the trace radius or determinant radius into FWHM in arcseconds? To do this, we need to quickly preview some functionality available from the calibrated image and source catalog.\n", - "\n", - "First, let's make sure that the PSF shape second moments from the PSF model in the calibrated exposure (calexp) match the entry for the source catalog for the PSF at the source position." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from lsst.afw.geom import Point2D\n", - "\n", - "# PSF shape at coordinates of first indexed souce\n", - "shape = psf.computeShape(Point2D(src['slot_Centroid_x'][0], src['slot_Centroid_y'][0])) \n", - "print(shape.getIxx())\n", - "# Compare to catalog entry for PSF shape\n", - "print(src['slot_PsfShape_xx'][0])" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Pixel scale in arcseconds\n", - "calexp.getWcs().getPixelScale().asArcseconds()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Units of the second moments\n", - "#src.schema.find('base_SdssShape_psf_xx').getField().getUnits()\n", - "src.schema.find('slot_PsfShape_xx').getField().getUnits()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We find that the trace radius and determinant radius are in units of pixels, and we found the pixel scale in arcseconds\n", - "\n", - "For a one-dimensional gaussian, the FWHM = $2 \\sqrt{2 \\ln 2}$ RMS = 2.35 RMS \n", - "\n", - "**THIS APPROXIMATION FOR THE FWHM SHOULD BE VERIFIED**" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Approximate FWHM in arcseconds\n", - "fwhm = 2 * np.sqrt(2. * np.log(2)) * shape.getTraceRadius() * calexp.getWcs().getPixelScale().asArcseconds()\n", - "print('FWHM = %.3f arcsec'%(fwhm))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### PSF Ellipticity\n", - "\n", - "**IMPORTANT:** Two conventions are commonly used in weak gravitational lensing to define ellipticity. The LSST Stack uses the \n", - "\"distortion\" convention as oppose to the \"shear\" convention. See https://en.wikipedia.org/wiki/Gravitational_lensing_formalism#Measures_of_ellipticity for both definitions.\n", - "\n", - "The Science Requirements Document (https://docushare.lsst.org/docushare/dsweb/Get/LPM-17) also uses the distortion convention.\n", - "\n", - "Formalism ellipticities in the distortion convention:\n", - "\n", - "$e_1 = \\frac{I_{xx} - I_{yy}}{I_{xx} + I_{yy}}$\n", - "\n", - "$e_2 = \\frac{2 I_{xy}}{I_{xx} + I_{yy}}$\n", - "\n", - "$\\tan(2 \\theta) = \\frac{2 I_{xy}}{I_{xx} - I_{yy}} $\n", - "\n", - "$e = \\sqrt{e_1^2 + e_2^2}$" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Put in a few assert statements to explicitly show what the Stack is doing" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Ellipticity\n", - "from lsst.afw.geom.ellipses import Quadrupole, SeparableDistortionTraceRadius\n", - "\n", - "q = Quadrupole(i_xx, i_yy, i_xy)\n", - "s = SeparableDistortionTraceRadius(q)\n", - "assert np.isclose(s.getE1(), (i_xx - i_yy) / (i_xx + i_yy)) # e1\n", - "print('e1 =', s.getE1())\n", - "assert np.isclose(s.getE2(), (2. * i_xy) / (i_xx + i_yy)) # e2\n", - "print('e2 =', s.getE2())\n", - "assert np.isclose(s.getEllipticity().getTheta(), np.arctan2(2. * i_xy, i_xx - i_yy) / 2.) # theta\n", - "print('theta =', s.getEllipticity().getTheta())\n", - "\n", - "# An alternative way to compute the angle\n", - "e1, e2 = s.getE1(), s.getE2() \n", - "assert np.allclose(np.arctan2(e2, e1) / 2., np.arctan2(2. * i_xy, i_xx - i_yy) / 2.)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "For visualization purposes, let's evaluate the PSF model at grid of points across the image" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "n = 100\n", - "x_array = np.arange(0, calexp.getDimensions()[0], 200)\n", - "y_array = np.arange(0, calexp.getDimensions()[1], 200)\n", - "xx, yy = np.meshgrid(x_array, y_array)\n", - "\n", - "print(calexp.getDimensions())\n", - "\n", - "size = []\n", - "i_xx = []\n", - "i_yy = []\n", - "i_xy = []\n", - "for x, y in zip(xx.flatten(), yy.flatten()):\n", - " point = Point2D(x, y)\n", - " shape = psf.computeShape(point)\n", - " size.append(shape.getTraceRadius())\n", - " i_xx.append(shape.getIxx())\n", - " i_yy.append(shape.getIyy())\n", - " i_xy.append(shape.getIxy())\n", - "size = np.reshape(size, xx.shape)\n", - "i_xx = np.reshape(i_xx, xx.shape)\n", - "i_yy = np.reshape(i_yy, xx.shape)\n", - "i_xy = np.reshape(i_xy, xx.shape)\n", - "\n", - "theta = np.arctan2(2. * i_xy, i_xx - i_yy) / 2.\n", - "e1 = (i_xx - i_yy) / (i_xx + i_yy)\n", - "e2 = (2. * i_xy) / (i_xx + i_yy)\n", - "theta_alternate = np.arctan2(e2, e1) / 2.\n", - "assert np.allclose(theta, theta_alternate)\n", - "e = np.sqrt(e1**2 + e2**2)\n", - "ex = e * np.cos(theta)\n", - "ey = e * np.sin(theta)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "and then plot the results in a few different ways" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plt.figure(figsize=(20, 6))\n", - "plt.subplots_adjust(wspace=0.5)\n", - "\n", - "plt.subplot(1, 4, 1)\n", - "plt.quiver(xx, yy, ex, ey, e, headlength=0., headwidth=1., pivot='mid', width=0.01)\n", - "#plt.quiver(xx, yy, scale=e, angles=phi, headlength=0., headwidth=1., pivot='mid', width=0.005)\n", - "#colorbar = plt.colorbar(label='r$\\sqrt(e1^{2} + e2^{2})$')\n", - "colorbar = plt.colorbar(label='e')\n", - "plt.clim(0., 0.05)\n", - "plt.xlabel('x')\n", - "plt.ylabel('y')\n", - "plt.title('Ellipticity Sticks')\n", - "\n", - "plt.subplot(1, 4, 2)\n", - "plt.scatter(xx, yy, c=e1, vmin=-0.05, vmax=0.05, cmap='bwr')\n", - "colorbar = plt.colorbar(label='e1')\n", - "plt.xlabel('x')\n", - "plt.ylabel('y')\n", - "plt.title('Value of e1')\n", - "\n", - "plt.subplot(1, 4, 3)\n", - "plt.scatter(xx, yy, c=e2, vmin=-0.05, vmax=0.05, cmap='bwr')\n", - "colorbar = plt.colorbar(label='e2')\n", - "plt.xlabel('x')\n", - "plt.ylabel('y')\n", - "plt.title('Value of e2')\n", - "\n", - "plt.subplot(1, 4, 4)\n", - "plt.scatter(xx, yy, c=size)\n", - "colorbar = plt.colorbar(label='Trace Radius')\n", - "plt.xlabel('x')\n", - "plt.ylabel('y')\n", - "plt.title('Ellipticity Modulus');" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Part 2: Source Catalog Shape Measurements\n", - "\n", - "The source catalogs include many measurements of source shapes and PSF model shapes. Let's look at the available columns relating to object shape:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# I'm still learning how to use regular expressions...\n", - "#import re\n", - "#list(filter(re.compile(\".*shape*\"), src.schema.getOrderedNames()))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Columns for shape measurements\n", - "for name in src.schema.getOrderedNames():\n", - " if 'shape' in name.lower():\n", - " print(name)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Columns for calibration measurements, including indicators of the stars using for PSF modeling and held in reserve for testing the PSF model." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Columns for calibration measurements\n", - "for name in src.schema.getOrderedNames():\n", - " if 'calib' in name.lower():\n", - " print(name)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "There are also \"slots\" defined for default measurements" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# How are the slots defined?\n", - "for slot in ['slot_PsfShape', 'slot_Shape']:\n", - " print('%s -> %s'%(slot, src.schema.getAliasMap()[slot]))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "To figure out what these columns mean, look at https://arxiv.org/abs/1705.06745\n", - "\n", - "**Non PSF-corrected shapes**\n", - "\n", - "SDSS algorithm:\n", - "\n", - "`base_SdssShape_xx/xy/yy` = Adaptive moments in arcsec$^2$, SDSS algorithm\n", - "\n", - "`base_SdssShape_psf_xx/xy/yy` = Adaptive moments of PSF evaluated at object position in arcsec$^2$, SDSS algorithm\n", - "\n", - "HSM algorithm: \n", - "\n", - "_Galaxy shapes are estimated on the coadded i-band images using the re-Gaussianization PSF correction method (Hirata & Seljak 2003)... In the course of the re-Gaussianization PSF-correction method, corrections are applied to account for dilution of the observed shape by the PSF, including the non-Gaussianity of both the PSF and the galaxy surface brightness profiles._\n", - "\n", - "`ext_shapeHSM_HsmPsfMoments_xx/xy/yy` = Adaptive moments of PSF evaluated at object position in arcsec$^2$, HSM algorithm\n", - "\n", - "`ext_shapeHSM_HsmSourceMoments_xx/xy/yy` = Adaptive moments in arcsec$^2$, not PSF-corrected , HSM algorithm\n", - "\n", - "`ext_shapeHSM_HsmSourceMomentsRound_xx/xy/yy` = ??\n", - "\n", - "**Regaussianization shapes based on data alone, PSF-corrected**\n", - "\n", - "`ext_shapeHSM_HsmShapeRegauss_e1/e2` = distortion in sky coordinates estimated by regaussianization method defined in distortion, HSM algorithm (**NEED TO CHECK IF THESE ARE SKY OR INSTRUMENT COORDINATES**)\n", - "\n", - "`ext_shapeHSM_HsmShapeRegauss_sigma` = non-calibrated shape measurement noise, HSM algorithm\n", - "\n", - "`ext_shapeHSM_HsmShapeRegauss_resolution` = resolution of galaxy image defined in equation (4) of https://arxiv.org/abs/1705.06745, HSM algorithm" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "With those definitions in mind, now verify that the PSF models from SDSS and HSM algorithms are consistent by comparing adaptive moments of the PSF. The plot below shows the ratio of adaptive moment sizes is very nearly unity. The plot also shows no trend in adaptive moments with respect to how well the source is resolved (i.e., source size), which is as expected since we are comparing PSF models as oppose to measurements for individual sources." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plt.figure()\n", - "plt.scatter(src['ext_shapeHSM_HsmShapeRegauss_resolution'], \n", - " src['base_SdssShape_psf_xx'] / src['ext_shapeHSM_HsmPsfMoments_xx'], s=1, label='xx')\n", - "plt.scatter(src['ext_shapeHSM_HsmShapeRegauss_resolution'], \n", - " src['base_SdssShape_psf_yy'] / src['ext_shapeHSM_HsmPsfMoments_yy'], s=1, label='yy')\n", - "plt.xlim(0., 1.)\n", - "plt.ylim(1. - 2.e-3, 1. + 2.e-3)\n", - "plt.legend(loc='upper left', markerscale=4)\n", - "plt.xlabel('ext_shapeHSM_HsmShapeRegauss_resolution')\n", - "plt.ylabel('base_SdssShape_psf / ext_shapeHSM_HsmPsfMoments');" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "To illustrate the difference between the PSF-corrected and non-PSF-corrected quantities within the HSM algorithm, now compare the ellipticity computed from the adaptive moments (not-PSF corrected) with the ellipticity including PSF-corrections. In the plot below, we see a clear trend with the amount of source resolution. For sources that are increasingly resolved, the PSF-corrected and non-PSF-corrected ellipticities converge. For sources that are not well resolved, the PSF-corrected and non-PSF-corrected ellipticities diverge." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "i_xx, i_yy, i_xy = src['ext_shapeHSM_HsmSourceMoments_xx'], src['ext_shapeHSM_HsmSourceMoments_yy'], src['ext_shapeHSM_HsmSourceMoments_xy']\n", - "e1_non_psf_corrected = (i_xx - i_yy) / (i_xx + i_yy)\n", - "ratio = e1_non_psf_corrected / src['ext_shapeHSM_HsmShapeRegauss_e1']\n", - "#trace_radius = np.sqrt((src['base_SdssShape_xx'] + src['base_SdssShape_yy']) / 2.)\n", - "\n", - "plt.figure()\n", - "plt.scatter(src['ext_shapeHSM_HsmShapeRegauss_resolution'], ratio, s=1)\n", - "plt.xlim(0., 1.)\n", - "plt.ylim(0., 1.2)\n", - "plt.xlabel('ext_shapeHSM_HsmShapeRegauss_resolution')\n", - "plt.ylabel('HsmSourceMoments_e1 / HsmShapeRegauss_e1');" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Intermission\n", - "\n", - "In the next section of the notebook, we'll use a full tract worth of data for some of the examples. It will take about 10 minutes to ingest the object catalogs from each patch in the tract. If you don't mind waiting, set the SHORTCUT variable below to `False`." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import pandas as pd" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Browse the deepCoadd_meas schema\n", - "dataid = {'filter':'HSC-I', 'tract':9936, 'patch':'0,0'}\n", - "coadd_meas = butler.get('deepCoadd_meas', dataId=dataid)\n", - "#coadd_meas.schema\n", - "#coadd_meas['slot_ModelFlux_fluxSigma']" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Set SHORTCUT = True for quick evaluation but lower statistics, \n", - "# SHORTCUT = False to get all the objects from all the patches in the tract (~10 mins)\n", - "# If you don't take the shortcut, youll need to do some uncommenting below as well.\n", - "\n", - "SHORTCUT = True\n", - "\n", - "skymap = butler.get('deepCoadd_skyMap')\n", - "\n", - "# Pick a random tract and collect all the patches\n", - "tract = 9936\n", - "patch_array = []\n", - "for ii in range(skymap.generateTract(tract).getNumPatches()[0]):\n", - " for jj in range(skymap.generateTract(tract).getNumPatches()[1]):\n", - " patch_array.append('%s,%s'%(ii, jj))\n", - "tract_array = np.tile(tract, len(patch_array))\n", - "\n", - "if SHORTCUT:\n", - " # Only get three patches\n", - " df_tract_patch = pd.DataFrame({'tract': [9936, 9936, 9936],\n", - " 'patch': ['0,0', '0,1', '0,2']})\n", - "else:\n", - " # Get all the object catalogs from one tract\n", - " df_tract_patch = pd.DataFrame({'tract': tract_array,\n", - " 'patch': patch_array})" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "selected_columns = ['id',\n", - " 'coord_ra', \n", - " 'coord_dec',\n", - " 'base_SdssCentroid_x',\n", - " 'base_SdssCentroid_y',\n", - " 'calib_psfCandidate',\n", - " 'calib_psfUsed',\n", - " 'base_ClassificationExtendedness_value',\n", - " 'base_ClassificationExtendedness_flag',\n", - " 'ext_shapeHSM_HsmPsfMoments_x',\n", - " 'ext_shapeHSM_HsmPsfMoments_y',\n", - " 'ext_shapeHSM_HsmPsfMoments_xx',\n", - " 'ext_shapeHSM_HsmPsfMoments_yy',\n", - " 'ext_shapeHSM_HsmPsfMoments_xy',\n", - " 'ext_shapeHSM_HsmPsfMoments_flag',\n", - " 'ext_shapeHSM_HsmPsfMoments_flag_no_pixels',\n", - " 'ext_shapeHSM_HsmPsfMoments_flag_not_contained',\n", - " 'ext_shapeHSM_HsmPsfMoments_flag_parent_source',\n", - " 'ext_shapeHSM_HsmShapeRegauss_e1',\n", - " 'ext_shapeHSM_HsmShapeRegauss_e2',\n", - " 'ext_shapeHSM_HsmShapeRegauss_sigma',\n", - " 'ext_shapeHSM_HsmShapeRegauss_resolution',\n", - " 'ext_shapeHSM_HsmShapeRegauss_flag',\n", - " 'ext_shapeHSM_HsmShapeRegauss_flag_no_pixels',\n", - " 'ext_shapeHSM_HsmShapeRegauss_flag_not_contained',\n", - " 'ext_shapeHSM_HsmShapeRegauss_flag_parent_source',\n", - " 'ext_shapeHSM_HsmShapeRegauss_flag_galsim',\n", - " 'ext_shapeHSM_HsmSourceMoments_x',\n", - " 'ext_shapeHSM_HsmSourceMoments_y',\n", - " 'ext_shapeHSM_HsmSourceMoments_xx',\n", - " 'ext_shapeHSM_HsmSourceMoments_yy',\n", - " 'ext_shapeHSM_HsmSourceMoments_xy',\n", - " 'ext_shapeHSM_HsmSourceMoments_flag',\n", - " 'ext_shapeHSM_HsmSourceMoments_flag_no_pixels',\n", - " 'ext_shapeHSM_HsmSourceMoments_flag_not_contained',\n", - " 'ext_shapeHSM_HsmSourceMoments_flag_parent_source',\n", - " 'slot_Centroid_x',\n", - " 'slot_Centroid_y',\n", - " 'slot_Shape_xx',\n", - " 'slot_Shape_yy',\n", - " 'slot_Shape_xy',\n", - " 'slot_PsfShape_xx',\n", - " 'slot_PsfShape_yy',\n", - " 'slot_PsfShape_xy',\n", - " 'slot_ModelFlux_flux',\n", - " 'slot_ModelFlux_fluxSigma',\n", - " 'slot_PsfFlux_flux',\n", - " 'slot_PsfFlux_fluxSigma']\n", - " " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "%%time\n", - "\n", - "coadd_array = []\n", - "for ii in range(0, len(df_tract_patch)):\n", - " tract, patch = df_tract_patch['tract'][ii], df_tract_patch['patch'][ii] \n", - " print(tract, patch)\n", - " dataid = {'filter':'HSC-I', 'tract':tract, 'patch':patch}\n", - " coadd_ref = butler.get('deepCoadd_ref', dataId=dataid)\n", - " coadd_meas = butler.get('deepCoadd_meas', dataId=dataid)\n", - " coadd_calib = butler.get('deepCoadd_calexp_calib', dataId=dataid)\n", - " \n", - " selected_rows = coadd_ref['detect_isPrimary'] \\\n", - " & ~coadd_meas['base_SdssCentroid_flag'] \\\n", - " & ~coadd_meas['base_PixelFlags_flag_interpolated'] \\\n", - " & ~coadd_meas['base_PixelFlags_flag_saturated'] \\\n", - " & ~coadd_meas['base_PsfFlux_flag'] \\\n", - " & ~coadd_meas['modelfit_CModel_flag']\n", - " \n", - " #print(len(selected_rows))\n", - " #print(coadd_meas.asAstropy().to_pandas().shape)\n", - " \n", - " coadd_array.append(coadd_meas.asAstropy().to_pandas().loc[selected_rows, selected_columns])\n", - " coadd_array[-1]['detect_isPrimary'] = coadd_ref['detect_isPrimary'][selected_rows]\n", - " \n", - " coadd_calib.setThrowOnNegativeFlux(False)\n", - " coadd_array[-1]['psf_mag'] = coadd_calib.getMagnitude(coadd_meas['base_PsfFlux_flux'][selected_rows])\n", - " coadd_array[-1]['cm_mag'] = coadd_calib.getMagnitude(coadd_meas['modelfit_CModel_flux'][selected_rows])\n", - " \n", - "df_coadd = pd.concat(coadd_array)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "If you just waited 10 minutes to get all the patches, it could be a good idea to write out your tract-level catalog for re-use later. \n", - "\n", - "> Note that you will need to do a `pip install --user tables` to write a pandas dataframe to HDF5 in the following cell." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "if not SHORTCUT:\n", - " import os\n", - " outfile = '%s/DATA/temp.h5'%(os.environ.get('HOME'))\n", - " print(outfile)\n", - " df_coadd.to_hdf(outfile, 'df')\n", - " \n", - "# Then, to read in this file, do:\n", - "\n", - "# df_coadd = pd.read_hdf(outfile, 'df')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "df_coadd.columns" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plt.figure()\n", - "plt.hexbin(np.degrees(df_coadd['coord_ra']), np.degrees(df_coadd['coord_dec']))\n", - "#plt.hexbin(df_coadd['base_SdssCentroid_x'], df_coadd['base_SdssCentroid_y'])\n", - "plt.xlim(plt.xlim()[::-1])\n", - "plt.xlabel('RA (deg)')\n", - "plt.ylabel('Dec (deg)');" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Part 3: Characterization\n", - "\n", - "Show a few examples of how to characterize shape measurements." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Ellipticity residuals for PSF calibration stars and reserve stars\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "def ellipticity(i_xx, i_yy, i_xy):\n", - " e1 = (i_xx - i_yy) / (i_xx + i_yy)\n", - " e2 = (2. * i_xy) / (i_xx + i_yy)\n", - " return e1, e2" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Ellipticity residuals\n", - "e1_star, e2_star = ellipticity(df_coadd['slot_Shape_xx'], df_coadd['slot_Shape_yy'], df_coadd['slot_Shape_xy'])\n", - "e1_psf, e2_psf = ellipticity(df_coadd['slot_PsfShape_xx'], df_coadd['slot_PsfShape_yy'], df_coadd['slot_PsfShape_xy'])\n", - "e1_res, e2_res = e1_star - e1_psf, e2_star - e2_psf " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "\"\"\"\n", - "# Plot to inform S/N selection\n", - "plt.figure()\n", - "plt.scatter(df_coadd['psf_mag'], df_coadd['psf_mag'] - df_coadd['cm_mag'], \n", - " c=df_coadd['slot_PsfFlux_flux'] / df_coadd['slot_PsfFlux_fluxSigma'], \n", - " vmin=1, vmax=100, s=1)\n", - "plt.colorbar()\n", - "plt.xlim(16., 30.)\n", - "plt.ylim(-0.1, 1.)\n", - "\"\"\";" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "cut_star = ((df_coadd['slot_PsfFlux_flux'] / df_coadd['slot_PsfFlux_fluxSigma']) > 100.) & (df_coadd['base_ClassificationExtendedness_value'] == 0)\n", - "\n", - "plt.figure()\n", - "#plt.yscale('log')\n", - "plt.hist(e1_res[cut_star] / e1_psf[cut_star], bins=np.linspace(-1., 1., 51), range=(-1., 1.), histtype='step', lw=2, density=True, label='e1')\n", - "plt.hist(e2_res[cut_star] / e2_psf[cut_star], bins=np.linspace(-1., 1., 51), range=(-1., 1.), histtype='step', lw=2, density=True, label='e2')\n", - "plt.xlim(-1., 1.)\n", - "plt.legend(loc='upper right')\n", - "plt.xlabel('e_star - e_psf / e_psf')\n", - "plt.ylabel('PDF');" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Repeat the exercise with the single-sensor, but compare the stars used for PSF modeling and those held in reserve for testing. In this case, there are only ~100 stars total on the single sensor so not enough statistics to see any more interesting trends. In practice, one would want to do this test with many more stars." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Ellipticity residuals\n", - "e1_star, e2_star = ellipticity(src['slot_Shape_xx'], src['slot_Shape_yy'], src['slot_Shape_xy'])\n", - "e1_psf, e2_psf = ellipticity(src['slot_PsfShape_xx'], src['slot_PsfShape_yy'], src['slot_PsfShape_xy'])\n", - "e1_res, e2_res = e1_star - e1_psf, e2_star - e2_psf " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "cut_psf_star = src['calib_psfUsed']\n", - "cut_reserved_star = src['calib_psf_reserved']\n", - "\n", - "print('Number of stars used for PSF modeling =', np.sum(cut_psf_star))\n", - "print('Number of reserved stars =', np.sum(cut_reserved_star)) # Should be about 20% of the sample used for PSF modeling\n", - "\n", - "plt.figure()\n", - "#plt.yscale('log')\n", - "plt.hist(e1_res[cut_psf_star] / e1_psf[cut_psf_star], bins=np.linspace(-1., 1., 21), range=(-1., 1.), histtype='step', lw=2, density=True, label='PSF')\n", - "plt.hist(e1_res[cut_reserved_star] / e1_psf[cut_reserved_star], bins=np.linspace(-1., 1., 21), range=(-1., 1.), histtype='step', lw=2, density=True, label='Reserved')\n", - "plt.xlim(-1., 1.)\n", - "plt.legend(loc='upper right')\n", - "plt.xlabel('e_star - e_psf / e_psf')\n", - "plt.ylabel('PDF');" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Compute aggregate statistics across the focal plane / tract\n", - "\n", - "This example shows how to use pandas to compute aggregate statistics across the focal plane. In this case, we plot the mean PSF model size (trace radius) computed in spatial bins over the sensor." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "#Convert the src object to a pandas dataframe\n", - "df_sensor = src.asAstropy().to_pandas()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Define a helper function to aggregate statistics in two-dimensions." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "def aggregateStatistics2D(df, x_name, x_bins, y_name, y_bins, z_name, operation):\n", - " \"\"\"\n", - " This function preserves empty bins, which get filled by nan values\n", - " \"\"\"\n", - " x_centers = x_bins[:-1] + 0.5 * np.diff(x_bins)\n", - " y_centers = y_bins[:-1] + 0.5 * np.diff(y_bins)\n", - " grouped = df.groupby([pd.cut(df[y_name], bins=y_bins),\n", - " pd.cut(df[x_name], bins=x_bins)])\n", - " xx, yy = np.meshgrid(x_centers, y_centers)\n", - " zz = grouped['id', z_name].agg(operation)[z_name].values.reshape(xx.shape)\n", - " return xx, yy, zz" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Define a new column for the PSF model trace radius\n", - "df_sensor['psf_trace_radius'] = np.sqrt((src['slot_PsfShape_xx'] + src['slot_PsfShape_yy']) / 2.)\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The next plot should look the same as the plot above that shows the PSF model size evaluated in a grid across the sensor" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "bin_size = 200 # pixels\n", - "x_bins = np.arange(0, calexp.getDimensions()[0] + bin_size, bin_size)\n", - "y_bins = np.arange(0, calexp.getDimensions()[1] + bin_size, bin_size)\n", - "\n", - "xx, yy, zz = aggregateStatistics2D(df_sensor, \n", - " 'slot_Centroid_x', x_bins, \n", - " 'slot_Centroid_y', y_bins, \n", - " 'psf_trace_radius', 'mean')\n", - "\n", - "plt.figure()\n", - "plt.pcolormesh(xx, yy, zz)\n", - "plt.colorbar(label='Trace Radius')\n", - "plt.xlabel('x')\n", - "plt.ylabel('y')\n", - "plt.xlim(0, calexp.getDimensions()[0])\n", - "plt.ylim(0, calexp.getDimensions()[1]);" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Now repeat the aggregaion exercise for the coadd tract." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Define a new column for the PSF model trace radius\n", - "df_coadd['psf_trace_radius'] = np.sqrt((df_coadd['ext_shapeHSM_HsmPsfMoments_xx'] + df_coadd['ext_shapeHSM_HsmPsfMoments_yy']) / 2.)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "bin_size = 200 # pixels\n", - "x_bins = np.arange(np.min(df_coadd.base_SdssCentroid_x), np.max(df_coadd.base_SdssCentroid_x) + bin_size, bin_size)\n", - "y_bins = np.arange(np.min(df_coadd.base_SdssCentroid_y), np.max(df_coadd.base_SdssCentroid_y) + bin_size, bin_size)\n", - "\n", - "xx, yy, zz = aggregateStatistics2D(df_coadd, \n", - " 'base_SdssCentroid_x', x_bins, \n", - " 'base_SdssCentroid_y', y_bins, \n", - " 'psf_trace_radius', 'mean')\n", - "\n", - "plt.figure()\n", - "plt.pcolormesh(xx, yy, zz)\n", - "plt.colorbar(label='Trace Radius')\n", - "plt.xlabel('x')\n", - "plt.ylabel('y')\n", - "plt.xlim(np.min(df_coadd.base_SdssCentroid_x), np.max(df_coadd.base_SdssCentroid_x))\n", - "plt.ylim(np.min(df_coadd.base_SdssCentroid_y), np.max(df_coadd.base_SdssCentroid_y));" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Ellipticity residual correlation functions\n", - "\n", - "Example of how to evaluate correlation functions, in this case, for ellipticity residuals on the coadd. This section reproduces the analysis of `validate_drp` (https://github.com/lsst/validate_drp) in a condensed form for easier reading." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import treecorr\n", - "import astropy.units as u" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Ellipticity residuals\n", - "e1_star, e2_star = ellipticity(df_coadd['slot_Shape_xx'], df_coadd['slot_Shape_yy'], df_coadd['slot_Shape_xy'])\n", - "e1_psf, e2_psf = ellipticity(df_coadd['slot_PsfShape_xx'], df_coadd['slot_PsfShape_yy'], df_coadd['slot_PsfShape_xy'])\n", - "e1_res, e2_res = e1_star - e1_psf, e2_star - e2_psf " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Select high S/N stars\n", - "cut_star = ((df_coadd['slot_PsfFlux_flux'] / df_coadd['slot_PsfFlux_fluxSigma']) > 100.) & (df_coadd['base_ClassificationExtendedness_value'] == 0)\n", - "cut_star = cut_star & ~np.isnan(e1_res) & ~np.isnan(e2_res)\n", - "\n", - "nbins=20\n", - "min_sep=0.25\n", - "max_sep=20\n", - "sep_units='arcmin'\n", - "verbose=False\n", - "\n", - "catTree = treecorr.Catalog(ra=df_coadd['coord_ra'][cut_star], dec=df_coadd['coord_dec'][cut_star], \n", - " g1=e1_res[cut_star], g2=e2_res[cut_star],\n", - " dec_units='radian', ra_units='radian')\n", - "gg = treecorr.GGCorrelation(nbins=nbins, min_sep=min_sep, max_sep=max_sep,\n", - " sep_units=sep_units,\n", - " verbose=verbose)\n", - "gg.process(catTree)\n", - "r = np.exp(gg.meanlogr) * u.arcmin\n", - "xip = gg.xip * u.Unit('')\n", - "xip_err = np.sqrt(gg.varxi) * u.Unit('')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plt.figure()\n", - "plt.xscale('log')\n", - "plt.yscale('log')\n", - "plt.errorbar(r.value, xip, yerr=xip_err)\n", - "plt.xlabel('Separation (arcmin)')\n", - "plt.ylabel('Median Residual Ellipticity Correlation');" - ] - } - ], - "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 -}