color postage

# python color-postage.py -o a.png 2-cutout-HSC-I-9589-pdr1_wide.fits 3-cutout-HSC-R-9589-pdr1_wide.fits 4-cutout-HSC-G-9589-pdr1_wide.fits

import pyfits
import argparse
import numpy
import PIL.Image


def main():
    parser = argparse.ArgumentParser()
    parser.add_argument('--out', '-o', required=True)
    parser.add_argument('fits', nargs=3)
    args = parser.parse_args()

    layers = []
    for f in args.fits:
        with pyfits.open(f) as hdul:
            x = scale(hdul[1].data, hdul[0].header['FLUXMAG0'])
            layers.append(x)

    layers = numpy.array(layers)
    layers[layers < 0] = 0
    layers[layers > 1] = 1
    layers = layers.transpose((1, 2, 0))[::-1, :, :]

    layers = numpy.array(255 * layers, dtype=numpy.uint8)
    img = PIL.Image.fromarray(layers)
    img.save(args.out)


def scale(x, fluxMag0):
    mag0 = 19
    scale = 10 ** (0.4 * mag0) / fluxMag0
    x *= scale

    u_min = -0.05
    u_max = 2. / 3.
    u_a = numpy.exp(10.)

    x = numpy.arcsinh(u_a*x) / numpy.arcsinh(u_a)
    x = (x - u_min) / (u_max - u_min)

    return x


if __name__ == '__main__':
    main()