Use miriad to regrid/reproject a set of images to make a nice stack.

multi.py

#! /usr/bin/env python 

from __future__ import print_function, division

import numpy as np
from subprocess import Popen
import os
import shutil


def regrid_wrapper(mirfile, outfile, crval, crpix, cdelt, naxis, proj="ZEA"):
    with open("regrid.log", "a+") as log:
        Popen("regrid axes=1,2 in={} out={} desc={},{},{},{},{},{},{},{} project={}".format(
              mirfile, outfile, crval[0], crpix[0], cdelt[0], naxis[0],
              crval[1], crpix[1], cdelt[1], naxis[1], proj), stderr=log, stdout=log, 
              shell=True).wait()

def reproject(ref, indir="./", coords=(None, None), dx=None, dy=None, \
	          remove_old=True, nopixel=False):
	"""Reproject/regrid all images to `ref`.

	All images will be in the same projection as `ref` and 
	will have the same pixel sizes. 
	"""


	dirs = ["mir/", "int/", "regrid/", "final/"]
	if not indir.endswith("/"): indir += "/"
	for d in dirs:
		if not os.path.exists(indir+d): 
			os.mkdir(indir+d)
		elif os.path.exists(indir+d) and remove_old: 
			shutil.rmtree(indir+d)
			os.mkdir(indir+d)





	# Convert `ref` file to MIRIAD format:
	Popen("fits in={0} out={1}mir/ref.tmp op=xyin".format(indir+ref, indir), \
		  shell=True).wait()




	for spec in os.listdir(indir):
		if spec.endswith(".fits"):

			print("Starting {0}...".format(spec))

			Popen("fits in={0} out={1} op=xyin".format(\
				  indir+spec, indir+"mir/"+spec), shell=True).wait()






			Popen("regrid in={0} out={1} axes=1,2 tin={2}/mir/ref.tmp".format(\
				  indir+"mir/"+spec, indir+"regrid/"+spec, \
				  indir), shell=True).wait()

			Popen("fits in={0} out={1} op=xyout".format(\
				  indir+"regrid/"+spec, indir+"int/"+spec), \
			      shell=True).wait()

	# shutil.copy2(indir+ref, indir+"int/"+ref)

	if None not in coords:

		from stamper.coord import hms_dms_dd

		try: coords = hms_dms_dd(coords[0], coords[1])
		except Exception: pass

		if dy is None: dy = dx
		if dx is None: raise ValueError("Please specify cut-out size.")

		for spec in os.listdir(indir+"int/"):

			Popen("mSubimage {0} {1} {2} {3} {4} {5}".format(\
				  indir+"int/"+spec, indir+"final/"+spec, coords[0], \
				  coords[1], dx, dy), shell=True).wait()





def main():
	"""
	"""

	from argparse import ArgumentParser

	ps = ArgumentParser()
	ps.add_argument("ref", type=str)
	ps.add_argument("-i", "--indir", default="./", type=str)
	ps.add_argument("--nopixel", action="store_true")

	args = ps.parse_args()

	reproject(args.ref, args.indir)


if __name__ == "__main__":
	main()