Measurement/DwarfGalaxySrcOverdensity.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Exploring source measurements in a dwarf galaxy\n",
    "<br>Owner(s): **Jeff Carlin** ([@jeffcarlin](https://github.com/LSSTScienceCollaborations/StackClub/issues/new?body=@jeffcarlin))\n",
    "<br>Last Verified to Run: **2020-12-18**\n",
    "<br>Verified Stack Release: **21.0**\n",
    "\n",
    "### Learning Objectives:\n",
    "\n",
    "After working through this tutorial you should be able to: \n",
    "1. Figure out how to find and extract data around a given RA, Dec position;\n",
    "2. Extract a cutout image of an object of interest, and a catalog of sources in the cutout footprint;\n",
    "3. Apply quality cuts to the source catalog, and plot CMDs and spatial maps of sources;\n",
    "4. Diagnose possible problems with the photometry in crowded regions;\n",
    "5. Apply your understanding of the root cause of this problem to re-run photometry with different pipeline processes set up;\n",
    "6. Extract well-measured photometry in the crowded region of interest, and proceed with using these data for science.\n",
    "\n",
    "### Logistics\n",
    "This notebook is intended to be runnable on `https://lsst-lsp-stable.ncsa.illinois.edu/` from a local git clone of https://github.com/LSSTScienceCollaborations/StackClub.\n",
    "\n",
    "## Set-up"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next few cells give you some options for your \"Set-up\" section - you may not need them all."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll need the `stackclub` package to be installed. If you are not developing this package, you can install it using `pip`, like this:\n",
    "```\n",
    "pip install git+git://github.com/LSSTScienceCollaborations/StackClub.git#egg=stackclub\n",
    "```\n",
    "If you are developing the `stackclub` package (e. g., by adding modules to it to support the Stack Club tutorial that you are writing), you'll need to make a local, editable installation. In the top level folder of the `StackClub` repo, do:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "! cd .. && python setup.py -q develop --user && cd -"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When editing the `stackclub` package files, we want the latest version to be imported when we re-run the import command. To enable this, we need the %autoreload magic command."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "You can find the Stack version that this notebook is running by using eups list -s on the terminal command line:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# What version of the Stack am I using?\n",
    "! echo $HOSTNAME\n",
    "! eups list -s | grep lsst_distrib"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "For this tutorial we'll need the following modules:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import glob\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from astropy.table import vstack\n",
    "from astropy.table import Table\n",
    "from astropy.coordinates import SkyCoord\n",
    "from astropy import units as u\n",
    "from astropy.visualization import ZScaleInterval\n",
    "from astropy.io import ascii"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import lsst.daf.persistence as dafPersist\n",
    "import lsst.geom as geom\n",
    "import lsst.afw.image as afwImage\n",
    "import lsst.afw.table as afwTable\n",
    "import lsst.afw.display as afwDisplay\n",
    "import lsst.pipe.base as pipeBase\n",
    "from lsst.meas.base import NoiseReplacer, NoiseReplacerConfig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up some plotting defaults:\n",
    "\n",
    "params = {\n",
    "   'axes.labelsize': 28,\n",
    "   'font.size': 24,\n",
    "   'legend.fontsize': 14,\n",
    "   'xtick.major.width': 3,\n",
    "   'xtick.minor.width': 2,\n",
    "   'xtick.major.size': 12,\n",
    "   'xtick.minor.size': 6,\n",
    "   'xtick.direction': 'in',\n",
    "   'xtick.top': True,\n",
    "   'lines.linewidth':3,\n",
    "   'axes.linewidth':3,\n",
    "   'axes.labelweight':3,\n",
    "   'axes.titleweight':3,\n",
    "   'ytick.major.width':3,\n",
    "   'ytick.minor.width':2,\n",
    "   'ytick.major.size': 12,\n",
    "   'ytick.minor.size': 6,\n",
    "   'ytick.direction': 'in',\n",
    "   'ytick.right': True,\n",
    "   'figure.figsize': [9, 8]\n",
    "   }\n",
    "\n",
    "plt.rcParams.update(params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are going to look at a known dwarf galaxy that is present in the HSC-SSP data release, but for which (as you'll see), the photometry is poor. This galaxy is called \"KKR 25\" -- here is more information about KKR 25 on NED: \n",
    "\n",
    "https://ned.ipac.caltech.edu/byname?objname=KKR%2025&hconst=67.8&omegam=0.308&omegav=0.692&wmap=4&corr_z=1\n",
    "\n",
    "To help with finding the correct tract, the tract IDs of HSC data can be found at: https://hsc-release.mtk.nao.ac.jp/doc/index.php/database/\n",
    "\n",
    "The following data discovery section follows this demo notebook: https://github.com/lsst-sqre/notebook-demo/blob/u/jcarlin/DM-15123/hsc_cosmos_comparison.ipynb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Find and read the data\n",
    "\n",
    "We will point to the shared repository containing processed HSC-SSP data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "depth = 'DEEP' #'WIDE' -- which survey depth (DEEP, WIDE, ULTRA-DEEP) do you want?\n",
    "butler = dafPersist.Butler('/datasets/hsc/repo/rerun/DM-13666/%s/'%(depth))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the Gen2 Butler, the easiest way to know what tracts and patches are present is to parse the filenames. This will be easier with the Gen3 Butler."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "infiles = glob.glob('/datasets/hsc/repo/rerun/DM-13666/%s/deepCoadd-results/HSC-I/*/*/calexp-*.fits'%(depth))\n",
    "# How many tract/patch files are in this repo?\n",
    "print(len(infiles))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to find which tracts/patches our target is in, we first loop over all the files in the repository and extract their central RA, Dec coordinates into a pandas dataframe.\n",
    "\n",
    "This involves reading the \"skyMap\" objects from the butler repository, and using the WCS and bounding box information to extract the coordinates for each patch."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def getPatchCenter(skymap, tract, patch):\n",
    "    # Find the center of the patch BBox in pixel coordinates and convert with tract WCS\n",
    "    patch = list(map(int, patch.split(',')))\n",
    "    wcs = skymap.generateTract(tract).getWcs()\n",
    "    bbox = skymap.generateTract(tract).getPatchInfo(patch).getInnerBBox()\n",
    "    sky = wcs.pixelToSky(np.mean((bbox.getBeginX(), bbox.getEndX())), \n",
    "                         np.mean((bbox.getBeginY(), bbox.getEndY())))\n",
    "    return sky[0].asDegrees(), sky[1].asDegrees()\n",
    "\n",
    "tract_array = []\n",
    "patch_array = []\n",
    "\n",
    "# Each line of \"infiles\" looks like this (infile[0]):\n",
    "# '/datasets/hsc/repo/rerun/DM-13666/DEEP/deepCoadd-results/HSC-I/17130/8,8/calexp-HSC-I-17130-8,8.fits'\n",
    "\n",
    "# Split on the '/', then parse the resulting array. Here's the result of infiles[0].split('/'):\n",
    "#['',\n",
    "# 'datasets',\n",
    "# 'hsc',\n",
    "# 'repo',\n",
    "# 'rerun',\n",
    "# 'DM-13666',\n",
    "# 'DEEP',\n",
    "# 'deepCoadd-results',\n",
    "# 'HSC-I',\n",
    "# '17130',\n",
    "# '8,8',\n",
    "# 'calexp-HSC-I-17130-8,8.fits']\n",
    "\n",
    "# So the tract number is third from the end, and patch is second from the end.\n",
    "\n",
    "for ii in range(0, len(infiles)):\n",
    "    parts = infiles[ii].split('/')\n",
    "    tract_array.append(int(parts[-3]))\n",
    "    patch_array.append(parts[-2])\n",
    "\n",
    "print('Found %i patches'%(len(patch_array)))\n",
    "    \n",
    "skymap = butler.get('deepCoadd_skyMap')\n",
    "    \n",
    "ra_center_array = np.tile(0., len(patch_array))\n",
    "dec_center_array = np.tile(0., len(patch_array))\n",
    "for ii, (tract, patch) in enumerate(zip(tract_array, patch_array)):\n",
    "    ra_center_array[ii], dec_center_array[ii] = getPatchCenter(skymap, tract, patch)\n",
    "    \n",
    "# Assemble in pandas data frame    \n",
    "data = {'tract': tract_array,\n",
    "        'patch': patch_array, \n",
    "        'ra_center': ra_center_array, \n",
    "        'dec_center': dec_center_array}\n",
    "df_tract_patch = pd.DataFrame(data=data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Where are these data on the sky? Let's plot a map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot RA, Dec of patch centers:\n",
    "\n",
    "fig=plt.figure()\n",
    "plt.scatter(df_tract_patch['ra_center'], df_tract_patch['dec_center'], marker='.', s=5)\n",
    "# From looking it up on NED, KKR 25 is in the \"ELAIS-N1\" field.\n",
    "# Uncomment the xlim, ylim below to zoom in on that region.\n",
    "#plt.xlim(238., 246.)\n",
    "#plt.ylim(52., 58.)\n",
    "plt.gca().invert_xaxis() # flip the x-axis so the plot is \"sky-right\"\n",
    "plt.xlabel('RA (deg)')\n",
    "plt.ylabel('Dec (deg)')\n",
    "plt.minorticks_on()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Locate KKR 25 within the patches with HSC data\n",
    "\n",
    "To figure out which patch(es) KKR 25 is located within, we will use Astropy \"SkyCoord\" objects. It is then straightforward to calculate the separation of KKR 25 from the centers of all the patches we extracted above. Then we'll find the closest patches, which (hopefully) contain KKR 25."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create \"SkyCoord\" objects for the central coordinate of KKR25 (or CetusIII, VirgoI?), and the patch centers of HSC:\n",
    "kkr25_pos = SkyCoord(243.448333*u.degree, 54.37111*u.degree, frame='icrs') # KKR25\n",
    "#cosmos_patch_cen = SkyCoord(180.04*u.degree, 0.68*u.degree, frame='icrs') # VirgoI\n",
    "hsc_patch_cen = SkyCoord(df_tract_patch['ra_center'].values*u.degree, df_tract_patch['dec_center'].values*u.degree, frame='icrs')\n",
    "\n",
    "# Calculate the separation (in degrees) of each HSC patch from KKR25:\n",
    "sep_patches = hsc_patch_cen.separation(kkr25_pos)\n",
    "\n",
    "# Pick patches within 'size_select' deg of KKR25:\n",
    "size_select = 0.25 # in degrees\n",
    "selection_patch = (sep_patches.degree < size_select) \n",
    "print('Number of patches selected: ',np.sum(selection_patch))\n",
    "\n",
    "# Add the separation (in arcminutes) to the dataframe:\n",
    "df_tract_patch['sep_from_kkr25'] = sep_patches.arcmin\n",
    "\n",
    "# Extract just the selected (nearby) patches:\n",
    "df_selected = df_tract_patch[selection_patch]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Print some info about the selected patches:\n",
    "df_selected"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From this table, it looks like KKR 25 is nearest the center of (tract, patch)=(17130, '1,4'). \n",
    "\n",
    "For context about the size of patches, HSC pixels are ~0.168\" per pixel, and a patch is 4k x 4k. Thus each patch is about (4k * 0.168\") = 672 arcsec = 11.2 arcmin = 0.1867 deg. on a side."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Confirm that KKR 25 is visible in the images of these patches\n",
    "\n",
    "To do this, we will create \"data references\" -- for coadd images, these are based on the filter, tract, and patch. These three dimensions specify a unique data ID, and are fed to the butler as a dictionary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataref14 = {'tract':17130, 'patch':'1,4', 'filter':'HSC-G'}\n",
    "\n",
    "# Note that we can extract the information directly from the dataframe we made above rather\n",
    "#   than explicity typing the values, like so:\n",
    "#dataref14 = {'tract':df_selected.iloc[4]['tract'],'patch':df_selected.iloc[4]['patch'],\n",
    "#             'filter':'HSC-G'}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you want to know what types of data *may be* present in the butler repository, you can uncomment the following line for a list. BEWARE: it's a long list!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#butler.getDatasetTypes()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, that only lists dataset types that are defined. It doesn't actually tell you if they exist *in this repository*. To find that out, you can use the \"datasetExists\" method. (For example, here we'll check whether the \"forced_src\" catalog exists for dataref14.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "butler.datasetExists('deepCoadd_forced_src', dataId=dataref14)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note, however, that in the Gen2 butler, this doesn't _guarantee_ that the catalog actually exists as a file on disk. It only tells us that some processing steps were executed that intended to create the catalog, and thus added its entry to the registry. This issue will be remedied in the Gen3 butler, so it is not important to dwell on it further."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get the images and catalogs for the two patches of interest:\n",
    "\n",
    "In particular, we'll use the forced photometry (\"deepCoadd_forced_src\") and the deep coadd images (\"deepCoadd_calexp\"). \n",
    "\n",
    "Notice that the *calexp* has a WCS and a \"PhotoCalib\" associated with it - we'll extract these for later use (if we want to convert XY pixel values to RA, Dec, or fluxes to magnitudes)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Full patch image\n",
    "image14 = butler.get('deepCoadd_calexp', dataId=dataref14)\n",
    "# Because an entire tract shares a WCS, the corner of the patch (or cutout) isn't necessarily at (X,Y)=(0,0). Get the XY0 pixel values:\n",
    "xy0_14 = image14.getXY0() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zscale = ZScaleInterval()\n",
    "\n",
    "plt.rcParams.update({'figure.figsize' : (9, 8)})\n",
    "\n",
    "vmin14, vmax14 = zscale.get_limits(image14.image.array)\n",
    "# Get the dimensions of the image so we can set plot limits\n",
    "imsize14 = image14.getDimensions()\n",
    "plt.imshow(image14.image.array, vmin=vmin14, vmax=vmax14, cmap='binary')\n",
    "\n",
    "# Set the plot range to the dimensions:\n",
    "plt.xlim(0,imsize14[0])\n",
    "plt.ylim(0,imsize14[1])\n",
    "plt.colorbar()\n",
    "plt.xlabel('X (pix)')\n",
    "plt.ylabel('Y (pix)')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "KKR 25 is the semi-resolved blob of stars at the top, partially obscured by a bright star."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create an image (from multiple patches) centered on KKR 25:\n",
    "(borrowing heavily from `ip_diffim`)\n",
    "\n",
    "Maybe instead of being limited by the patch boundaries, we would like an image that is centered on the dwarf. This will require putting together images from multiple patches"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# A useful function we'll use in a little bit.\n",
    "#   Calculate the number of pixels in a given N-arcmin area.\n",
    "def arcmin2pix(narcmin):\n",
    "    arcsecPerPix = 0.168  # HSC pixels are 0.168\"/pixel\n",
    "    arcminPerPix = arcsecPerPix / 60.0\n",
    "    npix = narcmin/arcminPerPix\n",
    "    return int(np.rint(npix))    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This function creates a cutout image centered on *raDec*, of size *size_arcmin*, with the center of the image in the patch denoted by *dataref*. It does the following:\n",
    "1. Extract the skymap, calexp, and WCS.\n",
    "2. Create an empty bounding box of the specified size, centered on the given position.\n",
    "3. Taking the corners of the bounding box, use the function defined earlier to look up which patches overlap that bounding box.\n",
    "4. Extract the subregions (using the \"deepCoadd_calexp_sub\" object) from the overlapping patches, and insert them into the bounding box we created.\n",
    "5. Return the cutout image as an ExposureF."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cutoutImage(raDec, size_arcmin, dataref):\n",
    "    '''\n",
    "    Function to extract an image cutout.\n",
    "    Inputs:\n",
    "      raDec: an afwGeom.SpherePoint object with desired RA, Dec center \n",
    "      size_arcmin: desired image size (in arcminutes)\n",
    "      filter: desired filter of extracted image (string)\n",
    "    '''\n",
    "    # Start with a center position and make a bbox centered there, and expand it to include data from neighboring patches:\n",
    "    skyMap =  butler.get(\"deepCoadd_skyMap\", immediate=True, dataId=dataref)\n",
    "    tmp_calexp = butler.get(\"deepCoadd_calexp\", immediate=True, dataId=dataref)\n",
    "    expWcs = tmp_calexp.getWcs()\n",
    "\n",
    "    size_pixels = arcmin2pix(size_arcmin)\n",
    "    print('npixels: ',size_pixels)\n",
    "    cutoutSize = geom.ExtentI(size_pixels, size_pixels) # size of cutout in pixels\n",
    "    print('cutoutSize: ',cutoutSize)\n",
    "    print(expWcs.skyToPixel(raDec))\n",
    "    xy = geom.Point2I(expWcs.skyToPixel(raDec)) # central XY coordinate of our star's RA, Dec position\n",
    "\n",
    "    # Create the bounding box:\n",
    "    bbox = geom.Box2I(xy - cutoutSize//2, cutoutSize)\n",
    "    bbox = geom.Box2D(bbox)\n",
    "\n",
    "    tractInfo = skyMap.findTract(raDec)\n",
    "    skyCorners = [expWcs.pixelToSky(pixPos) for pixPos in bbox.getCorners()]\n",
    "    patchList = tractInfo.findPatchList(skyCorners)\n",
    "\n",
    "    # compute coadd bbox\n",
    "    coaddWcs = tractInfo.getWcs()\n",
    "    coaddBBox = geom.Box2D()\n",
    "    for skyPos in skyCorners:\n",
    "        coaddBBox.include(coaddWcs.skyToPixel(skyPos))\n",
    "    coaddBBox = geom.Box2I(coaddBBox)\n",
    "\n",
    "    # assemble coadd exposure from subregions of patches\n",
    "    coaddExposure = afwImage.ExposureF(coaddBBox, coaddWcs)\n",
    "    coaddExposure.maskedImage.set(np.nan, afwImage.Mask.getPlaneBitMask(\"NO_DATA\"), np.nan)\n",
    "    nPatchesFound = 0\n",
    "    coaddFilter = None\n",
    "    for patchInfo in patchList:\n",
    "        patchSubBBox = patchInfo.getOuterBBox()\n",
    "        patchSubBBox.clip(coaddBBox)\n",
    "        patchArgDict = dict(\n",
    "                            datasetType=\"deepCoadd_calexp_sub\",\n",
    "                            bbox=patchSubBBox,\n",
    "                            tract=tractInfo.getId(),\n",
    "                            patch=\"%s,%s\" % (patchInfo.getIndex()[0], patchInfo.getIndex()[1]),\n",
    "                            filter = dataref['filter'],\n",
    "                            )\n",
    "        print(\"Reading patch %s\" % patchArgDict)\n",
    "        coaddPatch = butler.get(**patchArgDict)\n",
    "        nPatchesFound += 1\n",
    "        coaddExposure.maskedImage.assign(coaddPatch.maskedImage, coaddPatch.getBBox())\n",
    "        if coaddFilter is None:\n",
    "            coaddFilter = coaddPatch.getFilter()\n",
    "\n",
    "    if nPatchesFound == 0:\n",
    "        raise RuntimeError(\"No patches found!\")\n",
    "\n",
    "    coaddExposure.setFilter(coaddFilter)\n",
    "    return pipeBase.Struct(exposure=coaddExposure, sources=None)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now define a similar function, but creating a catalog of the sources in the cutout region."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cutoutCatalog(raDec, size_arcmin, dataref):\n",
    "    '''\n",
    "    Function to extract an catalog cutout.\n",
    "    Inputs:\n",
    "      raDec: an afwGeom.SpherePoint object with desired RA, Dec center \n",
    "      size_arcmin: desired image size (in arcminutes)\n",
    "      filter: desired filter of extracted image (string)\n",
    "    '''\n",
    "    # Start with a center position and make a bbox centered there, and expand it to include data from neighboring patches:\n",
    "    skyMap =  butler.get(\"deepCoadd_skyMap\", immediate=True, dataId=dataref)\n",
    "    tmp_calexp = butler.get(\"deepCoadd_calexp\", immediate=True, dataId=dataref)\n",
    "    expWcs = tmp_calexp.getWcs()\n",
    "\n",
    "    size_pixels = arcmin2pix(size_arcmin)\n",
    "    cutoutSize = geom.ExtentI(size_pixels, size_pixels) # size of cutout in pixels\n",
    "    xy = geom.Point2I(expWcs.skyToPixel(raDec)) # central XY coordinate of our star's RA, Dec position\n",
    "    \n",
    "    # Create the bounding box:\n",
    "    bbox = geom.Box2I(xy - cutoutSize//2, cutoutSize)\n",
    "    bbox = geom.Box2D(bbox)\n",
    "\n",
    "    tractInfo = skyMap.findTract(raDec)\n",
    "    skyCorners = [expWcs.pixelToSky(pixPos) for pixPos in bbox.getCorners()]\n",
    "    patchList = tractInfo.findPatchList(skyCorners)\n",
    "    print('Npatches: ',len(patchList))\n",
    "\n",
    "    nPatchesFound = 0\n",
    "    for patchInfo in patchList:\n",
    "        patchArgDict = dict(\n",
    "                            datasetType=\"deepCoadd_forced_src\",\n",
    "                            tract=tractInfo.getId(),\n",
    "                            patch=\"%s,%s\" % (patchInfo.getIndex()[0], patchInfo.getIndex()[1]),\n",
    "                            filter = dataref['filter'],\n",
    "                            flags=afwTable.SOURCE_IO_NO_FOOTPRINTS\n",
    "                            )\n",
    "        print(\"Reading patch %s\" % patchArgDict)\n",
    "        src = butler.get(**patchArgDict)\n",
    "        nPatchesFound += 1\n",
    "\n",
    "        patchArgDict_ref = dict(\n",
    "                                datasetType=\"deepCoadd_ref\",\n",
    "                                tract=tractInfo.getId(),\n",
    "                                patch=\"%s,%s\" % (patchInfo.getIndex()[0], patchInfo.getIndex()[1]),\n",
    "                                filter = dataref['filter'],\n",
    "                                flags=afwTable.SOURCE_IO_NO_FOOTPRINTS\n",
    "                                )\n",
    "        print(\"Reading ref cat for patch %s\" % patchArgDict_ref)\n",
    "        refcat = butler.get(**patchArgDict_ref)\n",
    "\n",
    "        patchArgDict_calib = dict(\n",
    "                                datasetType=\"deepCoadd_calexp_photoCalib\",\n",
    "                                tract=tractInfo.getId(),\n",
    "                                patch=\"%s,%s\" % (patchInfo.getIndex()[0], patchInfo.getIndex()[1]),\n",
    "                                filter = dataref['filter'],\n",
    "                                flags=afwTable.SOURCE_IO_NO_FOOTPRINTS\n",
    "                                )\n",
    "        print(\"Reading calib for patch %s\" % patchArgDict_calib)\n",
    "        calib = butler.get(**patchArgDict_calib)\n",
    "        psfmags = calib.instFluxToMagnitude(src, 'base_PsfFlux')\n",
    "        modelmags = calib.instFluxToMagnitude(src, 'modelfit_CModel')\n",
    "        \n",
    "        # Create an empty Astropy Table, then insert the magnitudes into it.\n",
    "        mags = Table()\n",
    "        mags['psfmag'] = psfmags\n",
    "        mags['modelmag'] = modelmags\n",
    "\n",
    "\n",
    "        if nPatchesFound == 1:\n",
    "            print('patch: 1')\n",
    "            tmp_src = src\n",
    "            tmp_ref = refcat\n",
    "            tmp_mags = mags\n",
    "            tmp_calib = calib\n",
    "        else:\n",
    "            print('patch: ',nPatchesFound)\n",
    "            tmp_src.extend(src)\n",
    "            tmp_ref.extend(refcat)\n",
    "            tmp_mags = vstack([tmp_mags,mags])\n",
    "            \n",
    "        print('length of tmp_src, tmp_ref, tmp_mags: ',len(tmp_src),len(tmp_ref), len(tmp_mags))\n",
    "\n",
    "\n",
    "    print('Making deep copies')\n",
    "    cat_final = tmp_src.copy(deep = True)\n",
    "    ref_final = tmp_ref.copy(deep = True)\n",
    "\n",
    "\n",
    "    if nPatchesFound == 0:\n",
    "        raise RuntimeError(\"No patches found!\")\n",
    "\n",
    "    return cat_final, ref_final, tmp_mags"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "size_arcmin = 3.0\n",
    "\n",
    "# Create a SpherePoint object at the center of KKR 25:\n",
    "raDec = geom.SpherePoint(kkr25_pos.ra.value, kkr25_pos.dec.value, geom.degrees)\n",
    "\n",
    "# Extract the cutout image:\n",
    "imgCutout = cutoutImage(raDec, size_arcmin, dataref14)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice from the outputs of the previous cell that two patches were needed to create the image. Now, how does it look? (Can we clearly see the dwarf galaxy in it?)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams.update({'figure.figsize' : (9, 8)})\n",
    "\n",
    "# Use lsst.afw.display with the matplotlib backend to display the image\n",
    "afwDisplay.setDefaultBackend('matplotlib') \n",
    "plt.figure()\n",
    "afw_display = afwDisplay.Display()\n",
    "afw_display.scale('linear', 'zscale')\n",
    "afw_display.mtv(imgCutout.exposure.image)\n",
    "plt.xticks([])\n",
    "plt.yticks([])\n",
    "# plt.savefig('kkr25_image.png') # Uncomment if you want to save the image.\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There's a beautiful, partially-resolved galaxy there. Unfortunately, there's a bright star right next to it. We'll have to see how that affects the measurements from the Stack.\n",
    "\n",
    "What does the mask plane look like in this region?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "afwDisplay.setDefaultBackend('matplotlib') \n",
    "plt.figure()\n",
    "afw_display = afwDisplay.Display()\n",
    "afw_display.scale('linear', 'zscale')\n",
    "afw_display.mtv(imgCutout.exposure.maskedImage)\n",
    "plt.xticks([])\n",
    "plt.yticks([])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The MaskedImage is overplotting the values of the mask plane. \n",
    "# Extract the mapping between mask color and mask bit directly from the mask object:\n",
    "\n",
    "for maskName, maskBit in imgCutout.exposure.mask.getMaskPlaneDict().items():\n",
    "    print('{}: {}'.format(afw_display.getMaskPlaneColor(maskName),maskName))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OK, so this is telling us that the stars in KKR 25 (because they are blue in the mask) are flagged as DETECTED, but are mostly beneath the BAD mask from the bright star. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extract concatenated source catalogs in the same footprint using the function we defined above. (Note: this returns catalog data for the entire patches that are part of the requested footprint, and not just within the radius specified.) We'll grab three bands (gri).\n",
    "\n",
    "The *cutoutCatalog* function returns the forced photometry (\"deepCoadd_forced_src\"), the reference catalog of combined information for all bandpasses (\"deepCoadd_ref\"), and the calibrated magnitudes (PSF and cmodel magnitudes).\n",
    "\n",
    "*The following cell may take a few minutes to execute.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dr14_i = dataref14.copy()\n",
    "dr14_i['filter'] = 'HSC-I'\n",
    "dr14_r = dataref14.copy()\n",
    "dr14_r['filter'] = 'HSC-R'\n",
    "dr14_g = dataref14.copy()\n",
    "dr14_g['filter'] = 'HSC-G'\n",
    "cat_g, refcat_g, mags_g = cutoutCatalog(raDec, size_arcmin, dr14_g)\n",
    "cat_r, refcat_r, mags_r = cutoutCatalog(raDec, size_arcmin, dr14_r)\n",
    "cat_i, refcat_i, mags_i = cutoutCatalog(raDec, size_arcmin, dr14_i)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Make some selections based on the data flags:\n",
    "\n",
    " - detect_isPrimary (see https://pipelines.lsst.io/getting-started/multiband-analysis.html?#filtering-for-unique-deblended-sources-with-the-detect-isprimary-flag for a description of what that flag means), \n",
    " - require PSF fluxes to be measured (i.e., not 'NaN'), \n",
    " - point sources (based on ClassificationExtendedness)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "is_primary = (refcat_g['detect_isPrimary'] == True) & (refcat_r['detect_isPrimary'] == True) &\\\n",
    "             (refcat_i['detect_isPrimary'] == True) \n",
    "is_finite = (np.isfinite(cat_g.getPsfInstFlux())) & (np.isfinite(cat_r.getPsfInstFlux())) & (np.isfinite(cat_i.getPsfInstFlux()))\n",
    "ptsource = (cat_g['base_ClassificationExtendedness_flag'] == False) & (cat_r['base_ClassificationExtendedness_flag'] == False) & (cat_i['base_ClassificationExtendedness_flag'] == False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Scatter plots of the source positions with various subsets selected:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams.update({'figure.figsize' : (14, 7)})\n",
    "\n",
    "plt.subplot(121)\n",
    "plt.plot(np.rad2deg(cat_g['coord_ra'][is_primary & ptsource]),np.rad2deg(cat_g['coord_dec'][is_primary & ptsource]),'k.',alpha=0.05)\n",
    "plt.xlim(243.65,243.29)\n",
    "plt.xlabel('RA (deg)')\n",
    "plt.ylabel('Dec (deg)')\n",
    "plt.minorticks_on()\n",
    "plt.title('Deblended point sources')\n",
    "\n",
    "plt.subplot(122)\n",
    "plt.plot(np.rad2deg(cat_g['coord_ra'][is_primary]),np.rad2deg(cat_g['coord_dec'][is_primary]),'k.',alpha=0.05)\n",
    "plt.xlim(243.65,243.29)\n",
    "plt.title('All sources')\n",
    "plt.xlabel('RA (deg)')\n",
    "plt.minorticks_on()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Good - we can see a concentration of sources in the center, which correspond to KKR 25."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot a color-magnitude diagram (CMD) of objects near KKR 25:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because the stack uses PanSTARRS-1 (PS1) for photometric calibration, we will overplot a PS1 isochrone for comparison. This was downloaded from the Padova isochrones web interface (http://stev.oapd.inaf.it/cgi-bin/cmd_3.3)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read in the isochrone file:\n",
    "iso = ascii.read('/project/jcarlin/isochrones/parsec12_panstarrs1/0_000500.dat',header_start = 13)\n",
    "\n",
    "dist = 1.7e6 # in pc # Roughly the distance to KKR 25 (see the NED link above)\n",
    "ebv = 0.008 # (E(B-V) = 0.008 according to Makarov+2012) -- You could also look this up using the central coordinates.\n",
    "dmod = 5.0*np.log10(dist) - 5.0 # convert distance to distance modulus\n",
    "\n",
    "# Raw PS1-band isochrone magnitudes\n",
    "g_iso_ps1 = iso['gP1']\n",
    "r_iso_ps1 = iso['rP1']\n",
    "i_iso_ps1 = iso['iP1']\n",
    "\n",
    "# Extinction coefficients from Schlafly+2011:\n",
    "A_g = 3.172 * ebv\n",
    "A_r = 2.271 * ebv\n",
    "A_i = 1.682 * ebv\n",
    "\n",
    "# NOTE: we are *adding* extinction to the isochrone, instead of looking up the values for all stars. We'll work with\n",
    "#   un-dereddened magnitudes throughout.\n",
    "g_iso = g_iso_ps1 + dmod + A_g\n",
    "r_iso = r_iso_ps1 + dmod + A_r\n",
    "i_iso = i_iso_ps1 + dmod + A_i "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Select a smaller region around KKR 25:\n",
    "coord_g = SkyCoord(cat_g['coord_ra']*u.radian, cat_g['coord_dec']*u.radian, frame='icrs')\n",
    "### Notice in the line above that coordinates in the catalog are in radians by default.\n",
    "\n",
    "# Calculate the separation (in degrees) of each HSC patch from KKR25:\n",
    "sep_kkr25 = coord_g.separation(kkr25_pos)\n",
    "\n",
    "# Select within some small radius:\n",
    "sel_rad = 2.5 # radius in arcminutes\n",
    "sel_cen = (sep_kkr25.arcmin < sel_rad)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following uses the difference in PSF and cmodel magnitudes to select things that are \"star-like\". See the Bosch+2018 paper on the HSC pipeline (https://ui.adsabs.harvard.edu/abs/2018PASJ...70S...5B/abstract) for info about this. The \"classification_extendedness\" flag should have taken care of this already, but this is a way of refining that if you want a stricter selection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# difference in mags (PSF - cmodel)\n",
    "magdiff_i = mags_i['psfmag'][:,0]-mags_i['modelmag'][:,0]\n",
    "starsel_i = (np.abs(magdiff_i) < 0.05)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams.update({'figure.figsize' : (14, 7)})\n",
    "\n",
    "plt.figure()\n",
    "\n",
    "# The following selects only up to the RGB tip (i.e., no post-RGB, AGB, etc.) from the isochrone, for clarity\n",
    "okstage = (iso['stage'] <= 3)\n",
    "\n",
    "plt.subplot(121)\n",
    "plt.plot(mags_g['psfmag'][is_primary & sel_cen][:,0]-mags_i['psfmag'][is_primary & sel_cen][:,0],\\\n",
    "         mags_g['psfmag'][is_primary & sel_cen][:,0], 'ko', ms=5, alpha=0.1, label='all sources')\n",
    "plt.plot(mags_g['psfmag'][is_primary & starsel_i & sel_cen][:,0]-mags_i['psfmag'][is_primary & starsel_i & sel_cen][:,0],\\\n",
    "         mags_g['psfmag'][is_primary & starsel_i & sel_cen][:,0], 'ro', ms=5, alpha=0.75, label='point sources')\n",
    "plt.plot(g_iso[okstage]-i_iso[okstage], g_iso[okstage], 'k', label='__no label__')#s', ms=1)\n",
    "plt.xlabel('g-i')\n",
    "plt.ylabel('g')\n",
    "plt.legend()\n",
    "plt.xlim(-1,3)\n",
    "plt.ylim(33,21)\n",
    "plt.title('within 2.5 arcmin')\n",
    "plt.minorticks_on()\n",
    "\n",
    "plt.subplot(122)\n",
    "plt.scatter(mags_g['psfmag'][is_primary & sel_cen][:,0]-mags_i['psfmag'][is_primary & sel_cen][:,0],\\\n",
    "            mags_g['psfmag'][is_primary & sel_cen][:,0], c=sep_kkr25[is_primary & sel_cen].arcmin, cmap='Spectral_r', s=5)\n",
    "plt.plot(g_iso[okstage]-i_iso[okstage], g_iso[okstage], 'k')#s', ms=1)\n",
    "plt.xlabel('g-i')\n",
    "plt.ylabel('g')\n",
    "plt.xlim(-1,3)\n",
    "plt.ylim(33,21)\n",
    "plt.colorbar(label='dist from KKR 25 (arcmin)')\n",
    "plt.minorticks_on()\n",
    "plt.show()\n",
    "#plt.savefig('cmds_kkr25_noisereplace_on.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, we can see that most of the stars nearest the center of KKR 25 (the blue/green points above) have \"bogus\" photometry (i.e., g > 28th mag)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just to confirm that it is the stars near KKR 25 that are getting set to unrealistically faint magnitudes, let's plot the positions of objects classified as point sources, color-coded by their g-band PSF magnitude."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams.update({'figure.figsize' : (9, 8)})\n",
    "\n",
    "gmags0 = mags_g['psfmag']\n",
    "gmags0[np.isnan(gmags0)] = 0.0\n",
    "plt.scatter(np.rad2deg(cat_g['coord_ra'][is_primary & ptsource]),\\\n",
    "            np.rad2deg(cat_g['coord_dec'][is_primary & ptsource]),\\\n",
    "            c=gmags0[is_primary & ptsource][:,0],s=2,vmax=30,vmin=18)\n",
    "# NOTE: in the above line, gmags0 must be subscripted to select only the magnitude column. gmags0 contains both the\n",
    "#       magnitudes and their errors, so to get mags only, use gmags0[:0] (for errors, use gmags0[:1])\n",
    "\n",
    "plt.xlim(np.max(np.rad2deg(cat_g['coord_ra'][is_primary & ptsource])),\\\n",
    "         np.min(np.rad2deg(cat_g['coord_ra'][is_primary & ptsource])))\n",
    "plt.ylim(np.min(np.rad2deg(cat_g['coord_dec'][is_primary & ptsource])),\\\n",
    "         np.max(np.rad2deg(cat_g['coord_dec'][is_primary & ptsource])))\n",
    "plt.colorbar(label='g mag')\n",
    "plt.xlabel('RA (deg)')\n",
    "plt.ylabel('Dec (deg)')\n",
    "plt.minorticks_on()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Most of the sources in KKR25, while detected, seem to have magnitudes >= 30(!). This is unrealistically faint - for example, in the CMDs above, you can see that the faintest reliably-measured objects are at g~27.\n",
    "\n",
    "Plot magnitude vs. position to see if it is true that the g>28 sources are concentrated in KKR25."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(np.rad2deg(cat_g['coord_ra'][is_primary & ptsource]),\\\n",
    "         gmags0[is_primary & ptsource][:,0],'b.',alpha=0.1, label='point source')\n",
    "plt.plot(np.rad2deg(cat_g['coord_ra'][is_primary & ptsource & sel_cen]),\\\n",
    "         gmags0[is_primary & ptsource & sel_cen][:,0],'r.',alpha=0.25, label='near KKR 25')\n",
    "plt.ylim(18,38)\n",
    "plt.xlim(np.max(np.rad2deg(cat_g['coord_ra'][is_primary & ptsource])),\\\n",
    "         np.min(np.rad2deg(cat_g['coord_ra'][is_primary & ptsource])))\n",
    "plt.legend()\n",
    "plt.xlabel('RA (deg)')\n",
    "plt.ylabel('g')\n",
    "plt.minorticks_on()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sure enough, there is a big blob of objects at RA~243.45, with magnitudes fainter than g>28. So, it looks like the Stack *detects* sources in the center, but when it measures them, something goes wrong.\n",
    "\n",
    "After some exploration, we have figured out what happens in cases like this. After detecting sources in an image, the measurements are done by *first replacing the flux in all neighboring footprints with noise*, then measuring the flux of the source of interest. In crowded regions such as KKR 25, the footprints of neighboring sources almost always overlap the source of interest, meaning the source itself gets replaced with noise before running the measurement task. This means that the measurement returns only background noise, or a flux consistent with zero."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Try using the noise replacer to replace the detection footprints with noise\n",
    "\n",
    "Just to confirm that this is the issue, let's try using the \"Noise Replacer\" to replace footprints with noise in the image, to see how the footprints in some small region look.\n",
    "\n",
    "For more details on footprints and the noise replacer, see the Stack Club \"Footprints\" notebook in the Source Detection section: https://stackclub.readthedocs.io/en/latest/notebooks.html#source-detection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Source catalog:\n",
    "cat14 = butler.get('deepCoadd_forced_src', dataId=dr14_g)\n",
    "\n",
    "# Reference object:\n",
    "ref14 = butler.get('deepCoadd_ref', dataId=dr14_g)\n",
    "\n",
    "# Calibrated exposure:\n",
    "calexp14 = butler.get('deepCoadd_calexp', dataId=dr14_g)\n",
    "# We are going to manipulate calexp14, so make a clone of the original:\n",
    "calexp14_orig = calexp14.clone()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# extract the footprints of all detected sources into a dict:\n",
    "fp_dict_g = {measRecord.getId(): (measRecord.getParent(), measRecord.getFootprint()) for measRecord in ref14}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Instantiate the NoiseReplacer and feed it the calexp image and the footprints:\n",
    "nr_config = NoiseReplacerConfig()\n",
    "noiseReplacer = NoiseReplacer(nr_config, calexp14, fp_dict_g)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display the original image:\n",
    "plt.figure()\n",
    "afw_display = afwDisplay.Display()\n",
    "afw_display.setMaskTransparency(100)\n",
    "afw_display.scale('linear', 'zscale')\n",
    "afw_display.mtv(calexp14_orig.getImage()[5700:6700,19100:20100])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display the image with the source footprints replaced by noise:\n",
    "plt.figure()\n",
    "afw_display = afwDisplay.Display()\n",
    "afw_display.setMaskTransparency(100)\n",
    "afw_display.scale('linear', 'zscale')\n",
    "afw_display.mtv(calexp14.getImage()[5700:6700,19100:20100])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The fact that they subtracted off nicely confirms that the sources in KKR 25 were *detected*. Now the question is whether one can get *measurements* for those detected objects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract a small cutout image, offset slightly from the center of KKR 25, for efficiency:\n",
    "cutout_pos = SkyCoord(kkr25_pos.ra+(0.5*u.arcmin), kkr25_pos.dec+(0.5*u.arcmin), frame='icrs') # KKR25\n",
    "raDec_cutout = geom.SpherePoint(cutout_pos.ra.value, cutout_pos.dec.value, geom.degrees)\n",
    "# Extract the cutout image:\n",
    "cutout_size = 0.4 # in arcmin\n",
    "gcutout = cutoutImage(raDec_cutout, cutout_size, dr14_g)\n",
    "icutout = cutoutImage(raDec_cutout, cutout_size, dr14_i)\n",
    "gcat, grefcat, gcutout_mags = cutoutCatalog(raDec_cutout, cutout_size, dr14_g)\n",
    "icat, irefcat, icutout_mags = cutoutCatalog(raDec_cutout, cutout_size, dr14_i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# data quality cuts as we did before\n",
    "is_primary_cutout = (grefcat['detect_isPrimary'] == True) &\\\n",
    "                    (irefcat['detect_isPrimary'] == True) \n",
    "is_finite_cutout = (np.isfinite(gcat.getPsfInstFlux())) & (np.isfinite(icat.getPsfInstFlux()))\n",
    "magcut_cutout = (gcutout_mags['psfmag'][:,0] < 28.0)\n",
    "ptsource_cutout = (gcat['base_ClassificationExtendedness_flag'] == False) &\\\n",
    "                  (icat['base_ClassificationExtendedness_flag'] == False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Cut down the list of footprints from the entire calexp to only those in the small subregion of interest:\n",
    "# X: 5970-6130, Y: 19775-19905\n",
    "xmin = 5970\n",
    "xmax = 6130\n",
    "ymin = 19775\n",
    "ymax = 19905\n",
    "keep14 = (ref14.getX() < xmax) & (ref14.getX() > xmin) & (ref14.getY() < ymax) & (ref14.getY() > ymin)\n",
    "ref14keep = ref14[keep14]\n",
    "fp_dict_g14 = {measRecord.getId(): (measRecord.getParent(), measRecord.getFootprint()) for measRecord in ref14keep}\n",
    "srcids14 = ref14keep['id']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have a list (srcids14) of the ids for all the source footprints in the region of interest."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make another clone of the calexp image, in which we will:\n",
    "# (a) Replace pixels in all source footprints with noise, then\n",
    "# (b) Add a single source footprint's back into the image.\n",
    "\n",
    "c14 = calexp14_orig.clone()\n",
    "\n",
    "# Instantiate the NoiseReplacer and feed it the calexp image and the footprints:\n",
    "nr_config = NoiseReplacerConfig()\n",
    "noiseReplacer = NoiseReplacer(nr_config, c14, fp_dict_g14)\n",
    "\n",
    "# The following line simply selects an index number for which source (from srcids14) to visualize.\n",
    "# (Change the integer to visualize a different source.)\n",
    "ind = 17 # 53, 55, 56 are huge!, 21 is crazy huge!\n",
    "noiseReplacer.insertSource(srcids14[ind])\n",
    "print('(X,Y): ', ref14keep[ind].getX(), ref14keep[ind].getY())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract the image array for the selected footprint:\n",
    "fp = fp_dict_g14[srcids14[ind]][1].getSpans().unflatten(fp_dict_g14[srcids14[ind]][1].getImageArray())\n",
    "# Extract the bounding box for the footprint:\n",
    "fp_bbox = fp_dict_g14[srcids14[ind]][1].getBBox()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.colors import LogNorm\n",
    "plt.imshow(fp,\n",
    "           cmap='viridis', origin='lower', norm=LogNorm(vmin=1e-6, vmax=fp.max()))\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams.update({'figure.figsize' : (14, 7)})\n",
    "\n",
    "f, ax = plt.subplots(nrows=1, ncols=2, sharey=True)\n",
    "\n",
    "axflat = ax.flatten()\n",
    "ax1 = axflat[0]\n",
    "ax2 = axflat[1]\n",
    "data = calexp14_orig.getMaskedImage()[xmin:xmax,ymin:ymax].getImage().array\n",
    "ax1.minorticks_on()\n",
    "ax1.imshow(data, origin='lower',vmin=-0.01, vmax=0.5)\n",
    "ax1.plot(ref14keep[ind].getX()-xmin, ref14keep[ind].getY()-ymin, 'o', ms=30, color='white', fillstyle='none', markeredgewidth=4)\n",
    "ax1.set_title('original image')\n",
    "\n",
    "data = c14.getMaskedImage()[xmin:xmax,ymin:ymax].getImage().array\n",
    "ax2.imshow(data, origin='lower',vmin=-0.01, vmax=0.5)\n",
    "ax2.plot(ref14keep[ind].getX()-xmin, ref14keep[ind].getY()-ymin, 'o', ms=30, color='white', fillstyle='none', markeredgewidth=4)\n",
    "ax2.minorticks_on()\n",
    "ax2.set_title('footprints replaced w/ noise')\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This one looks like a reasonable point source. However, if you try different footprints, you will see that some of them span *large* areas, such that if you subtract them off before doing measurements, you will have nothing left to measure. Let's look at one of those large footprints:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make another clone of the calexp image, in which we will:\n",
    "# (a) Replace pixels in all source footprints with noise, then\n",
    "# (b) Add a single source footprint's back into the image.\n",
    "\n",
    "c14 = calexp14_orig.clone()\n",
    "\n",
    "# Instantiate the NoiseReplacer and feed it the calexp image and the footprints:\n",
    "noiseReplacer = NoiseReplacer(nr_config, c14, fp_dict_g14)\n",
    "\n",
    "# The following line simply selects an index number for which source (from srcids14) to visualize.\n",
    "# (Change the integer to visualize a different source.)\n",
    "ind_large = 56\n",
    "noiseReplacer.insertSource(srcids14[ind_large])\n",
    "print('(X,Y): ', ref14keep[ind_large].getX(), ref14keep[ind_large].getY())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract the image array for the selected footprint:\n",
    "fp_large = fp_dict_g14[srcids14[ind_large]][1].getSpans().unflatten(fp_dict_g14[srcids14[ind_large]][1].getImageArray())\n",
    "# Extract the bounding box for the footprint:\n",
    "fp_bbox_large = fp_dict_g14[srcids14[ind_large]][1].getBBox()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(fp_large,\n",
    "           cmap='viridis', origin='lower', norm=LogNorm(vmin=1e-6, vmax=fp_large.max()))\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Whoa, this thing is huge. It has connected the pixels of many detected sources into one contiguous footprint. (Notice that this footprint's area is larger than the 160x130-pixel cutout image that we selected.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams.update({'figure.figsize' : (14, 7)})\n",
    "\n",
    "f, ax = plt.subplots(nrows=1, ncols=2, sharey=True)\n",
    "\n",
    "axflat = ax.flatten()\n",
    "ax1 = axflat[0]\n",
    "ax2 = axflat[1]\n",
    "data = calexp14_orig.getMaskedImage()[xmin:xmax,ymin:ymax].getImage().array\n",
    "ax1.minorticks_on()\n",
    "ax1.imshow(data, origin='lower',vmin=-0.01, vmax=0.5)\n",
    "ax1.plot(ref14keep[ind_large].getX()-xmin, ref14keep[ind_large].getY()-ymin, 'o', ms=30, color='white', fillstyle='none', markeredgewidth=4)\n",
    "ax1.set_title('original image')\n",
    "\n",
    "data = c14.getMaskedImage()[xmin:xmax,ymin:ymax].getImage().array\n",
    "ax2.imshow(data, origin='lower',vmin=-0.01, vmax=0.1)\n",
    "ax2.plot(ref14keep[ind_large].getX()-xmin, ref14keep[ind_large].getY()-ymin, 'o', ms=30, color='white', fillstyle='none', markeredgewidth=4)\n",
    "ax2.minorticks_on()\n",
    "ax2.set_title('footprints replaced w/ noise')\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see in the right panel that the selected footprint fills almost the entire sub-image. Also note that there are resolved sources visible *within* that footprint, so that masking this footprint would make those sources unmeasureable. Clearly, these oversized footprints are the reason we get ~zero fluxes for stars in the crowded parts of KKR 25.\n",
    "\n",
    "One last thing to check: we can extract the area of each footprint, so that we could imagine testing whether this issue had happened in a given region by examining footprint sizes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check the area of the footprints we have selected:\n",
    "srctmp = ref14keep[ind]\n",
    "fprtmp = srctmp.getFootprint()\n",
    "srctmp_large = ref14keep[ind_large]\n",
    "fprtmp_large = srctmp_large.getFootprint()\n",
    "print(fprtmp.getArea(),' pixels^2, ',fprtmp_large.getArea(),' pixels^2')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Re-run forced photometry, but without replacing all the footprints with noise:\n",
    "\n",
    "For this section, we will re-run the photometry, turning off the feature that replaces neighboring footprints with noise. Rather than doing this in the notebook, I instead re-ran the relevant tasks from the command line. The following gives a brief summary of how this was done:\n",
    "\n",
    "Some recent processing of HSC-PDR2 resides at: /datasets/hsc/repo/rerun/DM-23243\n",
    "\n",
    "```\n",
    "source /opt/lsst/software/stack/loadLSST.bash\n",
    "setup -t v20_0_0 lsst_distrib\n",
    "setup obs_subaru\n",
    "\n",
    "multiBandDriver.py /datasets/hsc/repo --rerun DM-23243/DEEP:/project/jcarlin/kkr25/rerun/noise_replace_offPDR2 --id tract=17130 filter=HSC-G^HSC-I^HSC-R --batch-type=slurm --mpiexec='-bind-to socket' --job mbKKR25 --time 30000 --nodes 2 --procs 23 —configfile 'mbconfig.py' \n",
    "```\n",
    "\n",
    "The file 'mbconfig.py' contains the following:\n",
    "```\n",
    "config.measureCoaddSources.measurement.doReplaceWithNoise=False\n",
    "config.forcedPhotCoadd.measurement.doReplaceWithNoise=False\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once that has finished running, we can examine the outputs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "butler = dafPersist.Butler('/project/jcarlin/kkr25/rerun/noise_replace_offPDR2')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extract catalogs from this new repository in a small region centered on KKR 25:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "size_arcmin = 3.0\n",
    "\n",
    "# Create a SpherePoint object at the center of KKR 25:\n",
    "raDec = geom.SpherePoint(kkr25_pos.ra.value, kkr25_pos.dec.value, geom.degrees)\n",
    "\n",
    "# Tract/patch where KKR 25's center appears\n",
    "tract = 17130\n",
    "patch = '1,4'\n",
    "dr14_g = {'tract':tract, 'patch':patch, 'filter':'HSC-G'}\n",
    "dr14_r = {'tract':tract, 'patch':patch, 'filter':'HSC-R'}\n",
    "dr14_i = {'tract':tract, 'patch':patch, 'filter':'HSC-I'}\n",
    "\n",
    "cat_g, refcat_g, mags_g = cutoutCatalog(raDec, size_arcmin, dr14_g)\n",
    "cat_r, refcat_r, mags_r = cutoutCatalog(raDec, size_arcmin, dr14_r)\n",
    "cat_i, refcat_i, mags_i = cutoutCatalog(raDec, size_arcmin, dr14_i)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Apply quality cuts to select \"good,\" point-source-like objects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "is_primary = (refcat_g['detect_isPrimary'] == True) & (refcat_r['detect_isPrimary'] == True) &\\\n",
    "             (refcat_i['detect_isPrimary'] == True) \n",
    "is_finite = (np.isfinite(cat_g.getPsfInstFlux())) & (np.isfinite(cat_r.getPsfInstFlux())) & (np.isfinite(cat_i.getPsfInstFlux()))\n",
    "deblended = (cat_g['deblend_nChild'] == 0) & (cat_r['deblend_nChild'] == 0) & (cat_i['deblend_nChild'] == 0)\n",
    "ptsource = (cat_g['base_ClassificationExtendedness_flag'] == False) & (cat_r['base_ClassificationExtendedness_flag'] == False) & (cat_i['base_ClassificationExtendedness_flag'] == False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell is not totally necessary, but is a convenient way to change the radius of the selection cone."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Select a smaller region around KKR 25:\n",
    "coord_g = SkyCoord(cat_g['coord_ra']*u.radian, cat_g['coord_dec']*u.radian, frame='icrs')\n",
    "\n",
    "# Calculate the separation (in degrees) of each HSC patch from KKR25:\n",
    "sep_kkr25 = coord_g.separation(kkr25_pos)\n",
    "\n",
    "# Select within some small radius:\n",
    "sel_rad = 2.5\n",
    "sel_cen = (sep_kkr25.arcmin < sel_rad)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As before, we will select a more reliable set of \"stars\" using cuts on the difference between model and PSF magnitudes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# difference in mags (PSF - cmodel)\n",
    "magdiff_i = mags_i['psfmag'][:,0]-mags_i['modelmag'][:,0]\n",
    "starsel_i = (np.abs(magdiff_i) < 0.05)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Drumroll please...\n",
    "\n",
    "Now we are ready to check the CMD of sources measured with noise replacement turned off. Hopefully we will see a well-defined red giant branch of KKR 25..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams.update({'figure.figsize' : (14, 7)})\n",
    "\n",
    "plt.figure()\n",
    "\n",
    "okstage = (iso['stage'] <= 3)\n",
    "\n",
    "plt.subplot(121)\n",
    "plt.plot(mags_g['psfmag'][is_primary & deblended & sel_cen][:,0]-mags_i['psfmag'][is_primary & deblended & sel_cen][:,0],\\\n",
    "         mags_g['psfmag'][is_primary & deblended & sel_cen][:,0], 'ko', ms=5, alpha=0.1, label='all sources')\n",
    "plt.plot(mags_g['psfmag'][is_primary & deblended & starsel_i & sel_cen][:,0]-mags_i['psfmag'][is_primary & deblended & starsel_i & sel_cen][:,0],\\\n",
    "         mags_g['psfmag'][is_primary & deblended & starsel_i & sel_cen][:,0], 'ro', ms=5, alpha=0.75, label='point sources')\n",
    "plt.plot(g_iso[okstage]-i_iso[okstage], g_iso[okstage], 'k', label='__no label__')#s', ms=1)\n",
    "plt.xlabel('g-i')\n",
    "plt.ylabel('g')\n",
    "plt.legend()\n",
    "plt.xlim(-1,3)\n",
    "plt.ylim(33,21)\n",
    "plt.title('within 2.5 arcmin')\n",
    "plt.minorticks_on()\n",
    "\n",
    "plt.subplot(122)\n",
    "plt.scatter(mags_g['psfmag'][is_primary & sel_cen][:,0]-mags_i['psfmag'][is_primary & sel_cen][:,0],\\\n",
    "            mags_g['psfmag'][is_primary & sel_cen][:,0], c=sep_kkr25[is_primary & sel_cen].arcmin, cmap='Spectral_r', s=5)\n",
    "plt.plot(g_iso[okstage]-i_iso[okstage], g_iso[okstage], 'k')#s', ms=1)\n",
    "plt.xlabel('g-i')\n",
    "plt.ylabel('g')\n",
    "plt.xlim(-1,3)\n",
    "plt.ylim(33,21)\n",
    "plt.colorbar(label='dist from KKR 25 (arcmin)')\n",
    "plt.minorticks_on()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Success! The red points in the left panel clearly define a red giant branch consistent with the known distance to KKR 25 (as evidenced by the fact that the isochrone matches the data pretty well).\n",
    "\n",
    "Notice also that, unlike before, there are very few points fainter than g=28. This is a good thing, since those were all junk. Let's zoom in on the CMD and bask in the warm feelings of excitement to do science with our well-measured photometry of dwarf galaxy KKR 25:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams.update({'figure.figsize' : (7, 7)})\n",
    "\n",
    "plt.figure()\n",
    "\n",
    "okstage = (iso['stage'] <= 3)\n",
    "\n",
    "plt.plot(mags_g['psfmag'][is_primary & deblended & sel_cen][:,0]-mags_i['psfmag'][is_primary & deblended & sel_cen][:,0],\\\n",
    "         mags_g['psfmag'][is_primary & deblended & sel_cen][:,0], 'ko', ms=5, alpha=0.1, label='all sources')\n",
    "plt.plot(mags_g['psfmag'][is_primary & deblended & starsel_i & sel_cen][:,0]-mags_i['psfmag'][is_primary & deblended & starsel_i & sel_cen][:,0],\\\n",
    "         mags_g['psfmag'][is_primary & deblended & starsel_i & sel_cen][:,0], 'ro', ms=5, alpha=0.75, label='point sources')\n",
    "plt.plot(g_iso[okstage]-i_iso[okstage], g_iso[okstage], 'k', label='__no label__')#s', ms=1)\n",
    "plt.xlabel('g-i')\n",
    "plt.ylabel('g')\n",
    "plt.legend()\n",
    "plt.xlim(-1.3,2.8)\n",
    "plt.ylim(27.8,22.3)\n",
    "plt.title('within 2.5 arcmin')\n",
    "plt.minorticks_on()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "In this notebook, we explored the photometry around a dwarf galaxy in the HSC dataset, with a particular interest in extracting scientifically useful photometry for stars in this nearby, partially-resolved galaxy.\n",
    "\n",
    "Among other things, we did the following: \n",
    "\n",
    "* Identified the tract and patch containing an object of known coordinates.\n",
    "\n",
    "* Extracted cutout images and catalogs centered on the identified position, including data from multiple patches if needed. There is a beautiful dwarf galaxy -- KKR 25 -- in the images, which we explored in more detail.\n",
    "\n",
    "* Plotted a color-magnitude diagram centered on the object of interest. In this particular case, the red giant branch of dwarf galaxy KKR 25 is not evident in the CMD, even though it looks like it _should_ be.\n",
    "\n",
    "* Used a variety of measured quantities from the source catalogs to (a) determine that the stars in KKR 25 are being detected, but just have poor measurements, and (b) to explore why their measurements are all consistent with having zero flux.\n",
    "\n",
    "* Following the Stack Club \"Footprints\" notebook in the Source Detection section: https://stackclub.readthedocs.io/en/latest/notebooks.html#source-detection, we looked at source footprints within the KKR 25 region and determined that the issue is that some footprints are _huge_, and that the default behavior of the measurement algorithm(s) to replace neighboring footprints with noise was causing all of the sources to get replaced with noise, too. \n",
    "\n",
    "* Reprocessed a small chunk of data with the flag `doReplaceWithNoise` set to `False`. This allows sources to be measured even in the crowded regions of KKR 25.\n",
    "\n",
    "* We checked the CMD of this newly-processed dataset, and found that the red giant branch of KKR 25 is clearly visible, and the photometry looks good! From here, we could proceed with our scientific exploration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.8"
  },
  "livereveal": {
   "scroll": true,
   "start_slideshow_at": "selected"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}