-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathalign_image_ciao.py
executable file
·176 lines (144 loc) · 5.31 KB
/
align_image_ciao.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
#! /usr/astro/ciao-4.6/bin/python
#Author: Kari A. Frank
#Date: November 7, 2014
#Purpose: Align a deconvolved, smoothed, ACIS image with a reference image.
#Usage: align_image.py targetimg refimg targetfpars reffpars [--niterations niteraions] [--burnin burnin] [--clobber clobber]
#
#Input:
#
# targetimg = name of fits file containing the image to be shifted
# refimg = name of fits file containing the reference image
#
# targetfpars = name of fpars file for target image
# reffpars = name of fpars file for reference image
#
# niterations = number of iterations to use for the MCMC (default = 1000)
# burnin = number of iterations to disregard (default = 500)
#
# clobber = specifies whether files should be overwritten
# if they already exist (same as in ciao tools,
# default = 'no')
#
#Output:
# - a translated copy of targetimg that should be aligned with refimg
#
#Usage Notes:
#---------------------------------------
# Import Modules
#---------------------------------------
import argparse,os,sys
import pyfits as fits
import ciao_contrib.runtool as crt #all functions should be prefixed with crt
from home_grown import ls_to_list
import numpy as np
import fitting
#import home_grown as hg #import my own custom functions
#---------------------------------------
# Parse arguments and set defaults
#---------------------------------------
parser = argparse.ArgumentParser(description='Align the target image with the reference image.')
pwd = os.getcwd()
parser.add_argument('targetimg',help='Target fits image.')
parser.add_argument('refimg',help='Reference fits image.')
parser.add_argument('targetfpars',help='Target fpars file.')
parser.add_argument('reffpars',help='Reference fpars file.')
parser.add_argument('--niterations',help='Number of MCMC iterations',default=1000)
parser.add_argument('--burnin',help='Number of MCMC iterations to disregard',default=500)
parser.add_argument('--clobber',help='Overwrite existing files.',default='no')
args = parser.parse_args()
args.burnin = int(args.burnin)
args.niterations = int(args.niterations)
#----Set file paths----
#---------------------------------------
# Read in fpars files and set priors
#---------------------------------------
# read files and headers
target_fpars_file = fits.open(args.targetfpars)
ref_fpars_file = fits.open(args.reffpars)
target_head = fits.getheader(args.targetfpars)
ref_head = fits.getheader(args.reffpars)
# get tables
target_table = target_fpars_file[1].data
ref_table = ref_fpars_file[1].data
# get values
target_x0 = target_table.field('X0')[0]
target_x0err = target_table.field('X0ERR')[0]
target_y0 = target_table.field('Y0')[0]
target_y0err = target_table.field('Y0ERR')[0]
ref_x0 = ref_table.field('X0')[0]
ref_x0err = ref_table.field('X0ERR')[0]
ref_y0 = ref_table.field('Y0')[0]
ref_y0err = ref_table.field('Y0ERR')[0]
# calculate deltaX and deltaY and associated sigmas
dx = target_x0 - ref_x0
#dy = target_y0 - ref_y0
dy = ref_y0 - target_y0
sigma_dx = 5.0*(target_x0err**2.0 + ref_x0err**2.0)**0.5
sigma_dy = 5.0*(target_y0err**2.0 + ref_y0err**2.0)**0.5
print dx,dy,sigma_dx,sigma_dy
#---------------------------------------
# Prep Image Files
#---------------------------------------
# get obsids
target_obsid = crt.dmkeypar(args.targetimg,"OBS_ID")
ref_obsid = crt.dmkeypar(args.refimg,"OBS_ID")
# Create normalized images
# Target image
# get normalization constant
targetimgfile = fits.open(args.targetimg)
targetimg = targetimgfile[0].data
targetimgfile.close()
targetimgtotal = np.sum(targetimg)
# divide image
crt.dmimgcalc.punlearn()
crt.dmimgcalc.infile = args.targetimg
crt.dmimgcalc.infile2 = 'none'
crt.dmimgcalc.outfile = args.targetimg+'_normalized'
crt.dmimgcalc.op = 'add'
crt.dmimgcalc.weight = 1.0/targetimgtotal
crt.dmimgcalc.clobber = args.clobber
crt.dmimgcalc()
# Reference image
# get normalization constant
refimgfile = fits.open(args.refimg)
refimg = refimgfile[0].data
refimgfile.close()
refimgtotal = np.sum(refimg)
# divide image
crt.dmimgcalc.punlearn()
crt.dmimgcalc.infile = args.refimg
crt.dmimgcalc.infile2 = 'none'
crt.dmimgcalc.outfile = args.refimg+'_normalized'
crt.dmimgcalc.op = 'add'
crt.dmimgcalc.weight = 1.0/refimgtotal
crt.dmimgcalc.clobber = args.clobber
crt.dmimgcalc()
#---------------------------------------
# Run MCMC
#---------------------------------------
mcmc_pars = fitting.mcmc_align_img(args.targetimg+'_normalized',args.refimg+'_normalized',[dx,dy],[sigma_dx,sigma_dy],niterations=args.niterations)
final_dx = np.median(mcmc_pars[args.burnin:,1])
final_dy = np.median(mcmc_pars[args.burnin:,2])
print 'final dx,dy = ',final_dx,final_dy
sigmadx = np.std(mcmc_pars[args.burnin:,1])
sigmady = np.std(mcmc_pars[args.burnin:,2])
print 'sigmadx,sigmady = ',sigmadx,sigmady
#---------------------------------------
# Translate Image with Best dx, dy
#---------------------------------------
# copy image
crt.dmcopy.punlearn()
crt.dmcopy.infile = args.targetimg
crt.dmcopy.outfile = args.targetimg+'_trans'
crt.dmcopy.clobber = args.clobber
crt.dmcopy()
# shift image
crt.wcs_update.punlearn()
crt.wcs_update.infile = args.targetimg+'_trans'
crt.wcs_update.wcsfile = args.targetimg+'_trans'
crt.wcs_update.deltax = str(final_dx)
crt.wcs_update.deltay = str(final_dy)
crt.wcs_update.rotang = 0.0
crt.wcs_update.scalefac = 1.0
crt.wcs_update.clobber = 'yes'
crt.wcs_update()