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
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
|
# WCSINFO --
procedure wcsinfo (image)
string image { prompt="Image name" }
bool verbose = no { prompt="Verbose?" }
string pos = "" { prompt="POS string" }
real size = 0.0 { prompt="SIZE value" }
real ra = 0.0 { prompt="RA value" }
real dec = 0.0 { prompt="DEC value" }
real llx = 0.0 { prompt="LLX wcs corner" }
real lly = 0.0 { prompt="LLY wcs corner" }
real urx = 0.0 { prompt="URX wcs corner" }
real ury = 0.0 { prompt="URY wcs corner" }
real midx = 0.0 { prompt="X wcs midpoint" }
real midy = 0.0 { prompt="Y wcs midpoint" }
real scale = 0.0 { prompt="Scale (\"/pix)" }
real rot = 0.0 { prompt="Rotation (deg)" }
real width = 0.0 { prompt="Plt width (arcmin)" }
real height = 0.0 { prompt="Plt height (arcmin)" }
string axmap = "" { prompt="Axis mapping" }
real epoch = 0.0 { prompt="Epoch" }
real equinox = 0.0 { prompt="Equinox" }
string ctype1 = "" { prompt="CTYPE1" }
string ctype2 = "" { prompt="CTYPE2" }
real crval1 = 0.0 { prompt="CRVAL1" }
real crval2 = 0.0 { prompt="CRVAL2" }
real crpix1 = 0.0 { prompt="CRPIX1" }
real crpix2 = 0.0 { prompt="CRPIX2" }
real cd11 = 0.0 { prompt="CD1_1" }
real cd12 = 0.0 { prompt="CD1_2" }
real cd21 = 0.0 { prompt="CD2_1" }
real cd22 = 0.0 { prompt="CD2_2" }
bool ispix = no { prompt="Pixel only coords?" }
begin
string img, lpos
real cdelt1, cdelt2 # old wcs/pix increment
real crota1, crota2 # old rotation keywords
real sys_epoch, sys_equinox # epoch info
real xrot, yrot # plate scale and rotation
real x1, x2, y1, y2 # corner pos
real lllx, llly, lurx, lury, lmx, lmy
int nx, ny, cx, cy
string imtmp, img_orig, cfile, ofile
bool verb
# Check for proper packages and reset environment.
reset imtype = "fits"
#flpr 0
# Get the task parameters.
img = image
verb = verbose
# Initialize.
cd11 = 1.0 ; cd12 = 0.0
cd21 = 0.0 ; cd22 = 1.0
crpix1 = 0.0 ; crpix2 = 0.0
crval1 = 0.0 ; crval2 = 0.0
cdelt1 = 0.0 ; cdelt2 = 0.0
ctype1 = "" ; ctype2 = ""
crota1 = 0.0 ; crota2 = 90.0 # assume normal x/y axes and
scale = 1.0 ; rot = 0.0 # default to 1"/pix scale
xrot = 0.0 ; yrot = 0.0 ; rot = 0.0
sys_epoch = 0.0
sys_equinox = 0.0
# Work on a copy of the image, not the original.
imtmp = mktemp ("/tmp/w")
imcopy (img, imtmp, verbose-)
img_orig = img
img = imtmp
# Get the size and midpoint of the image (in pixels).
hselect (img, "i_naxis1,i_naxis2", yes) | scan (nx, ny)
cx = nx / 2
cy = ny / 2
# First, see whether we have FITS standard keywords or a plate
# solution header (e.g. from DSS). If the latter, fix it.
hselect (img, "cd1_1,cd2_2", yes) | scan (x, y)
if (nscan() < 2) {
hselect (img, "amdx1,amdy1", yes) | scan (x, y)
if (nscan() >= 1) {
if (verb)
print ("Converting platesol header to std WCS... ")
makewcs (img, verbose-)
}
}
# Now convert the image explicitly to FK5. This will also give us a
# standard CD matrix instead of the CROTA/CDELT keywords
imcctran (img, "FK5", nx=20, ny=20, longpole-, verbose-, update+)
hselect (img, "crval1,crval2", yes) | scan (crval1, crval2)
if (nscan() != 2) {
if (verb)
print ("No CRVAL keywords")
}
hselect (img, "crpix1,crpix2", yes) | scan (crpix1, crpix2)
if (nscan() != 2) {
if (verb)
print ("No CRPIX keywords")
}
hselect (img, "ctype1,ctype2", yes) | scan (ctype1, ctype2)
if (nscan() == 2) {
s1 = strlwr (ctype1)
if (strstr ("dec", s1) == 1 || strstr ("lat", s1) == 1) {
axmap = "2 1"
} else if (strstr ("ra", s1) == 1 || strstr ("lon", s1) == 1) {
axmap = "1 2"
} else {
axmap = "1 2"
ispix = yes
}
;
if (strstr ("-sip", strlwr (s1)) > 0) {
ctype1 = substr (ctype1, 1, (strlen(ctype1)-4))
ctype2 = substr (ctype2, 1, (strlen(ctype2)-4))
hedit (img, "ctype1", ctype1, verify-,update+,show-, >& "dev$null")
hedit (img, "ctype2", ctype2, verify-,update+,show-, >& "dev$null")
hedit (img, "WAT0_001", "", verify-,del+,update+,show-,>&"dev$null")
hedit (img, "WAT1_001", "", verify-,del+,update+,show-,>&"dev$null")
hedit (img, "WAT2_001", "", verify-,del+,update+,show-,>&"dev$null")
}
;
# Check for an invalid WAT keyword suggesting physical coords. If we
# have them and CTYPE keywords, delete the WAT values.
hselect (img, "WAT0_001", "yes") | scan (s1)
if (s1 == "system=physical") {
hedit (img, "WAT0_001", del+, verify-, show-, update+)
hedit (img, "WAT1_001", del+, verify-, show-, update+)
hedit (img, "WAT2_001", del+, verify-, show-, update+)
hedit (img, "LTM1_1", del+, verify-, show-, update+)
hedit (img, "LTM2_2", del+, verify-, show-, update+)
hedit (img, "LTV1", del+, verify-, show-, update+)
hedit (img, "LTV2", del+, verify-, show-, update+)
}
;
} else if (verb) {
print ("No CTYPE keywords")
}
if (ispix) {
# See if we have an RA/DEC keyword and put it at the middle of the image.
hselect (img, "ra,dec", yes) | scan (x, y)
if (nscan() == 2) {
crval1 = x * 15.
if (x > 24)
crval1 = x
crval2 = y
crpix1 = cx
crpix2 = cy
scale = 1.
rot = 0.
}
} else {
# If we have CD matrix keywords use 'em, otherwise try to fall back
# to the older CDELT keywords.
hselect (img, "cd1_1,cd2_2", yes) | scan (x, y)
if (nscan() >= 1) {
# Get the CD matrix from the image. We initialize the values above
# so if they don't exist in the header we won't change the value here
# (but we need to read them one at a time).
hselect (img, "cd1_1", yes) | scan (cd11)
hselect (img, "cd1_2", yes) | scan (cd12)
hselect (img, "cd2_1", yes) | scan (cd21)
hselect (img, "cd2_2", yes) | scan (cd22)
# Compute the plate scale (arcsec/pixel) for the image.
scale = 3600. * sqrt ((cd11**2+cd21**2+cd12**2+cd22**2)/2.)
xrot = abs (datan ( cd21, cd11))
yrot = abs (datan (-cd12, cd22))
rot = (xrot + yrot) / 2.0
} else {
hselect (img, "cdelt1,cdelt2", yes) | scan (cdelt1,cdelt2)
if (nscan() == 2)
scale = 3600. * sqrt ((cdelt1**2 + cdelt2**2) / 2.)
hselect (img, "crota1", yes) | scan (crota1)
if (nscan() == 1)
xrot = crota1
hselect (img, "crota2", yes) | scan (crota2)
if (nscan() == 1)
yrot = crota2
rot = yrot
}
}
# Get the epoch and equinox
x = -1.0 ; y = -1.0
hselect (img, "epoch", yes) | scan (x)
if (nscan() == 1 && x > 0.0)
sys_epoch = x
else
sys_epoch = 2000.0
hselect (img, "equinox", yes) | scan (x)
if (nscan() == 1 && y > 0.0)
sys_equinox = y
else
sys_equinox = 2000.0
epoch = sys_epoch
equinox = sys_equinox
if (ispix) {
lllx = (crval1 - (cx/3600.))/15.; # assumes hours of RA
llly = (crval2 - (cy/3600.));
lurx = (crval1 + (cx/3600.))/15.; # assumes hours of RA
lury = (crval2 + (cy/3600.));
lmx = crval1 / 15.
lmy = crval2
llx = lllx * 15.0 ; lly = llly
urx = lurx * 15.0 ; ury = lury
midx = lmx * 15.0 ; midy = lmy
} else {
cfile = mktemp ("/tmp/cf")
ofile = mktemp ("/tmp/of")
print ("1.0 1.0", > cfile)
print (nx // " " // ny, >> cfile)
print (cx // " " // cy, >> cfile)
skyctran (cfile, ofile, img, img//" world", >& "dev$null")
tail (ofile, nl=1) | head ("STDIN", nl=1) | scan (x, y, lmx, lmy )
tail (ofile, nl=2) | head ("STDIN", nl=1) | scan (x, y, lurx, lury)
tail (ofile, nl=3) | head ("STDIN", nl=1) | scan (x, y, lllx, llly)
# Update the task parameters.
if (axmap == "1 2") {
llx = lllx * 15.0 ; lly = llly
urx = lurx * 15.0 ; ury = lury
midx = lmx * 15.0 ; midy = lmy
ra = midx ; dec = midy
} else {
lly = lllx * 15.0 ; llx = llly
ury = lurx * 15.0 ; urx = lury
midy = lmx * 15.0 ; midx = lmy
ra = midy ; dec = midx
}
printf ("%g,%g\n", ra, dec) | scan (lpos)
pos = lpos
size = max ((scale * nx), (scale * ny)) / 3600. # in degrees
delete (cfile, verify-, >& "dev$null")
delete (ofile, verify-, >& "dev$null")
}
width = real (scale * nx) / 60. # plate width in arcmin
height = real (scale * ny) / 60. # plate height in arcmin
if (verb) {
print ("CRPIX1 = " // crpix1)
print ("CRPIX2 = " // crpix2)
print ("CTYPE1 = " // ctype1)
print ("CTYPE2 = " // ctype2)
printf ("CRVAL1 = %g\t/ X wcs = %.3H\n", crval1, crval1)
printf ("CRVAL2 = %g\t/ Y wcs = %.3h\n", crval2, crval2)
print ("CD1_1 = " // cd11)
print ("CD1_2 = " // cd12)
print ("CD2_1 = " // cd21)
print ("CD2_2 = " // cd22)
print ("CDELT1 = " // cdelt1) ; print ("CDELT2 = " // cdelt2)
print ("CROTA1 = " // crota1) ; print ("CROTA2 = " // crota2)
print ("EPOCH = " // sys_epoch)
print ("EQUINOX = " // sys_equinox)
print ("PLTSCALE= " // scale // " / arc/pix")
print ("PLT_ROT = " // rot // " / degrees")
printf ("PLT_WID = %-14.9g\t/ Field width (arcmin)\n", width);
printf ("PLT_HGT = %-14.9g\t/ Field height (arcmin)\n", height);
printf ("PLT_LLX = %-14.9g\t/ LL X wcs = %.3h\n", lllx*15., lllx)
printf ("PLT_URX = %-14.9g\t/ UR X wcs = %.3h\n", lurx*15., lurx)
printf ("PLT_MX = %-14.9g\t/ Mid X wcs = %.3h\n", lmx*15., lmx)
printf ("PLT_LLY = %-14.9g\t/ LL Y wcs = %.3h\n", llly, llly)
printf ("PLT_URY = %-14.9g\t/ UR Y wcs = %.3h\n", lury, lury)
printf ("PLT_MY = %-14.9g\t/ Mid Y wcs = %.3h\n", lmy, lmy)
printf ("AXMAP = %s\t\t\t/ Axis mapping\n", axmap)
printf ("ISPIX = %b\t\t\t/ Pixel mapping?\n", ispix)
}
# Clean up.
if (imaccess (imtmp) == yes)
delete (imtmp//".fits", verify-, >& "dev$null")
end
|
