Show Tract.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assemble and display a tract"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import collections\n",
    "import numpy as np\n",
    "\n",
    "import lsst.daf.persistence as dafPersist\n",
    "import lsst.afw.image as afwImage\n",
    "import lsst.afw.detection as afwDetect\n",
    "import lsst.afw.geom as afwGeom\n",
    "import lsst.afw.cameraGeom.utils as cameraGeomUtils\n",
    "import lsst.afw.display.ds9 as ds9\n",
    "import lsst.afw.display.utils as displayUtils\n",
    "import lsst.afw.display.rgb as afwRgb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import os\n",
    "os.environ[\"SUPRIME_DATA_DIR\"] = \"/work/yasuda/rerun/yasuda/deep2\"\n",
    "refDates = {\"HSC-G\": \"2014-10-01\",\n",
    "            \"HSC-R\": \"2014-09-22\",\n",
    "            \"HSC-I\": \"2014-09-22\",\n",
    "            \"HSC-Z\": \"2014-10-01\",\n",
    "            \"HSC-Y\": \"2014-09-18\"}\n",
    "imgDates = {'HSC-G': '2015-07-14',\n",
    "            'HSC-R': '2015-07-15',\n",
    "            'HSC-I': '2015-07-20',\n",
    "            'HSC-Z': '2015-07-22',\n",
    "            'HSC-Y': '2015-07-23'}\n",
    "\n",
    "butlers = {}\n",
    "for f in imgDates:\n",
    "    rerun = imgDates[f]\n",
    "    dataPath = os.path.join(os.environ[\"SUPRIME_DATA_DIR\"], \"rerun\", rerun)\n",
    "\n",
    "    butlers[f] = dafPersist.Butler(dataPath)\n",
    "skymap = butlers.values()[0].get(\"deepCoadd_skyMap\")\n",
    "\n",
    "# butler type for difference images\n",
    "if False:\n",
    "    deepCoadd_diff = \"deepCoadd_tempExp_diff\"\n",
    "    deepCoadd_direct = \"deepCoadd_tempExp\"  \n",
    "else:\n",
    "    deepCoadd_diff = \"deepCoadd_diff\"\n",
    "    deepCoadd_direct = \"deepCoadd_calexp\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def getButler(dataId):\n",
    "    return butlers[dataId[\"filter\"]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def assembleTiles(images):\n",
    "    \"\"\"Assemble a list of tiles according to their XY0 values\"\"\"\n",
    "    bigBBox = afwGeom.BoxI()\n",
    "\n",
    "    for im in images:\n",
    "        bigBBox.include(im.getBBox(afwImage.PARENT))\n",
    "\n",
    "    bigIm = afwImage.MaskedImageF(bigBBox)\n",
    "    for im in images:\n",
    "        if True:\n",
    "            sub = bigIm.Factory(bigIm, im.getBBox(afwImage.PARENT), afwImage.PARENT)\n",
    "            sub <<= im.getMaskedImage()\n",
    "            del sub\n",
    "        else:\n",
    "            bigIm[im.getBBox(afwImage.PARENT)] = im\n",
    "\n",
    "    return bigIm\n",
    "\n",
    "#-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-\n",
    "\n",
    "def assemblePatches(patches, patchInfo):\n",
    "    \"\"\"Assemble a list of tiles according to their XY0 values\"\"\"\n",
    "    bigBBox = afwGeom.BoxI()\n",
    "    \n",
    "    for p in patches:\n",
    "        pi = patchInfo[p]\n",
    "        bigBBox.include(pi.getInnerBBox())\n",
    "\n",
    "    bigIm = afwImage.MaskedImageF(bigBBox)\n",
    "    for p in patches:\n",
    "        im = patches[p].getMaskedImage()\n",
    "        pi = patchInfo[p]\n",
    "\n",
    "        if True:\n",
    "            sub = bigIm.Factory(bigIm, pi.getInnerBBox(), afwImage.PARENT)\n",
    "            sub <<= im.Factory(im, pi.getInnerBBox(), afwImage.PARENT)\n",
    "            del sub\n",
    "        else:\n",
    "            bigIm[pi.getInnerBBox()] = im[pi.getInnerBBox()]\n",
    "\n",
    "    wcs = patches[p].getWcs().clone()\n",
    "    wcs.shiftReferencePixel(afwGeom.ExtentD(bigIm.getXY0() - patches[p].getXY0()))\n",
    "\n",
    "    return afwImage.makeExposure(bigIm, wcs)\n",
    "\n",
    "#-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-\n",
    "\n",
    "def getTract(butler, dataId, pps=None, butlerType=\"deepCoadd_calexp\"):\n",
    "    skymap = butler.get(\"deepCoadd_skyMap\")\n",
    "    patchInfo = dict(('%d,%d' % _.getIndex(), _) for _ in skymap[dataId[\"tract\"]])\n",
    "\n",
    "    if pps is None:\n",
    "        pps = list(patchInfo)\n",
    "        \n",
    "    filterName = dataId[\"filter\"]\n",
    "    if len(filterName) == 1:\n",
    "        filterName = \"HSC-%s\" % filterName.upper()\n",
    "        \n",
    "    patches = {}\n",
    "    for pp in pps:\n",
    "        patches[pp] = butler.get(butlerType, dataId, filter=filterName, patch=pp)\n",
    "\n",
    "    return assemblePatches(patches, patchInfo)\n",
    "\n",
    "def getPatchInfo(butler, dataId, pps=None, butlerType=\"deepCoadd_calexp\"):\n",
    "    skymap = butler.get(\"deepCoadd_skyMap\")\n",
    "    patchInfo = dict(('%d,%d' % _.getIndex(), _) for _ in skymap[dataId[\"tract\"]])\n",
    "    \n",
    "    if pps is None:\n",
    "        pps = list(patchInfo)\n",
    "\n",
    "    filterName = dataId[\"filter\"]\n",
    "    if len(filterName) == 1:\n",
    "        filterName = \"HSC-%s\" % filterName.upper()\n",
    "        \n",
    "    info = collections.OrderedDict()\n",
    "    for pp in pps:\n",
    "        patch = butler.get(butlerType, dataId, filter=filterName, patch=pp)\n",
    "        info[pp] = patch.getInfo()\n",
    "\n",
    "    return info"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dataId = dict(tract=9463, filter='HSC-G', patch='3,6')\n",
    "if False:\n",
    "    p = getButler(dataId).get(deepCoadd_direct, dataId, immediate=True)\n",
    "    ds9.mtv(p, frame=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3,6\n",
      "\t34490  0.121\n",
      "\t34494  0.154\n",
      "\t34498  0.122\n",
      "\t34506  0.123\n",
      "\t34510  0.154\n",
      "\t34514  0.117\n",
      "\t34518  0.119\n",
      "\t34520  0.090\n"
     ]
    }
   ],
   "source": [
    "if True:\n",
    "    info = getPatchInfo(getButler(dataId), dataId=dataId, pps=['3,6'])\n",
    "\n",
    "for pp in info:\n",
    "    print pp\n",
    "    inputs = info[pp].getCoaddInputs()\n",
    "    w = 0\n",
    "    for v in inputs.visits:\n",
    "        w += v[\"weight\"]\n",
    "    for v in inputs.visits:\n",
    "        print \"\\t%5d  %.3f\" % (v.getId(), v[\"weight\"]/w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "if False:\n",
    "    d = getButler(dataId).get(deepCoadd_diff, dataId, visit=24164, immediate=True)\n",
    "    ds9.mtv(d, frame=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "pps = [\"%d,%d\" % (i , j) for i in range(3, 6) for j in range(3, 6)]\n",
    "#pps = [\"%d,%d\" % (i , j) for i in range(5, 6) for j in range(5, 6)]\n",
    "pps = ['2,6', '3,6']\n",
    "visits = [-1]   # we're subtracting coadds from each other\n",
    "#visits = [24164, 24196, 24228]\n",
    "\n",
    "fullPatchIm = getTract(getButler(dataId), dataId=dataId, pps=pps,\n",
    "                       butlerType=\"deepCoadd_calexp\")\n",
    "patchIm = fullPatchIm\n",
    "\n",
    "if False:\n",
    "    ds9.mtv(patchIm, frame=1, title=\"%d %s\" % (dataId[\"tract\"], str(pps)[1:-1].replace(\"'\", \"\")))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dataId.update(filter='HSC-R')\n",
    "diffIms = collections.OrderedDict()\n",
    "for frame, v in enumerate(visits, 2):\n",
    "    if v in diffIms:\n",
    "        continue\n",
    "\n",
    "    diffIms[v] = getTract(getButler(dataId), dataId=dict(visit=v, **dataId), pps=pps,\n",
    "                          butlerType=deepCoadd_diff)\n",
    "    if False:\n",
    "        ds9.mtv(diffIms[v], frame=frame, title=v if v > 0 else dataId[\"filter\"])\n",
    "\n",
    "visitIms = collections.OrderedDict()\n",
    "for frame, v in enumerate(diffIms.keys(), 2 + len(visits)):\n",
    "    if v in visitIms:\n",
    "        continue\n",
    "\n",
    "    visitIms[v] = getTract(getButler(dataId), dataId=dict(visit=v, **dataId), pps=pps,\n",
    "                           butlerType=deepCoadd_direct)\n",
    "    if False:\n",
    "        ds9.mtv(visitIms[v], frame=frame, title=v if v > 0 else dataId[\"filter\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "diffImsC = collections.OrderedDict()\n",
    "visitImsC = collections.OrderedDict()\n",
    "\n",
    "filters = ['HSC-G', 'HSC-R', 'HSC-I']\n",
    "for frame, f in enumerate(filters, 2):\n",
    "    dataId['filter'] = f\n",
    "\n",
    "    if f not in diffImsC:\n",
    "        diffImsC[f] = getTract(getButler(dataId), dataId=dict(**dataId), pps=pps,\n",
    "                              butlerType=deepCoadd_diff)\n",
    "    if True:\n",
    "        ds9.mtv(diffImsC[f], frame=frame, title=dataId[\"filter\"])\n",
    "\n",
    "    if f not in visitImsC:\n",
    "        visitImsC[f] = getTract(getButler(dataId), dataId=dict(**dataId), pps=pps,\n",
    "                               butlerType=deepCoadd_direct)\n",
    "    if True:\n",
    "        ds9.mtv(visitImsC[f], frame=frame + len(filters), title=dataId[\"filter\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now some code to mask the artefacts in a difference image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def maskArtefacts(patchIm, diffIms, threshold=2, nGrow=15, brightThreshold=np.inf,\n",
    "                  clippedNpixMin=1000, includeChipGaps=True):\n",
    "    \"\"\"\n",
    "    threshold : nSigma (but the data's correlated, so this isn't directly interpretable)\n",
    "    nGrow     : pixels; saturated footprints only\n",
    "    \"\"\"\n",
    "    fs = afwDetect.FootprintSet(patchIm.getMaskedImage(), \n",
    "                                afwDetect.Threshold(threshold, afwDetect.Threshold.PIXEL_STDEV))\n",
    "\n",
    "    x0, y0 = patchIm.getXY0()\n",
    "    mask = patchIm.getMaskedImage().getMask()\n",
    "    image = patchIm.getMaskedImage().getImage()\n",
    "\n",
    "    CLIPPED = 1 << mask.getMaskPlane(\"CLIPPED\")\n",
    "    SAT = 1 << mask.getMaskPlane(\"SAT\")\n",
    "    INTRP = 1 << mask.getMaskPlane(\"INTRP\")\n",
    "\n",
    "    maskedPixels = (mask.getArray() & SAT)\n",
    "    height, width = maskedPixels.shape\n",
    "\n",
    "    for foot in fs.getFootprints():\n",
    "        peak = foot.getPeaks()[0]\n",
    "        x, y, = peak.getIx(), peak.getIy()\n",
    "        setBits = foot.overlapsMask(mask)\n",
    "        if image.get(x - x0, y - y0) > brightThreshold or (setBits & SAT) or \\\n",
    "           ((setBits & CLIPPED) and foot.getNpix() >= clippedNpixMin):\n",
    "            if (setBits & SAT) and nGrow:\n",
    "                foot = afwDetect.growFootprint(foot, nGrow)\n",
    "            for s in foot.getSpans():\n",
    "                sy = s.getY() - y0\n",
    "                sx0, sx1 = s.getX0() - x0, s.getX1() - x0\n",
    "\n",
    "                try:\n",
    "                    maskedPixels[sy, sx0:sx1+1] = 1\n",
    "                except IndexError:\n",
    "                    if sy < 0 or sy >= height:\n",
    "                        continue\n",
    "                    if sx0 < 0:\n",
    "                        sx0 = 0\n",
    "                    if sx1 >= width:\n",
    "                        sx1 = width - 1\n",
    "                    maskedPixels[sy, sx0:sx1+1] = 1\n",
    "\n",
    "    if includeChipGaps:             # include union of chip gaps in maskedPixels\n",
    "        for v in diffIms:\n",
    "            maskedPixels[np.isnan(diffIms[v][bbox].getMaskedImage().getImage().getArray())] = 1\n",
    "\n",
    "    cleaned = collections.OrderedDict()\n",
    "    for v in diffIms:\n",
    "        cleaned[v] = diffIms[v][bbox].clone()\n",
    "        imArr = cleaned[v].getMaskedImage().getImage().getArray()\n",
    "        mskArr = cleaned[v].getMaskedImage().getMask().getArray()\n",
    "        mskArr &= CLIPPED\n",
    "\n",
    "        bad = np.logical_or(np.isnan(imArr), maskedPixels > 0)        \n",
    "        good = np.logical_and(np.isfinite(imArr), maskedPixels == 0)\n",
    "\n",
    "        q25, q50, q75 = np.percentile(imArr[good], [25, 50, 75])\n",
    "        imArr[bad] = np.random.normal(q50, 0.741*(q75 - q25), size=np.sum(bad))\n",
    "        mskArr[bad] |= INTRP\n",
    "\n",
    "    return cleaned"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "if True:\n",
    "    bbox = fullPatchIm.getBBox()\n",
    "else:\n",
    "    bbox = afwGeom.BoxI(afwGeom.PointI(9000, 9000), afwGeom.PointI(12000 - 1, 12000 - 1))\n",
    "\n",
    "patchIm = fullPatchIm[bbox]\n",
    "\n",
    "if False:\n",
    "    ds9.mtv(patchIm, frame=1, title=\"%d %s\" % (tract, str(pps)[1:-1].replace(\"'\", \"\")))\n",
    "\n",
    "cleaned = maskArtefacts(patchIm, diffIms if False else diffImsC)\n",
    "\n",
    "if False:\n",
    "    frame0 = 2\n",
    "    for frame, v in enumerate(diffIms, frame0):\n",
    "        title = v if len(diffIms) > 12000 else dataId[\"filter\"]\n",
    "\n",
    "        if True:\n",
    "            ds9.mtv(cleaned[v], frame=frame, title=title)\n",
    "        if True:\n",
    "            ds9.mtv(visitImsC[v][bbox], frame=frame + len(diffIms), title=title)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 121,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "if len(cleaned) == 1:\n",
    "    R, G, B = [im.getMaskedImage().getImage().clone() for im in 3*cleaned.values()]\n",
    "    G *= 0\n",
    "    B *= 0\n",
    "else:\n",
    "    R, G, B = reversed([im.getMaskedImage().getImage().clone() for im in cleaned.values()])\n",
    "\n",
    "    # make the blue a bit brighter by mixing in some green\n",
    "    if isinstance(cleaned.keys()[0], int):\n",
    "        G.getArray()[:] += 0.5*B.getArray()\n",
    "\n",
    "# Set the pixel masks\n",
    "mskVal = 0.025\n",
    "for im in [G, B]:\n",
    "    im.getArray()[bad] += mskVal\n",
    "    \n",
    "clippedPixels = (patchIm.getMaskedImage().getMask().getArray() & CLIPPED) == CLIPPED\n",
    "for im in [R, B]:\n",
    "    im.getArray()[clippedPixels] += mskVal\n",
    "\n",
    "#\n",
    "# Add in the ghost of the direct image\n",
    "#\n",
    "addGhost = False\n",
    "if addGhost:\n",
    "    for im in [R, G, B]:\n",
    "        im.getArray()[:] += 0.2*patchIm.getMaskedImage().getImage().getArray()\n",
    "    \n",
    "rgb = afwRgb.makeRGB(R, G, B, \n",
    "                     minimum=-0.1, range=1, Q=4\n",
    "                     #minimum=-.2, range=.4, Q=20\n",
    "                    )\n",
    "afwRgb.writeRGB(\"tract-%d-patches-diff%s.png\" % (dataId[\"tract\"], \"-ghost\" if addGhost else \"\"), rgb)\n",
    "\n",
    "if False:\n",
    "    afwRgb.displayRGB(rgb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rgb = afwRgb.makeRGB(*reversed([im.getMaskedImage().getImage() for im in visitImsC.values()]), \n",
    "                     minimum=-0.1, range=1, Q=4)\n",
    "afwRgb.writeRGB(\"tract-%d-patches-%s.png\" % (dataId[\"tract\"], \"gri\"), rgb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "gray = afwRgb.makeRGB(*3*[patchIm[bbox].getMaskedImage()], minimum=0, range=4, Q=8)\n",
    "afwRgb.writeRGB(\"tract-%d-patches-%s.png\" % (dataId[\"tract\"], \"gray\"), gray)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "if False:\n",
    "    ds9.mtv(butler.get(\"deepCoadd_tempExp\", tract=0, patch='3,3', filter='HSC-G', visit=24228))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[195.5876425659588, 4.032620171976911]\n"
     ]
    }
   ],
   "source": [
    "import lsst.afw.cameraGeom as afwCamGeom\n",
    "boreSight = calexp.getDetector().getCenter()   # pixels in FPA coordinates\n",
    "pointing = calexp.getWcs().pixelToSky(boreSight.getPixels(1)) # where the boresight's pointing\n",
    "print [_.asDegrees() for _ in pointing]"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}