Tips for image data

Create one fits data to entire observed region

Because of the large FOV of HSC, coadd data is output per patch. However, in some cases, one large fits data may be used. The following python script will be useful to generate entire large.fits.

# Example of python script execution
# Specify all coadd images you want to combine.
# hscPipe 3.x
python make-large-image.py $home/hsc/rerun/dith_16h_test/deepCoadd/HSC-I/0/*.fits

# hscPipe 4.x
python make-large-image.py $home/hsc/rerun/dith_16h_test/deepCoadd/HSC-I/0/*/calexp*.fits

The content of the script is shown below;

#!/usr/bin/env python
# encoding: utf-8

# Library about image in AFW (Application FrameWork)
import lsst.afw.image as afwImage
# Library about geometry in AFW
import lsst.afw.geom  as afwGeom

import sys

def main():

        if len(sys.argv) <= 1:
                print "Select more than 1 patch"
                return 1

        largeImage = makeLargeImageFromPatches(sys.argv[1:])

        print "Saving data (Wait a while...)"
        largeImage.writeFits("large.fits")


def makeLargeImageFromPatches(filenames):

        # Get metric information of image
        bboxes = []
        for filename in filenames:
                print "Get bounding box:", filename, "..."
                bboxes.append(
                        # Get bounding box of this patch in parent(=tract).
                        afwImage.ExposureF(filename).getBBox(afwImage.PARENT)
                )

        # Derive a BBox that includes all BBox
        #     Return value of getMin/Max methods is inclusive
        #     (Region = [min, max]) cf. getBegin/End (Use below)
        x0 = min(bbox.getMinX() for bbox in bboxes)
        y0 = min(bbox.getMinY() for bbox in bboxes)
        x1 = max(bbox.getMaxX() for bbox in bboxes)
        y1 = max(bbox.getMaxY() for bbox in bboxes)

        width  = x1 - x0 + 1
        height = y1 - y0 + 1

        # Prepare large image
        print "Prepare for image: ({}, {})".format(width, height)
        exposure = afwImage.ExposureF(width, height)

        # Set position of image in the tract
        exposure.setXY0(afwGeom.Point2I(x0, y0))

        # Extract am image part from Exposure
        largeImage = exposure.getMaskedImage()

        # Paste patches to largeImage
        for filename in filenames:
                print "Get patch:", filename, "..."

        # Load patches again
                patch = afwImage.ExposureF(filename)
                bbox = patch.getBBox(afwImage.PARENT)

                # Position of this patch in the largeImage
                #     Return value of getBegin/End methods is one past the end
                #     (Region = [begin,end) ) cf. getMin/Max (Used above)
                xp0 = bbox.getBeginX() - x0
                yp0 = bbox.getBeginY() - y0
                xp1 = bbox.getEndX  () - x0
                yp1 = bbox.getEndY  () - y0

                # Pase patches (Please refer to numpy documets for the syntax of [])
                largeImage[xp0:xp1, yp0:yp1] = patch.getMaskedImage()[:]

                # Assume that all patches have the same information (e.g. WCS),
                # Set one of them to large image.
                if not exposure.hasWcs() and patch.hasWcs():
                        exposure.setCalib   (patch.getCalib   ())
                        exposure.setFilter  (patch.getFilter  ())
                        exposure.setMetadata(patch.getMetadata())
                        exposure.setWcs     (patch.getWcs     ())

        return exposure


if __name__ == "__main__":
        main()

Figure 1 shows the image of large.fits made with this script via ds9.

Fig 1. Image of large.fits