#!/usr/bin/env python3

import os
import sys
import astropy.io.fits as pyfits
import numpy

if __name__ == "__main__":

    filtername = sys.argv[1]

    _x = sys.argv[2].split("-")
    x1 = int(_x[0])
    x2 = int(_x[1])

    _y = sys.argv[3].split("-")
    y1 = int(_x[0])
    y2 = int(_x[1])

    out_fn = sys.argv[4]

    # rebin = 1
    try:
        rebin = int(sys.argv[5])
        print("Applying on-the-fly binning x %d" % (rebin))
    except:
        rebin = 1

    nx = x2-x1+1
    ny = y2-y1+1

    pad = 100
    # generate large output image
    img_out = numpy.zeros(((ny*4000+2*pad)//rebin, (nx*4000+2*pad)//rebin))
    print("Final output image size: %d x %d pixels" % (img_out.shape[1], img_out.shape[0]))

    prim_hdr = None
    img_hdr = None
    for x in range(nx): #x1,x2+1):
        for y in range(ny): #y1, y2+1):
            fn_dict = dict(filter=filtername, x=x+x1, y=y+y1)
            in_fn = "%(filter)s/0/%(x)d,%(y)d/calexp-%(filter)s-0-%(x)d,%(y)d.fits" % fn_dict
            if (not os.path.isfile(in_fn)):
                continue
            in_hdu = pyfits.open(in_fn)
            print("Reading image from %s" % (in_fn))

            img = in_hdu[1].data

            if (rebin == 1):
                img = in_hdu[1].data
            else:
                _img = in_hdu[1].data
                out_size_x = _img.shape[1] // rebin
                out_size_y = _img.shape[0] // rebin
                operation = numpy.nansum
                rb1 = numpy.array(numpy.reshape(_img, (out_size_x, rebin, out_size_y, rebin)), dtype=numpy.float32)
                rb2 = operation(rb1, axis=-1)
                img = operation(rb2, axis=1)

            _x1 = (x * 4000) // rebin
            _y1 = (y * 4000) // rebin
            _x2 = _x1 + img.shape[1]
            _y2 = _y1 + img.shape[0]
            img_out[_y1:_y2, _x1:_x2] = img

            if (prim_hdr is None):
                prim_hdr = in_hdu[0].header
            if (img_hdr is None):
                img_hdr = in_hdu[1].header

    # adjust WCS for binning
    for h in ['CRPIX1', 'CRPIX2']:
        img_hdr[h] /= rebin
    for h in ['CD1_1', 'CD1_2', 'CD2_1', 'CD2_2']:
        img_hdr[h] *= rebin

    out_hdulist = pyfits.HDUList([
        pyfits.PrimaryHDU(header=prim_hdr),
        pyfits.ImageHDU(header=img_hdr, data=img_out)
    ])
    out_hdulist.writeto(out_fn, overwrite=True)