1#!/usr/bin/env python 2"""Module to difference two FITS images.""" 3 4import sys 5import os 6import numpy as np 7import pyfits 8 9 10def save_fits_image(filename, data, header, img1=None, img2=None): 11 """Save a FITS image.""" 12 # Reshape to add the frequency axis 13 data = np.reshape(data, (1, 1, data.shape[0], data.shape[1])) 14 new_hdr = pyfits.header.Header() 15 for i, item in enumerate(header.items()): 16 if item[0] != 'HISTORY': 17 new_hdr.append(item) 18 if img1 and img2: 19 new_hdr.append(('HISTORY', '' * 60)) 20 new_hdr.append(('HISTORY', '-' * 60)) 21 new_hdr.append(('HISTORY', 'Diff created from image1 - image2:')) 22 new_hdr.append(('HISTORY', '- image1 : %s' % img1)) 23 new_hdr.append(('HISTORY', '- image2 : %s' % img2)) 24 new_hdr.append(('HISTORY', '-' * 60)) 25 new_hdr.append(('HISTORY', ' - Max : % .3e' % np.max(data))) 26 new_hdr.append(('HISTORY', ' - Min : % .3e' % np.min(data))) 27 new_hdr.append(('HISTORY', ' - Mean : % .3e' % np.mean(data))) 28 new_hdr.append(('HISTORY', ' - STD : % .3e' % np.std(data))) 29 new_hdr.append(('HISTORY', ' - RMS : % .3e' % 30 np.sqrt(np.mean(data**2)))) 31 new_hdr.append(('HISTORY', '-' * 60)) 32 new_hdr.append(('HISTORY', '' * 60)) 33 34 if (os.path.exists(filename)): 35 print '+ WARNING, output FITS file already exsits, overwriting.' 36 os.remove(filename) 37 pyfits.writeto(filename, data, new_hdr) 38 39 40def load_fits_image(filename): 41 """Load a FITS image.""" 42 43 hdulist = pyfits.open(filename) 44 data = hdulist[0].data 45 hdr = hdulist[0].header 46 return np.squeeze(data), hdr 47 48 49def fits_diff(outname, file1, file2): 50 if not os.path.exists(file1): 51 print 'ERROR: missing input file', file1 52 return 53 if not os.path.exists(file2): 54 print 'ERROR: missing input file', file2 55 return 56 if os.path.exists(outname): 57 return 58 59 f1, h1 = load_fits_image(file1) 60 f2, h2 = load_fits_image(file2) 61 diff = f1 - f2 62 print '-' * 80 63 print '+ Image size : %i x %i' % (f1.shape[0], f1.shape[1]) 64 print '+ File 1 : %s' % file1 65 print '+ File 2 : %s' % file2 66 print '+ Diff : file1 - file2' 67 print '+ Output name : %s' % (outname) 68 print '+ Diff stats:' 69 print ' - Max : % .3e' % np.max(diff) 70 print ' - Min : % .3e' % np.min(diff) 71 print ' - Mean : % .3e' % np.mean(diff) 72 print ' - STD : % .3e' % np.std(diff) 73 print ' - RMS : % .3e' % np.sqrt(np.mean(diff**2)) 74 print '-' * 80 75 save_fits_image(outname, diff, h1, file1, file2) 76 77 78if __name__ == '__main__': 79 80 if len(sys.argv) - 1 != 3: 81 print 'Usage: fits_diff.py <diff name> <file1> <file2>' 82 print '' 83 print 'Performs: diff = file1 - file2' 84 sys.exit(1) 85 86 outname = sys.argv[1] 87 file1 = sys.argv[2] 88 file2 = sys.argv[3] 89 90 f1, h1 = load_fits_image(file1) 91 f2, h2 = load_fits_image(file2) 92 93 diff = f1 - f2 94 95 print '-' * 80 96 print '+ Image size : %i x %i' % (f1.shape[0], f1.shape[1]) 97 print '+ File 1 : %s' % file1 98 print '+ File 2 : %s' % file2 99 print '+ Diff : file1 - file2' 100 print '+ Output name : %s' % (outname) 101 print '+ Diff stats:' 102 print ' - Max : % .3e' % np.max(diff) 103 print ' - Min : % .3e' % np.min(diff) 104 print ' - Mean : % .3e' % np.mean(diff) 105 print ' - STD : % .3e' % np.std(diff) 106 print ' - RMS : % .3e' % np.sqrt(np.mean(diff**2)) 107 print '-' * 80 108 109 save_fits_image(outname, diff, h1, file1, file2) 110