surface_photometry.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Doing quick surface photometry of a galaxy given its coordinates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Owners: **Aaron Watkins, Nushkia Chamba**  \n",
    "Last Verified to Run: **14 August 2020**  \n",
    "Verified Stack release: **20.0.0+7**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Objectives: \n",
    "\n",
    "This worksheet is a demonstration of how to do the following:\n",
    "1. Locate a specific object from available databases in the Stack using its right ascension and declination coordinates\n",
    "2. Obtain a thumbnail of the object for image analysis\n",
    "3. Create nice RGB images from multiband exposure\n",
    "4. Mask unwanted sources surrounding and overlapping the object using source detection footprints and Holoviews\n",
    "5. Do quick surface photometry of the object on the masked image, with a local background estimation\n",
    "\n",
    "The basic goal is to provide a demonstration of one way to perform one's own image analysis within the Stack, on present and future imaging data releases.  It can be used as a stepping off point for development of one's own software.\n",
    "\n",
    "### Acknowledgements \n",
    "\n",
    "We would like to thank Alex Drlica-Wagner and Greg Madejski for giving us the opportunity to participate in the online (8 week) LSST Stack Club Course 2020, which went ahead in May and June during the peak of the COVID-19 lockdown. \n",
    "\n",
    "Special thanks to Alex (again!), Yusra, Remy, Keith, Fred and Cristina for providing help and advice while we put this notebook together. We also extend our gratitude to the rest of the Stack Club team for welcoming us warmly into the weekly Stack Club sessions and for building a friendly community to work with :) \n",
    "\n",
    "If you use the techniques presented here, please cite:\n",
    "[Watkins et al. (2016)](https://ui.adsabs.harvard.edu/abs/2016ApJ...826...59W/abstract) and\n",
    "[Chamba et al. (2020)](https://ui.adsabs.harvard.edu/abs/2020A%26A...633L...3C/abstract)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **Initial Setup**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# What version of the Stack am I using?\n",
    "! echo $HOSTNAME\n",
    "! eups list lsst_distrib -s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **Packages and dependencies**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Image retrieval tasks\n",
    "from lsst.obs.hsc.hscMapper import HscMapper\n",
    "import lsst.daf.persistence as dafPersist\n",
    "from lsst.daf.persistence import Butler\n",
    "import lsst.afw.table as afwTable\n",
    "import lsst.afw.display as afwDisplay\n",
    "import lsst.afw.detection as afwDetection\n",
    "import lsst.geom as geom\n",
    "from lsst.geom import Box2D\n",
    "from lsst.afw.image import Exposure, Image, PARENT, MultibandExposure, MultibandImage\n",
    "\n",
    "# Detection and deblending tasks\n",
    "from lsst.meas.algorithms import SourceDetectionTask\n",
    "from lsst.meas.base import SingleFrameMeasurementTask\n",
    "from lsst.afw.table import SourceCatalog\n",
    "import scarlet\n",
    "import scarlet.display\n",
    "import lsst.meas.extensions.scarlet as mes\n",
    "from scarlet_extensions.initialization import detection as det\n",
    "\n",
    "# All other (non-lsst) tasks\n",
    "import os\n",
    "import sys\n",
    "import glob\n",
    "from copy import copy\n",
    "import warnings\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.colors as colors\n",
    "from matplotlib.patches import Ellipse\n",
    "from IPython.display import display, Markdown\n",
    "import numpy as np\n",
    "from astropy.table import Table\n",
    "import astropy.coordinates as coord\n",
    "from astropy.visualization import make_lupton_rgb\n",
    "\n",
    "afwDisplay.setDefaultBackend('matplotlib')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One final dependency is the following:\n",
    "https://github.com/herjy/scarlet_extensions\n",
    "\n",
    "Git clone this, then move the scarlet_extensions subdirectory into this notebook's directory to perform the deblending aspects of this notebook.\n",
    "\n",
    "This package requires SEP (SExtractor for Python).  In an LSST terminal, do `pip install sep`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **Functions**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will define some functions, in advance, to do the actual surface photometry.  We will be focusing on galaxies, hence require a way to transform the thumbnail images into in-plane galaxy coordinates and a way to measure the light profiles in those coordinates.\n",
    "\n",
    "To create the RGB image itself, we use of the `make_lupton_rb()` function in Astropy. However, the RGB image from this function consistently makes galaxies more red than their global colors suggest. Therefore, in this notebook we manually scale the individual rgb channels prior to input and thus improve the quality of the resulting RGB image. \n",
    "We add a 'scale' parameter to do this in the following function called `createRGB()`.\n",
    "\n",
    "Inspiration for exploring a 'scale' parameter stemmed from Lee Kelvin's [plotfits](https://rdrr.io/cran/astro/man/plotfits.html) code in the R programming language (Thanks Javier Román for pointing NC to this reference!). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mk_ell(imdata, xcen, ycen, pa, ell):    \n",
    "    '''\n",
    "    For reprojecting into elliptical coordinates.  \n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    imdata: ndarray \n",
    "            Image data array\n",
    "    xcen: float \n",
    "            Center x coordinate (pixels)\n",
    "    ycen: float \n",
    "            Center y coordinate (pixels)\n",
    "    pa: float   \n",
    "            Position angle (degrees)\n",
    "    ell: float   \n",
    "            Ellipticity (1-axial ratio)\n",
    "            \n",
    "    Returns\n",
    "    -------\n",
    "    ellipse: dictionary object\n",
    "            {'xEll': x coordinates in inclined galaxy plane\n",
    "            'yEll': y coordinates in inclined galaxy plane\n",
    "            'ellRad': radius in inclined galaxy plane\n",
    "            'theta': angle in galaxy plane, starting 0 at PA and increasing CCW (radians)\n",
    "            }\n",
    "            \n",
    "    For example, to plot the deprojected image, do e.g.:\n",
    "        plt.pcolormesh(xEll, yEll, image)\n",
    "        plt.imshow(ellRad, origin='bottom')\n",
    "    '''\n",
    "    dimx = imdata.shape[1]\n",
    "    dimy = imdata.shape[0]\n",
    "    if (dimy > dimx):\n",
    "        x = np.arange(1, dimy+1)\n",
    "        y = x.reshape(-1, 1)\n",
    "    else:\n",
    "        x = np.arange(1, dimx+1)\n",
    "        y = x.reshape(-1, 1)\n",
    "    # Currently written to start at PA and increase CCW (due north is 0 degrees)\n",
    "    xEll = (x - xcen)*np.cos(np.radians(pa+90)) + \\\n",
    "                (y - ycen)*np.sin(np.radians(pa+90))\n",
    "    yEll = -(x - xcen)*np.sin(np.radians(pa+90)) + \\\n",
    "                 (y - ycen)*np.cos(np.radians(pa+90))\n",
    "\n",
    "    ellRad = np.sqrt(xEll**2 + (yEll/(1-ell))**2)\n",
    "    theta = np.arctan2(yEll, xEll)\n",
    "    theta[theta < 0.] = theta[theta < 0.]+2.0*np.pi # Removing negative values\n",
    "    \n",
    "    ellipse = {'xEll': xEll,\n",
    "               'yEll': yEll,\n",
    "               'ellRad': ellRad,\n",
    "               'theta': theta}\n",
    "    \n",
    "    return ellipse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def surfbprof(imdata, maskdata, cen, pa, ell, bwidth, maxbin, iflag):\n",
    "    '''\n",
    "    Quickly does isophotal profiles of galaxy using image sections.  \n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    imdata: ndarray \n",
    "            Image data array\n",
    "    maskdata: ndarray (same shape as imdata)\n",
    "            Mask data array (0 = not masked, > 0 is masked)*\n",
    "    cen: geom.PointI object\n",
    "            Central coordinates\n",
    "    pa: float   \n",
    "            Position angle (degrees)\n",
    "    ell: float   \n",
    "            Ellipticity (1-axial ratio)\n",
    "    bwidth: int \n",
    "            Bin width (must be pixels--works off image array)\n",
    "    maxbin: int\n",
    "            Maximum radius (must be pixels--works off image array)\n",
    "    iflag: int \n",
    "            Flag for whether or not to correct for inclination (mu[i]*(1-ell))\n",
    "            -- 0 = don't correct, != 0 means correct\n",
    "            \n",
    "    Returns \n",
    "    -------\n",
    "    photo_dict: dictionary object \n",
    "            {'rad': array\n",
    "                    Radius array (in pixels),\n",
    "            'mean_I': array (len(rad))\n",
    "                    Mean surface brightness within elliptical annuli\n",
    "            'med_I': array (len(rad))\n",
    "                    Median surface brightness within elliptical annuli\n",
    "            }\n",
    "\n",
    "    *NOTE: an alternative is to simply pass this function an array with NaN values\n",
    "    where the masked pixels should be.  We do this in this notebook, but having a \n",
    "    mask array separate may be more convenient for some users, so we leave this \n",
    "    function unaltered.\n",
    "    '''\n",
    "    xcen_i = int(np.round(cen[0], 0)) # Converting to integer value\n",
    "    ycen_i = int(np.round(cen[1], 0)) # Converting to integer value\n",
    "    \n",
    "    if len(np.array([pa]))==1:\n",
    "        ellipse = mk_ell(imdata, xcen_i, ycen_i, pa, ell)\n",
    "        ellRad = ellipse['ellRad']\n",
    "        ellRad = ellRad[ : imdata.shape[0], : imdata.shape[1]]\n",
    "    \n",
    "    bins = np.arange(0.5*bwidth, maxbin+0.5*bwidth, bwidth)\n",
    "    phot_dict = {}\n",
    "    mean_i = np.zeros(len(bins))-999\n",
    "    med_i = np.zeros(len(bins))-999\n",
    "    for i in range(len(bins)):\n",
    "        # First shrinking image array to only needed size\n",
    "        minbin = bins[i]-0.5*bwidth\n",
    "        maxbin = bins[i]+0.5*bwidth\n",
    "        offset = int(np.round(maxbin+5)) # how far to grow the image box\n",
    "        imsec = imdata[ycen_i-offset : ycen_i+offset,\n",
    "                       xcen_i-offset : xcen_i+offset]\n",
    "        masksec = maskdata[ycen_i-offset : ycen_i+offset,\n",
    "                           xcen_i-offset : xcen_i+offset]\n",
    "        \n",
    "        if len(np.array([pa]))==1:\n",
    "            ellRadsec = ellRad[ycen_i-offset : ycen_i+offset,\n",
    "                               xcen_i-offset : xcen_i+offset]\n",
    "            want_sb = (ellRadsec >= minbin) & (ellRadsec < maxbin) & (masksec == 0)\n",
    "        else:\n",
    "            ellipse = mk_ell(imsec, offset, offset, pa[i], ell[i])\n",
    "            ellRad = ellipse['ellRad']\n",
    "            want_sb = (ellRad >= minbin) & (ellRad < maxbin) & (masksec == 0)\n",
    "            \n",
    "        # Ignore empty bins\n",
    "        if len(want_sb[want_sb]) == 0:\n",
    "            mean_i[i] = -999\n",
    "            med_i[i] = -999\n",
    "        else:\n",
    "            if iflag:\n",
    "                mean_i[i] = np.nanmean(imsec[want_sb]) * (1-ell)\n",
    "                med_i[i] = np.nanmedian(imsec[want_sb]) * (1-ell)\n",
    "            else:\n",
    "                mean_i[i] = np.nanmean(imsec[want_sb])\n",
    "                med_i[i] = np.nanmedian(imsec[want_sb])\n",
    "    \n",
    "    phot_dict['rad'] = bins\n",
    "    phot_dict['mean_I'] = mean_i\n",
    "    phot_dict['med_I'] = med_i\n",
    "    \n",
    "    return phot_dict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def createRGB(image, bgr=\"gri\", stretch=1, Q=10, scale=None):\n",
    "    \"\"\"\n",
    "    Create an RGB color composite image.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    image : `MultibandImage`\n",
    "        `MultibandImage` to display.\n",
    "    bgr : sequence\n",
    "        A 3-element sequence of filter names (i.e. keys of the exps dict) indicating what band\n",
    "        to use for each channel. If `image` only has three filters then this parameter is ignored\n",
    "        and the filters in the image are used.\n",
    "    stretch: int\n",
    "        The linear stretch of the image.\n",
    "    Q: int\n",
    "        The Asinh softening parameter.\n",
    "    scale: list of 3 floats, each less than 1. (default: None)\n",
    "        Re-scales rgb channels  \n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    rgb: ndarray \n",
    "        RGB (integer, 8-bits per channel) colour image as an NxNx3 numpy array.\n",
    "    \"\"\"\n",
    "    \n",
    "    # If the image only has 3 bands, reverse the order of the bands to produce the RGB image\n",
    "    if len(image) == 3:\n",
    "        bgr = image.filters\n",
    "    # Extract the primary image component of each Exposure with the .image property, and use \n",
    "    # .array to get a NumPy array view.\n",
    "\n",
    "    if scale is None:\n",
    "        \n",
    "        r_im = image[bgr[2]].array # numpy array for the r channel\n",
    "        g_im = image[bgr[1]].array # numpy array for the g channel\n",
    "        b_im = image[bgr[0]].array # numpy array for the b channel\n",
    "    \n",
    "    else:\n",
    "        ### manually re-scaling the images here \n",
    "        r_im = image[bgr[2]].array*scale[0]\n",
    "        g_im = image[bgr[1]].array*scale[1]\n",
    "        b_im = image[bgr[0]].array*scale[2]     \n",
    "    \n",
    "    rgb = make_lupton_rgb(image_r=r_im,  \n",
    "                          image_g=g_im,  \n",
    "                          image_b=b_im,  \n",
    "                          stretch=stretch, Q=Q)  # parameters used to stretch and scale the pixel values\n",
    "    \n",
    "    return rgb          "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **Part 1: Selecting a galaxy by coordinates**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As detailed in other tutorials, first we need to set up a Butler based on the parameters we obtained in the previous section.  We do this by galaxy coordinates, looping through available imaging repositories (DEEP which includes DEEP and ULTRA-DEEP tracts and WIDE). We selected this repository as it is the most complete dataset of HSC PDR2 in the LSP at the time of writing this notebook [(more info)]( https://confluence.lsstcorp.org/display/DM/S20+HSC+PDR2+Reprocessing). Thank you Yusra AlSayaad for the much needed help with this!\n",
    "\n",
    "NOTE: currently this is built for HSC imaging data only.  Also, the current notebook requires Butler gen. 2; this will malfunction with Butler gen. 3 due to the HSC repository we are accessing here (/rerun/DM-23243).  For gen. 3, use /rerun/RC/w_2020_26/DM-*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The coordinates of a nice barred spiral galaxy in COSMOS\n",
    "ra = 150.746062\n",
    "dec = 2.343061"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pos = geom.SpherePoint(ra, dec, geom.degrees) # Using geom package for the galaxy coordinates\n",
    "hsc_repos = '/datasets/hsc/repo/rerun/' # Access HSC PDR2 re-processing repositories\n",
    "repo_dir_list = glob.glob(os.path.join(hsc_repos, 'DM-23243', 'OBJECT', '*'))\n",
    "print('Search repos...', repo_dir_list)\n",
    "\n",
    "for repo in repo_dir_list:\n",
    "    # First set up the Butler for the data repositories\n",
    "    butler = Butler(repo)\n",
    "    # The skyMap is used to locate image tracts and patches\n",
    "    skyMap = butler.get('deepCoadd_skyMap')    \n",
    "    tractList = skyMap.findAllTracts(pos)\n",
    "    # Only accept tracts that contain the galaxy coordinates\n",
    "    if not tractList:\n",
    "        continue\n",
    "    else:\n",
    "        tractInfo = tractList[0]\n",
    "        patchInfo = tractInfo.findPatch(pos)    \n",
    "        tract = tractInfo.getId() # This is now the tract number we want\n",
    "        patch = patchInfo.getIndex() # And this is the patch in the tract that we want\n",
    "        patchid = str(patch[0])+','+str(patch[1])\n",
    "        print('which repo? ', repo)\n",
    "        print('tract:', tract, ' patchID:', patchid)\n",
    "        break"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **Part 2: Loading in the galaxy thumbnails for analysis**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Prior to analysis, we first need to load in all of the relevant images.  In this case, we will get thumbnails for all five HSC photometric bands and place them into a dictionary, accessible by filter name."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NOTE: regarding photometric bands, the repository this notebook is using requires the photometric bands R2 and I2 (which are new filters that were used during HSC PDR2).  If you choose to use a different HSC repository and it throws an error, try switching these to filters 'HSC-I' and 'HSC-R' as in the commented out line below.  When running this notebook using Butler gen. 3 (see above), use 'HSC-I', and 'HSC-R' (if not from PDR2)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bands = ['HSC-G', 'HSC-R2', 'HSC-I2', 'HSC-Z', 'HSC-Y']\n",
    "#bands = ['HSC-G', 'HSC-R', 'HSC-I', 'HSC-Z', 'HSC-Y']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will load in the images by taking cutouts.  This is done by establishing a bounding box of a given size, then using the Butler to access the 'deepCoadd_sub' object for the given tract, patch, and filters.  For this demonstration the bounding box size is hard-coded; you can use detection catalogues to estimate this size for a more general-purpose project."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bwidth = 950\n",
    "cutout_size = geom.ExtentI(bwidth, bwidth)\n",
    "calexps = {}\n",
    "# First obtain position information regarding the bounding box. We'll use only the first band for convenience.\n",
    "wcs = butler.get(\"deepCoadd_wcs\", immediate=True,\n",
    "                 tract=tract, patch=patchid, filter=bands[0])\n",
    "xy = geom.PointI(wcs.skyToPixel(pos))\n",
    "bbox = geom.BoxI(xy - cutout_size//2, cutout_size)\n",
    "    \n",
    "# Now assigning the exposure objects into a dictionary, accessed by photometric band.\n",
    "for band in bands:\n",
    "    calexp = butler.get(\"deepCoadd_sub\",\n",
    "                        tract = tract, \n",
    "                        patch = patchid, \n",
    "                        filter = band,\n",
    "                        bbox = bbox, \n",
    "                        immediate = True)\n",
    "    calexps[band] = calexp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now use the methods attached to these as usual.  For example, we can grab the photometric zeropoint for our later surface photometry, and check that the photometric zeropoints are the same for all filters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "magZp = calexps['HSC-I2'].getPhotoCalib().instFluxToMagnitude(1)\n",
    "for band in bands:\n",
    "    print('ZP for '+band+': ',calexps[band].getPhotoCalib().instFluxToMagnitude(1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or we can access the image arrays and check that we found our target galaxy.\n",
    "\n",
    "(Here we are using the [afwDisplay package](https://github.com/LSSTScienceCollaborations/StackClub/blob/master/Visualization/AFW_Display_Demo.ipynb).)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,5, figsize=(15, 5))\n",
    "display = []\n",
    "for i in range(len(bands)):\n",
    "    plt.sca(ax[i])\n",
    "    afw_display = afwDisplay.Display(frame=fig)\n",
    "    display.append(afw_display)\n",
    "    display[i].scale(\"asinh\", \"zscale\")\n",
    "    display[i].setMaskTransparency(100)\n",
    "    display[i].mtv(calexps[bands[i]])\n",
    "    ax[i].set_title(bands[i])\n",
    "    ax[i].set_xticks([])\n",
    "    ax[i].set_yticks([])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **Making RGB image**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can make a nice RGB image of the galaxy! This is not only for aesthetic reasons, but you will see later that RGB images are very useful to visualise a segmentation map or mask to verify the detected sources in the image. It can also help in the analysis process to manually mask regions of the image (if needed), directly identify which regions of the mask belong to the galaxy, and also make a correspondance between the galaxy image and its structure in the surface brightness profile. \n",
    "\n",
    "In this notebook, we will compare the original RGB output from the default `make_lupton_rgb()` function and the scaled RGB version explored here. This function requires specifically a multiband exposure object to work with, so first we construct one of those."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Multiband exposures need a list of images and filters\n",
    "# Using only the first three bands\n",
    "coadds = [calexps[bnd] for bnd in bands[:3]]\n",
    "coadds = MultibandExposure.fromExposures(bands[:3], coadds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(20, 20), nrows=1, ncols=2)\n",
    "\n",
    "##original make_lupton_rgb without any scaling \n",
    "rgb_original = createRGB(coadds.image, bgr=bands, scale=None)\n",
    "ax[0].imshow(rgb_original, origin='lower')\n",
    "ax[0].set_title('ORIGINAL', fontsize=30)\n",
    "\n",
    "##make_lupton_rgb with scaled rgb channels (this notebook) \n",
    "ax[1].set_title('RE-SCALED (this notebook)', fontsize=30)\n",
    "rgb_scaled = createRGB(coadds.image, bgr=bands, scale=[0.5, 0.7, 1.0])\n",
    "ax[1].imshow(rgb_scaled, origin='lower')\n",
    "\n",
    "ax[0].set_axis_off()\n",
    "ax[1].set_axis_off()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In comparison, the re-scaled RGB image produces a visual representation of the spiral galaxy that is closer to the expected, blue global colour of the galaxy. Therefore, here onwards we use the re-scaled RGB image to visualise our results. \n",
    "\n",
    "Now we can prepare the images for photometry. First we need to elliptical parameters of the galaxy.  We'll do this visually, so our initial guess is position angle of 150 degrees, and axial ratio of 0.9:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We will use the I2 band, which looks comparitively deep\n",
    "image_I2 = calexps['HSC-I2'].image.array\n",
    "xcen = image_I2.shape[1]//2\n",
    "ycen = image_I2.shape[0]//2\n",
    "cen = (xcen, ycen)\n",
    "PA = 150\n",
    "q = 0.9\n",
    "size = 450 \n",
    "\n",
    "ellipse = Ellipse(xy=cen, width=size, height=size*q, angle= PA,\\\n",
    "                  edgecolor='deeppink', fc='None', lw=3, zorder = 2)\n",
    "\n",
    "# Here we add a luminance channel to view the outskirts of the object in grey scale more clearly \n",
    "rgb_alpha = np.stack([rgb_scaled[:,:,0], \n",
    "                      rgb_scaled[:,:,1], \n",
    "                      rgb_scaled[:,:,2],\n",
    "                      np.sum(rgb_scaled, axis=2)*10], \n",
    "                     axis=2)\n",
    "\n",
    "plt.figure(figsize=(10,10))\n",
    "plt.imshow(rgb_alpha, origin='bottom')\n",
    "plt.gca().add_patch(ellipse)\n",
    "\n",
    "plt.axis('off')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **Part 3: Masking unwanted sources**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we will set up our deblender and detection tasks.  This is using the [lsst_stack_deblender tutorial](https://github.com/LSSTScienceCollaborations/StackClub/blob/master/Deblending/lsst_stack_deblender.ipynb) as a baseline, so see that page for more detailed information.  We'll use the I2-band for our detection, which in our images above (in gray-scale) looks to have a deep sensitivity, and falls in a red wavelength regime that should help pick out more high-redshift objects as well.  The goal of detection is to make an initial, harsh mask to block out the galaxy and any external objects.  The goal of the deblending will be to pick out objects embedded within the galaxy for more careful photometry (that excludes blended sources)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **Detection**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "schema = SourceCatalog.Table.makeMinimalSchema() # Set up an empty schema\n",
    "\n",
    "# The source detection\n",
    "detectionConfig = SourceDetectionTask.ConfigClass()\n",
    "detectionConfig.thresholdValue = 1  # This is our detection threshold\n",
    "detectionConfig.thresholdType = \"stdev\" # Specifying we want it in standard deviations\n",
    "detectionTask = SourceDetectionTask(schema=schema) # Set up a detection task for the previous schema"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll run the detection task on the whole cutout to create the initial mask (galaxy + external objects).  Again, we'll use the I2-band given its high quality.  Note that detection works using a single photometric band, whereas the Scarlet deblender (which we make use of later) uses multiband images."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "detectionTable = SourceCatalog.Table.make(schema)\n",
    "detectionResult = detectionTask.run(detectionTable, calexps['HSC-I2']) # Using only I-band\n",
    "detectionCatalog = detectionResult.sources"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll use the coordinates of the calexp bounding box, along with the indexes and values of the footprints, to pinpoint the locations of the footprints on the image thumbnails.  The bounding box mininum is the lower-left corner of the thumbnail (indexes 0,0), so we just subtract that value from the footprint indices to find our coordinates in thumbnail image space.\n",
    "\n",
    "We'll also identify and then separate out the footprint that contains the galaxy we're interested in.  We do want to preserve this mask for a later step, that being the local sky background estimation (Part 4).  Finally, even at 1 sigma sensitivity, the detection task creates footprints that leave a lot of diffuse light leaking out through the resulting masks, so we use the `.dilate()` function contained within them to improve the masks on the fly.  We increase the galaxy mask size substantially to aid in the sky estimation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xy_begin = bbox.getMin()\n",
    "mask = np.zeros(image_I2.shape) # Masking all detections, except galaxy of interest\n",
    "sky_mask = np.zeros(image_I2.shape) # Mask all detections, including galaxy of interest, for sky estimation\n",
    "atol = 50 # Requires a high number: the footprint centroid is not the galaxy centroid\n",
    "\n",
    "for j in range(len(detectionCatalog['id'])):\n",
    "    centroid = detectionCatalog[j].getFootprint().getCentroid()\n",
    "    \n",
    "    # For identifying the target galaxy; is the footprint centroid within atol of the galaxy center\n",
    "    want1 = np.isclose(centroid[0], xy_begin[0]+bwidth//2, atol=atol)\n",
    "    want2 = np.isclose(centroid[1], xy_begin[1]+bwidth//2, atol=atol)\n",
    "    \n",
    "    #yind, xind = detectionCatalog[j].getFootprint().spans.indices()\n",
    "    fp = detectionCatalog[j].getFootprint()\n",
    "    if want1 & want2:\n",
    "        fp.dilate(75) # Increasing the size of the galaxy mask substantially\n",
    "    else:\n",
    "        fp.dilate(5) # Increasing other masks, but not as much\n",
    "    yind, xind = fp.spans.indices()\n",
    "    new_yind = np.array(yind, dtype=int) - xy_begin[1] # yind, xind are lists, so we need to convert to do math\n",
    "    new_xind = np.array(xind, dtype=int) - xy_begin[0] \n",
    "    # After dilation, some masks exceed the image boundary\n",
    "    in_im_y = (new_yind > 0) & (new_yind < 950)\n",
    "    in_im_x = (new_xind > 0) & (new_xind < 950)\n",
    "    sky_mask[new_yind[in_im_y & in_im_x], new_xind[in_im_y & in_im_x]] = np.nan\n",
    "    \n",
    "    if want1 & want2:\n",
    "        obj_ind = j\n",
    "        continue\n",
    "        \n",
    "    # Now setting all mask pixels to NaN\n",
    "    # Alternative with our software is to assign values > 0 to a blank image and use as mask array\n",
    "    mask[new_yind[in_im_y & in_im_x], new_xind[in_im_y & in_im_x]] = np.nan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Checking that our resulting mask makes sense\n",
    "fig, ax = plt.subplots(figsize=(21, 7), nrows=1, ncols=3)\n",
    "ax[0].imshow(rgb_alpha, origin='bottom')\n",
    "ax[0].set_title('RGBA IMAGE', fontsize=25)\n",
    "ax[0].set_xlabel('X [pixels]', fontsize=20)\n",
    "ax[0].set_ylabel('Y [pixels]', fontsize=20)\n",
    "\n",
    "# Copy the RGB image \n",
    "rgb_alpha_copy = copy(rgb_alpha)\n",
    "\n",
    "for k in range(0, 4):\n",
    "    rgb_alpha_copy[:,:,k] = rgb_alpha_copy[:,:,k] + mask \n",
    "\n",
    "ax[1].imshow(rgb_alpha_copy,  origin='bottom')\n",
    "ax[1].set_title('MASK', fontsize=25)\n",
    "ax[1].set_xlabel('X [pixels]', fontsize=20)\n",
    "\n",
    "for k in range(0, 4):\n",
    "    rgb_alpha_copy[:,:,k] = rgb_alpha_copy[:,:,k] + sky_mask\n",
    "\n",
    "ax[2].imshow(rgb_alpha_copy,  origin='bottom')\n",
    "ax[2].set_title('SKY MASK', fontsize=25)\n",
    "ax[2].set_xlabel('X [pixels]', fontsize=20)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This shows that we can pick out pixels that belong to everything except the central galaxy and whatever surrounding galaxies are part of the same footprint (right panel) or that we can include the target galaxy mask as well (for background sky estimation, shown in the next section)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **Deblending**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course, for a more correct mask we also want to include masks around sources overlapping the galaxy.  Because the current implementation of the deblender in the Stack is out of date, we're going to use the following package by Remy Joseph:\n",
    "\n",
    "https://github.com/herjy/scarlet_extensions\n",
    "\n",
    "To run this notebook from this point onward, you will need to do a git clone on this package in your working directory, then move the scarlet_extensions folder out from the newly created repository.  Also, this package requires SEP, a SExtractor protocol for Python.  This can be installed quite easily: `pip install sep` (however, make sure you do it in an LSST environment)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Scarlet deblending requires specifically a multiband exposure object to work with. We already created this for the `createRGB()` function, using the first three filters only."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "images = coadds.image.array\n",
    "seg_catalog, bg_rms = det.makeCatalog(images)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The object `seg_catalog` is a SExtactor table (use `seg_catalog.dtype` to see the names of the table elements), including position, position angle, and ellipticity estimates for the deblended sources within the target galaxy. So we can easily check the output of this segmentation by plotting the output ellipses on the galaxy image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(12,12))\n",
    "\n",
    "# Re-scaled color image from before \n",
    "im = ax.imshow(rgb_alpha, origin='lower')\n",
    "\n",
    "# Array to assign unique ID for deblended source \n",
    "objID = np.arange(start=0, stop=len(seg_catalog), dtype=np.int)\n",
    "ellipsemask_wIDs = np.zeros(image_I2.shape)\n",
    "\n",
    "#axis-ratio of objects\n",
    "q_objs = seg_catalog['b']/seg_catalog['a']\n",
    "\n",
    "#Choice here: objects with axis-ratio lower than 0.4, plot with smaller ellipse \n",
    "factors = np.ones(len(seg_catalog))*10\n",
    "factors[q_objs < 0.4] = 6\n",
    "\n",
    "# Then plot an ellipse for each deblender detection\n",
    "# Colors are randomized to aid in visibility\n",
    "for i in range(len(seg_catalog)):\n",
    "    \n",
    "    q_obj = q_objs[i]\n",
    "    \n",
    "    e = Ellipse(xy=(seg_catalog['x'][i], seg_catalog['y'][i]),\n",
    "                width=2.*factors[i]*seg_catalog['a'][i],\n",
    "                height=2.*factors[i]*seg_catalog['b'][i],\n",
    "                angle=seg_catalog['theta'][i] * 180. / np.pi, \n",
    "                lw=2)\n",
    "    \n",
    "    # Create ellipse around each source to assign objID\n",
    "    ellipse_dict = mk_ell(image_I2, \n",
    "                          seg_catalog['x'][i], \n",
    "                          seg_catalog['y'][i], \n",
    "                          seg_catalog['theta'][i] * 180. / np.pi - 90., \n",
    "                          1.-q_obj)\n",
    "\n",
    "    ellipsemask_wIDs[(ellipse_dict['ellRad'] < factors[i]*seg_catalog['a'][i])] = objID[i]\n",
    "    \n",
    "    color = list(np.random.uniform(size=3))\n",
    "    color.append(0.5)\n",
    "    e.set_facecolor(color)\n",
    "    e.set_edgecolor(color)\n",
    "    ax.add_artist(e)\n",
    "\n",
    "ax.set_xlabel('X [pixels]', fontsize=20)\n",
    "ax.set_ylabel('Y [pixels]', fontsize=20)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **Manual masking with Holoviews**\n",
    "\n",
    "We take advantage of Holoviews to (manually and interactively) select deblending sources to be later masked/unmasked. The following function was directly taken from the StackClub notebook by Keith Bechtol: Visualization/bokeh_holoviews_datashader.ipynb [github link](https://github.com/LSSTScienceCollaborations/StackClub/blob/project/bokeh_holoviews_datashader/bechtol/Visualization/bokeh_holoviews_datashader.ipynb) with a few additions to plot the ellipses of the detected objects above.\n",
    "\n",
    "We import the relevant packages below, in case the user wants to skip this step for whatever reason. (Though we do recommend careful manual adjustment of/additions to the default output.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import holoviews as hv\n",
    "from holoviews import streams\n",
    "from holoviews.operation.datashader import datashade, dynspread, rasterize\n",
    "hv.extension('bokeh')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In our case, we need to unmask the deblended central regions of the galaxy as well as expand the mask of the object at (250, 700) that has not been completely masked with the fixed ellipse size selected above. For this reason, we create two Holoviews plots so that the user may run seperate cells for masking and unmasking regions.\n",
    "\n",
    "In this step, we encourage the user to spend some time playing/exploring the manual masking process.  It is easy to redo anything done in the following two steps. You can reset all your selections by running the first cell below at any time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Array with all deblended objects detected as NaN\n",
    "# Re-run this cell if you want to reset all your masking selections\n",
    "mask_objects = np.zeros(ellipsemask_wIDs.shape)\n",
    "mask_objects[ellipsemask_wIDs > 0.] = np.nan\n",
    "\n",
    "# Initialize master mask (automated mask + manual mask)\n",
    "master_mask = np.zeros(mask.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Unmasking regions with Holoviews\n",
    "\n",
    "Using the left panel of the following plot, you may draw one region to select points (corresponding to the object masks above) that you want to unmask.\n",
    "\n",
    "NOTE: with this method, you can only select one region at a time.  See below for how to add multiple selections to the final mask."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%output size=150\n",
    "%%opts Points [tools=['box_select', 'lasso_select']]\n",
    "\n",
    "# Declare catalog points \n",
    "points = hv.Points((seg_catalog['x'], seg_catalog['y']))\n",
    "points.opts(color='cyan', size=5)\n",
    "\n",
    "# Draw the ellipses of the sources \n",
    "clusters = hv.Polygons([hv.Ellipse(seg_catalog['x'][i], seg_catalog['y'][i], (2.*factors[i]*seg_catalog['a'][i], 2.*factors[i]*seg_catalog['b'][i]), orientation= seg_catalog['theta'][i]) for i in range(0, len(seg_catalog))])\n",
    "clusters.opts(color='grey')\n",
    "\n",
    "# Draw the ellipse enclosing the galaxy\n",
    "galaxy_ellipse = hv.Ellipse(xcen, ycen, spec=(2*250, 2*250*q), orientation= np.radians(90. - PA)) \n",
    "galaxy_ellipse.opts(color='deeppink', line_width=3)\n",
    " \n",
    "# Declare points as source of selection stream\n",
    "selection_unmask = streams.Selection1D(source=points)\n",
    "\n",
    "# Write function that uses the selection indices to slice points and compute stats\n",
    "def selected_info(index):\n",
    "    selected = points.iloc[index]\n",
    "    if index:\n",
    "        label = 'Mean x, y: %.3f, %.3f' % tuple(selected.array().mean(axis=0))\n",
    "    else:\n",
    "        label = 'No selection'\n",
    "    return selected.relabel(label).options(color='red')\n",
    "\n",
    "# Combine points and DynamicMap\n",
    "clusters*galaxy_ellipse*points + hv.DynamicMap(selected_info, streams=[selection_unmask])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To UNMASK a selected region, you can run the cell below. Afterwards, you can make another selection above and repeat this step to unmask other regions if you wish."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "unmask_objID_selected = np.array(selection_unmask.index)\n",
    "\n",
    "# If selected index is objID then unmask\n",
    "\n",
    "for obj_sel in unmask_objID_selected:\n",
    "    mask_objects[ellipsemask_wIDs == obj_sel] = 0. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's overplot this new mask on the image to see if we have unmasked what we intended."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8,8))\n",
    "\n",
    "rgb_alpha_masked = copy(rgb_alpha)\n",
    "\n",
    "master_mask = mask_objects + mask # This is our updated mask with the unmasked region \n",
    "\n",
    "for k in range(0, 4):\n",
    "    rgb_alpha_masked[:,:,k] = rgb_alpha_masked[:,:,k] + master_mask\n",
    "    \n",
    "im = ax.imshow(rgb_alpha_masked, origin='bottom')\n",
    "\n",
    "ax.set_xlabel('X [pixels]', fontsize=20)\n",
    "ax.set_ylabel('Y [pixels]', fontsize=20)\n",
    "plt.title('RGB + UNMASKED region')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Masking regions with Holoview\n",
    "\n",
    "Using the left panel of the following plot, you may draw a region to select points (corresponding to the object masks above) that you want to mask, in the same manner as the previous unmasking step. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%output size=150\n",
    "%%opts Points [tools=['box_select', 'lasso_select']]\n",
    "\n",
    "# Declare points as source of selection stream\n",
    "selection_mask = streams.Selection1D(source=points)\n",
    "\n",
    "# Combine points and DynamicMap\n",
    "clusters*galaxy_ellipse*points + hv.DynamicMap(selected_info, streams=[selection_mask])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To MASK the selected region, you may run the cell below. For this particular example, both the UNMASK and MASK cells here need to be run to mask the sources surrounding the object of interest appropriately. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask_objID_selected = selection_mask.index\n",
    "\n",
    "# We intend to expand the ellipses of the sources selected by a user specified factor value 'fac'\n",
    "# Feel free to adjust this size to your preference\n",
    "fac= 30 \n",
    "\n",
    "# Create copy of ellipse mask to allow for re-setting above \n",
    "ellipsemask_wIDs_copy = copy(ellipsemask_wIDs)\n",
    "\n",
    "# Create bigger ellipse around each selected source above\n",
    "\n",
    "for i in mask_objID_selected:\n",
    "    ellipse_dict = mk_ell(image_I2, \n",
    "                          seg_catalog['x'][i], \n",
    "                          seg_catalog['y'][i], \n",
    "                          seg_catalog['theta'][i] * 180. / np.pi - 90., \n",
    "                          1.-q_obj)\n",
    "\n",
    "    ellipsemask_wIDs_copy[(ellipse_dict['ellRad'] < fac*seg_catalog['a'][i])] = i\n",
    "    \n",
    "    # Update mask_objects \n",
    "    mask_objects[(ellipsemask_wIDs_copy == np.float(i))] = np.nan"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot the new, updated mask to see if we have managed to mask the source at (250, 700)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8,8))\n",
    "\n",
    "rgb_alpha_masked = copy(rgb_alpha)\n",
    "\n",
    "master_mask = mask_objects + mask # This should now be our final mask \n",
    "\n",
    "for k in range(0, 4):\n",
    "    rgb_alpha_masked[:,:,k] = rgb_alpha_masked[:,:,k] + master_mask \n",
    "    \n",
    "im = ax.imshow(rgb_alpha_masked, origin='bottom')\n",
    "\n",
    "ax.set_xlabel('X [pixels]', fontsize=20)\n",
    "ax.set_ylabel('Y [pixels]', fontsize=20)\n",
    "plt.title('RGBA + UNMASKED + MASKED', fontsize=25)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once you visually verify that the obvious sources surrounding the galaxy have been masked, we can proceed to the next step!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **Part 4: Sky Background Estimation**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to do accurate surface photometry, we need an estimate of the local background, to correct for over- or under-subtraction of the sky local to the galaxy.  We will do this using a fairly traditional method, by determining the mean values in a sample of equal-area boxes containing unmasked pixels.  We distribute these boxes randomly in position about the cutout, rejecting boxes if they contain more than 50% masked pixels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.stats import sigma_clipped_stats\n",
    "meanflux_I2 = np.zeros(1000)\n",
    "sky_pixels = sky_mask + image_I2\n",
    "\n",
    "# Setting up random state for box distribution\n",
    "local_state = np.random.RandomState(2)\n",
    "dim2, dim1 = sky_pixels.shape # Need to keep boxes within the image cutout\n",
    "# Sky box size will be 100 x 100 pixels\n",
    "length = 100\n",
    "height = 100\n",
    "for k in range(0, 1000):    \n",
    "    xc, yc = local_state.randint(low=101, high=dim1-101, size=2)\n",
    "    # Ensuring the boxes land fully within the image cutout\n",
    "    x0 = np.round(xc-length//2)\n",
    "    x1 = np.round(xc+length//2)\n",
    "    y0 = np.round(yc-height//2)\n",
    "    y1 = np.round(yc+height//2) \n",
    "    \n",
    "    if x0 < 0: x0 = 0\n",
    "    if x1 > dim1: x1 = dim1-1\n",
    "    if y0 < 0: y0 = 0\n",
    "    if y1 > dim2: y1 = dim2-1\n",
    "        \n",
    "    pixels = sky_pixels[y0:y1,x0:x1]    \n",
    "    if len(np.isfinite(pixels)[np.isfinite(pixels)]) < 0.5*len(pixels):\n",
    "        meanflux_I2[k] = np.nan\n",
    "        continue    \n",
    "    meanflux_I2[k] = np.nanmean(np.ravel(pixels))\n",
    "\n",
    "stats = sigma_clipped_stats(meanflux_I2[~np.isnan(meanflux_I2)], sigma=3.0)\n",
    "mean_sky_I2 = stats[1]\n",
    "sigma_I2_sky = stats[2]\n",
    "print('Mean background value: {:.3f}'.format(mean_sky_I2))\n",
    "print('Standard Deviation: {:.3f}'.format(sigma_I2_sky ))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From this, we can see that the background value beyond the galaxy is close to 0 in these units.  From our quick check above, we can see what the standard deviation is in surface brightness, which will give us a rough estimate of the depth of this image (in the I-band) down to the noise limit.  We'll use our magnitude zeropoint and pixel scale that we grabbed much earlier."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're going to grab the pixel scale in arcseconds for the surface photometry."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pxScale = wcs.getPixelScale().asArcseconds()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu_lim_I2 = -2.5*np.log10(sigma_I2_sky) + magZp + 2.5*np.log10(pxScale**2)\n",
    "print('Limiting surface brightness in I-band is: {:.2f} mag/arcsec^2'.format(mu_lim_I2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **Part 5: Surface brightness profiles**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have appropriate masks, and now that we know the background sky value, we can create surface brightness profiles of the galaxy. \n",
    "\n",
    "We need the measured elliptical parameters for the galaxy (to characterise its global 2D shape) to compute the average brightness within elliptical annuli increasing in distance from the center of the galaxy. We visualise how this is computed below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image_source_and_sky = image_I2 + master_mask\n",
    "\n",
    "def to_mags(phot_dict, magZp=magZp):\n",
    "    scale_term = 2.5*np.log10(pxScale**2)\n",
    "    intens = phot_dict['mean_I'] # Using mean for now; median likely similar\n",
    "    mu = -2.5*np.log10(intens) + magZp + scale_term\n",
    "    rad = phot_dict['rad']*pxScale # Radius returned in arcseconds\n",
    "    \n",
    "    mu[~np.isfinite(mu)] = 999\n",
    "    \n",
    "    return mu, rad\n",
    "    \n",
    "photometry = {}\n",
    "phot_dict = surfbprof(image_source_and_sky - mean_sky_I2, # Using only the model+sky subtraction\n",
    "                      image_I2*0.0, \n",
    "                      cen, \n",
    "                      PA, \n",
    "                      1-q,\n",
    "                      2, image_I2.shape[0]//2-1, iflag=0)\n",
    "\n",
    "photometry['HSC-I2'] = phot_dict\n",
    "\n",
    "photometry['HSC-I2']['mu'], photometry['HSC-I2']['rad_as'] = to_mags(photometry['HSC-I2'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "photo_I2 = photometry['HSC-I2']\n",
    "rad = photo_I2['rad']\n",
    "sb = photo_I2['mu']\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(21,7), nrows=1, ncols=3)\n",
    "\n",
    "ax[0].imshow(rgb_alpha, norm=colors.LogNorm(), origin='bottom')\n",
    "ax[0].set_xlim((200, 750))\n",
    "ax[0].set_ylim((220, 700))\n",
    "ax[0].set_xlabel('X [pixels]', fontsize=20)\n",
    "ax[0].set_ylabel('Y [pixels]', fontsize=20)\n",
    "ax[0].set_title('RBGA IMAGE', fontsize=25)\n",
    "    \n",
    "ax[1].imshow(rgb_alpha_masked, norm=colors.LogNorm(), origin='bottom')\n",
    "ax[1].set_xlim((200, 750))\n",
    "ax[1].set_ylim((220, 700))\n",
    "ax[1].set_xlabel('X [pixels]', fontsize=20)\n",
    "ax[1].set_ylabel('Y [pixels]', fontsize=20)\n",
    "ax[1].set_title('MASTER MASK', fontsize=25)\n",
    "\n",
    "# Now plot a handful of bins to demonstrate isophotes\n",
    "bins = [75, 150, 220]\n",
    "cols = ['cyan', 'orchid', 'deeppink']\n",
    "lstyle = ['dotted', 'dashed', 'solid']\n",
    "\n",
    "for k in range(0, len(bins)):\n",
    "    # Plot ellipse at each radius in red\n",
    "    ellipse = Ellipse(xy=cen, width=2.* bins[k], height=2.* bins[k]*q, angle= PA,\\\n",
    "                      edgecolor=cols[k], fc='None', lw=3, ls=lstyle[k])\n",
    "\n",
    "    ax[1].add_patch(ellipse)\n",
    "    \n",
    "# Finally, plotting resulting surface brightness profile \n",
    "mu_bins = np.linspace(16, 30, num=10)\n",
    "\n",
    "ax[2].plot(rad, sb, 'k.')\n",
    "ax[2].set_ylim(16, 30)\n",
    "ax[2].set_xlim(0, 300)\n",
    "\n",
    "for i in range(0, len(bins)):\n",
    "    ax[2].plot(np.ones(len(mu_bins))*bins[i], mu_bins, cols[i], ls=lstyle[i], lw=3)\n",
    "\n",
    "ax[2].invert_yaxis()\n",
    "ax[2].set_ylabel(r'$\\mu_{I2}$ [mag/arcsec$^2$]', fontsize=20)\n",
    "ax[2].set_xlabel('radius [pixels]', fontsize=20)\n",
    "ax[2].set_title('SB PROFILE', fontsize=25)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The left panel shows the galaxy in a log scale, with the color image embedded at high surface brightness.  The central panel shows the same, but with a handful of ellipses marking the isophote shapes we are using to calculate the surface brightness profile. The over-plotted ellipses mark distinctive structures of the spiral galaxy: the bar (cyan), outer spiral arms (light purple) and disk (dark pink). Finally, the actual surface brightness profile is shown in the rightmost panel, with vertical lines showing the locations of the ellipses in the center panel.\n",
    "\n",
    "We can see it has a fairly typical barred galaxy profile: a bulge in the center, a slight depression around 30 pixels where the bar's unstable Lagrange points have cleared out the stellar populations, humps demarcating the inner and outer rings, and a generally exponential decline.  We note that around 230 pixels in radius, however, there is a drastic decline in surface brightness.  This occurs well above the noise limit we calculated earlier (closer to 27 than to the ~28.5 limit from the background scatter), implying either a sharp disk truncation, or a systematic oversubtraction during the image reduction.\n",
    "\n",
    "This was done using the PDR2-equivalent pipeline, suggesting that this pipeline's sky subtraction still manages to under-estimate flux for large angular-sized objects and low surface brightness objects to some extent."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "LSST",
   "language": "python",
   "name": "lsst"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}