1
2
3 """
4 Creates Source Extractor catalogues from deep/mosaic images in CASU-like
5 format.
6
7 @author: N.J. Cross
8 @org: WFAU, IfA, University of Edinburgh
9
10 @todo: Calculate seeing,avStellEll, go deeper with classification
11 @todo: import self from CU13 - use isTrialRun, redo ...
12 @todo: use SExtractor 286 - no longer need to break out extensions.
13
14 @requires: PyFITS, Numpy, Scientific, SExtractor
15 """
16
17 import math
18 import numpy
19 import os
20 import pyfits
21 from Scientific.Functions import LeastSquares as LSQ
22 import shutil
23
24 import wsatools.DbConnect.DbConstants as dbc
25 import wsatools.FitsUtils as fits
26 from wsatools.Logger import Logger
27 import wsatools.Statistics as stats
28 import wsatools.Utilities as utils
29 import wsatools.ExternalProcess as extp
30
31
33 """ Provides an interface between CU13 and functions to process an
34 image with SExtractor
35 """
36 comment = "Running SExtractor"
37
38
40 self.archive = archive
41 self.sysc = archive.sysc
42 self.bestImage = bestImageObj
43 self.dualImage = dualImageObj
44 self.isDualImage = dualImageObj is not None
45
46
47
49 """
50 """
51 maxSize = 60
52 donotChangeList = ['NAXIS', 'PCOUNT', 'GCOUNT', 'BITPIX', 'SIMPLE', 'EXTEND',
53 'NAXIS1', 'NAXIS2', 'XTENSION']
54
55 utils.ParsedFile.stripEOLComments = True
56 inparDict = dict(line.split(None, 1)
57 for line in utils.ParsedFile(self.inparFileName))
58
59 hdulist = fits.open(self.bestImage.catName(), 'update')
60 imghdu = fits.open(self.bestImage.fileName(), 'update')
61
62 pHdrCards = imghdu[0].header.ascard
63 if self.bestImage.isMef():
64 sHdrCards = [hdu.header.ascard
65 for det, hdu in enumerate(imghdu) if det > 0]
66 else:
67 sHdrCards = [imghdu[0].header.ascard]
68
69 for det, hdu in enumerate(hdulist):
70 if det > 0:
71 nxout = imghdu[det].header['NAXIS1']
72 nyout = imghdu[det].header['NAXIS2']
73
74 try:
75 sNoise = hdu.data.field('Sky_rms')
76 sLevel = hdu.data.field('Sky_level')
77 except:
78 sNoise = hdu.data.field('THRESHOLD')
79 sLevel = hdu.data.field('BACKGROUND')
80 try:
81 gaussSigm = hdu.data.field('Gaussian_sigma')
82 except:
83 gaussSigm = hdu.data.field('IMAGE_FWHM') / 2.355
84
85
86
87
88 fwhm = 2.355 * gaussSigm
89 ell = hdu.data.field('Ellipticity')
90 skynoise = stats.clippedMADHistogram(sNoise)
91 skylevel = stats.clippedMADHistogram(sLevel)
92 starMask = self.selectStars(hdu)
93 seeing, avStellEll = self.calcSeeingEll(fwhm, ell, starMask)
94 hdu.header.add_history(utils.makeDateTime().strftime("%Y%m%d")
95 + ' Adding in primary header from Image.')
96
97
98 for card in pHdrCards:
99 if card.key not in donotChangeList:
100 keyword = ('HIERARCH ' + card.key
101 if (len(card.key) > 8 and not
102 card.key.startswith('HIERARCH'))
103 else card.key)
104 hdu.header.update(keyword, card.value, card.comment)
105
106
107 sHCards = sHdrCards[det - 1]
108
109
110 hdu.header.add_history(utils.makeDateTime().strftime("%Y%m%d")
111 + ' Adding in extension header from Image.')
112 for card in sHCards:
113 if card.key not in donotChangeList:
114 keyword = ('HIERARCH ' + card.key
115 if (len(card.key) > 8 and not
116 card.key.startswith('HIERARCH'))
117 else card.key)
118 hdu.header.update(keyword, card.value, card.comment)
119 try:
120 hdu.header.update('CRPIX1',
121 imghdu[det].header['CRPIX1'],
122 'Dither offset X')
123 hdu.header.update('CRPIX2',
124 imghdu[det].header['CRPIX2'],
125 'Dither offset Y')
126 hdu.header.update(
127 'CRVAL1', imghdu[det].header['CRVAL1'],
128 '[deg] Right ascension at the reference pixel')
129 hdu.header.update(
130 'CRVAL2', imghdu[det].header['CRVAL2'],
131 '[deg] Declination at the reference pixel')
132 except:
133 raise ValueError(
134 "Missing CRPIX values in %s, det %s" %
135 (self.bestImage.imageName(), (det)))
136 if self.bestImage.sysc.isVSA():
137
138 hdu.header.update('TCRVL3', hdu.header['CRVAL1'],
139 '[deg] Right ascension at the reference pixel')
140 hdu.header.update('TCRVL5', hdu.header['CRVAL2'],
141 '[deg] Declination at the reference pixel')
142 hdu.header.update('TCTYP3', hdu.header['CTYPE1'],
143 'Algorithm type for axis 1')
144 hdu.header.update('TCTYP5', hdu.header['CTYPE2'],
145 'Algorithm type for axis 2')
146 hdu.header.update('TCRPX3', hdu.header['CRPIX1'],
147 'Dither offset X')
148 hdu.header.update('TCRPX5', hdu.header['CRPIX2'],
149 'Dither offset Y')
150 hdu.header.update('TC3_3', hdu.header['CD1_1'],
151 'Transformation matrix element')
152 hdu.header.update('TC3_5', hdu.header['CD1_2'],
153 'Transformation matrix element')
154 hdu.header.update('TC5_3', hdu.header['CD2_1'],
155 'Transformation matrix element')
156 hdu.header.update('TC5_5', hdu.header['CD2_2'],
157 'Transformation matrix element')
158
159
160
161
162
163
164
165 hdu.header.update('NXOUT', nxout, 'X-axis size of image')
166 hdu.header.update('NYOUT', nyout, 'Y-axis size of image')
167
168
169 thrCount = float(inparDict['DETECT_THRESH']) * skynoise[0]
170 hdu.header.update('THRESHOL', round(thrCount, 3),
171 'Isophotal analysis threshold (counts)')
172 hdu.header.update('SKYLEVEL', round(skylevel[0], 3),
173 'Median sky brightness (counts/pixel)')
174 hdu.header.update('SKYNOISE', round(skynoise[0], 3),
175 'Pixel noise at sky level (counts)')
176 hdu.header.update('SEEING', round(seeing, 3),
177 'Average FWHM (pixels)')
178 hdu.header.update('ELLIPTIC', round(avStellEll, 3),
179 'Ellipticity estimate from SEXWRAP')
180
181 if self.bestImage.isMef():
182 imghdu[det].header.update('SKYLEVEL', round(skylevel[0], 3),
183 'Median sky brightness (counts/pixel)')
184 imghdu[det].header.update('SKYNOISE', round(skynoise[0], 3),
185 'Pixel noise at sky level (counts)')
186 else:
187 imghdu[0].header.update('SKYLEVEL', round(skylevel[0], 3),
188 'Median sky brightness (counts/pixel)')
189 imghdu[0].header.update('SKYNOISE', round(skynoise[0], 3),
190 'Pixel noise at sky level (counts)')
191
192 imghdu[0].header.update('EXTINCT', 0.0, 'Extinction')
193 imghdu[0].header.update(self.bestImage.sysc.amStart, 1.0,
194 'Airmass at start of observation')
195 imghdu[0].header.update(self.bestImage.sysc.amEnd, 1.0,
196 'Airmass at end of observation')
197
198 hdu.header.update('MINPIX', float(inparDict['DETECT_MINAREA']),
199 'Minimum size for images (pixels)')
200 hdu.header.update('RCORE', 1. / float(inparDict['PIXEL_SCALE']),
201 'Core radius for default profile fit (pixels)')
202
203
204
205 hdu.header.update('CLASSIFD', True,
206 'Class flag: -1 stellar, 1 non-stellar, 0 noise')
207 hdu.header.update('APCORPK', 0.,
208 'Stellar aperture correction - peak height')
209 if self.bestImage.sysc.isWSA():
210 hdu.header.update('APCOR1', 0.,
211 'Stellar aperture correction - 1/2x core flux')
212 hdu.header.update('APCOR2', 0.,
213 'Stellar aperture correction - core/sqrt(2) flux')
214 hdu.header.update('APCOR3', 0.,
215 'Stellar aperture correction - core flux')
216 hdu.header.update('APCOR4', 0.,
217 'Stellar aperture correction - 1.41x core flux')
218 hdu.header.update('APCOR5', 0.,
219 'Stellar aperture correction - 2x core flux')
220 hdu.header.update('APCOR6', 0.,
221 'Stellar aperture correction - 2.82x core flux')
222 hdu.header.update('APCOR7', 0.,
223 'Stellar aperture correction - 4x core flux')
224 comment = 'SExtractor Param'
225 for kWrd in inparDict:
226 Val = inparDict[kWrd]
227 if len(kWrd) > 8:
228 kWrd = kWrd[:8]
229 if kWrd == 'FILTER':
230 kWrd = 'CNFILTER'
231 if len(Val) > maxSize - len(comment):
232 if os.sep in Val:
233 Val = os.path.join(os.path.basename(
234 os.path.dirname(Val)), os.path.basename(Val))
235 if len(Val) > maxSize - len(comment):
236 Val = '++' + Val[len(comment) - maxSize - 2:]
237 else:
238 Val = '++' + Val[len(comment) - maxSize - 2:]
239 hdulist[det].header.update(kWrd, Val, comment)
240
241 imghdu.close()
242 hdulist.close()
243
244
245
247 """
248 """
249 compAper = 8
250 apertures = range(1, 8, 1)
251 aperSize = ['1/2x core', 'core/sqrt(2)', 'core', '1.41x core', '2x core',
252 '2.82x core', '4x core']
253 hdulist = fits.open(self.bestImage.catName(), 'update')
254 for ii, hdu in enumerate(hdulist):
255 if ii > 0:
256 tbdata = hdu.data
257 avStellEll = hdu.header['ELLIPTIC']
258 seeing = hdu.header['SEEING']
259 aperMag3 = self.zeropoint + tbdata.field('Aper_flux_3')
260 try:
261 fwhm = tbdata.field('Gaussian_sigma') * 2.355
262 except:
263 fwhm = tbdata.field('IMAGE_FWHM')
264
265 apFluxMsk = numpy.where(numpy.greater(tbdata.field('Aper_flux_3'), 0.),
266 tbdata.field('Aper_flux_3'), 1.)
267 mask = numpy.where(numpy.greater(tbdata.field('Aper_flux_3'), 0.),
268 numpy.logical_and(numpy.logical_and(numpy.logical_and(
269 numpy.logical_and(
270 numpy.logical_and(numpy.greater(aperMag3, (self.sLimit + 2)),
271 numpy.less(numpy.divide(
272 tbdata.field('Aper_flux_3_err'), apFluxMsk), 0.01)),
273 fwhm < (1.2 * seeing)),
274 numpy.less_equal(tbdata.field('Ellipticity'), (1.2 * avStellEll))),
275 numpy.logical_or(numpy.equal(tbdata.field('Classification'), -2),
276 numpy.equal(tbdata.field('Classification'), -1))),
277 numpy.equal(tbdata.field('Error_bit_flag'), 0)), 0)
278 for aper in apertures:
279 headerAttr = ('WAPCOR%d' % aper
280 if self.bestImage.sysc.isWSA() else
281 'APCOR%d' % aper)
282 apFluxMsk = numpy.where(numpy.greater(tbdata.field('Aper_flux_%d' % compAper), 0.),
283 tbdata.field('Aper_flux_%d' % compAper), 1.)
284 fluxRatios = numpy.where(
285 numpy.greater(tbdata.field('Aper_flux_%d' % compAper), 0.),
286 numpy.multiply(mask, numpy.divide(
287 tbdata.field('Aper_flux_%d' % aper),
288 apFluxMsk)), 0.)
289 indices = numpy.nonzero(numpy.multiply(
290 fluxRatios, numpy.greater(fluxRatios, 0)))
291 apCors = numpy.multiply(-2.5, numpy.log10(
292 fluxRatios[indices]))
293 apCor, _psd = stats.clippedMADHistogram(apCors, minSd=0.0001)
294 hdu.header.update(headerAttr, round(apCor, 4),
295 'Stellar aperture correction - %s flux'
296 % aperSize[aper - 1])
297 hdulist.close()
298
299
300
302 """Use segmentation and weight image to calculate average confidence
303 """
304 hdulist = fits.open(self.bestImage.getSegmName())
305 segdata = hdulist[0].data
306 hdulist.close()
307 hdulist = fits.open(self.bestImage.confName())
308 confdata = hdulist[0].data
309 hdulist.close()
310 prod = segdata * confdata
311 avConf = numpy.zeros([max(segdata)])
312 for objID in range(1, max(segdata) + 1, 1):
313 numer = numpy.sum(numpy.where(numpy.equal(segdata, objID), prod, 0.))
314 denom = numpy.sum(numpy.where(numpy.equal(segdata, objID), objID, 0.))
315 if denom != 0.:
316 avConf[objID - 1] = int(numer / denom)
317 return avConf
318
319
320
322 """Calculates average seeing and ellipticity
323 """
324
325 fwhmStar = fwhm[numpy.nonzero(mask)]
326 ellStar = ell[numpy.nonzero(mask)]
327
328 medSeeing, _psdSeeing = stats.clippedMADHistogram(fwhmStar, 0.)
329 medStellEll, _psdStellEll = stats.clippedMADHistogram(ellStar, 0.)
330 return medSeeing, medStellEll
331
332
333
336 """Function to calculate a discrete classifier based on
337 magnitude, magnitude error and stat.
338 -1 stellar,+1 non-stellar,0 noise,-2 borderline stellar,
339 -3 borderline non-stellar,-9 saturated.
340 """
341 nObj = len(mag)
342 disClass = numpy.zeros([nObj])
343
344
345 self.getLimits(mag, magErr, errorLimit)
346
347 disClass = numpy.add(disClass, numpy.multiply(numpy.logical_and(
348 numpy.less(mag, self.mLimit), numpy.greater_equal(
349 stat, brightSGCut)), -1))
350
351 disClass = numpy.add(disClass, numpy.multiply(numpy.logical_and(
352 numpy.less(mag, self.mLimit), numpy.less_equal(stat, brightSGCut)), 1))
353
354 disClass = numpy.add(disClass, numpy.multiply(numpy.logical_and(
355 numpy.greater_equal(mag, self.mLimit), numpy.greater(
356 stat, faintSGCut)), -2))
357
358 disClass = numpy.add(disClass, numpy.multiply(numpy.logical_and(
359 numpy.greater_equal(mag, self.mLimit), numpy.less_equal(
360 stat, faintSGCut)), -3))
361
362 disClass = numpy.multiply(disClass, numpy.logical_not(numpy.logical_and(
363 numpy.less(mag, self.sLimit), numpy.greater(stat, satSGCut))))
364 disClass = numpy.add(disClass, numpy.multiply(numpy.logical_and(
365 numpy.less(mag, self.sLimit), numpy.greater(stat, satSGCut)), -9))
366
367 return numpy.add(disClass, numpy.multiply(numpy.equal(disClass, 0), -2))
368
369
370
371
372
374 """This function reads in data from the ascii file and converts
375 to a fits file.
376 """
377
378 orderedListWSA = ['Sequence_number', 'Isophotal_flux', 'X_coordinate',
379 'X_coordinate_err', 'Y_coordinate', 'Y_coordinate_err',
380 'Gaussian_sigma', 'Ellipticity', 'Position_angle',
381 'Areal_1_profile', 'Areal_2_profile', 'Areal_3_profile',
382 'Areal_4_profile', 'Areal_5_profile', 'Areal_6_profile',
383 'Areal_7_profile', 'Areal_8_profile', 'Peak_height',
384 'Peak_height_err', 'Aper_flux_1', 'Aper_flux_1_err',
385 'Aper_flux_2', 'Aper_flux_2_err', 'Aper_flux_3',
386 'Aper_flux_3_err', 'Aper_flux_4', 'Aper_flux_4_err',
387 'Aper_flux_5', 'Aper_flux_5_err', 'Aper_flux_6',
388 'Aper_flux_6_err', 'Aper_flux_7', 'Aper_flux_7_err',
389 'Aper_flux_8', 'Aper_flux_8_err', 'Aper_flux_9',
390 'Aper_flux_9_err', 'Aper_flux_10', 'Aper_flux_10_err',
391 'Aper_flux_11', 'Aper_flux_11_err', 'Aper_flux_12',
392 'Aper_flux_12_err', 'Aper_flux_13', 'Aper_flux_13_err',
393 'Petr_radius', 'Kron_radius', 'Hall_radius', 'Petr_flux',
394 'Petr_flux_err', 'Kron_flux', 'Kron_flux_err', 'Hall_flux',
395 'Hall_flux_err', 'Error_bit_flag', 'Sky_level', 'Sky_rms',
396 'Parent_or_child', 'RA', 'DEC', 'Classification', 'Statistic',
397 'Blank63', 'Blank64', 'Blank65', 'Blank66', 'Blank67',
398 'Blank68', 'Blank69', 'Blank70', 'Blank71', 'Blank72',
399 'Blank73', 'Blank74', 'Blank75', 'Blank76', 'Blank77',
400 'Blank78', 'Blank79', 'Blank80']
401 finalUnitsWSA = ['Number', 'ADU', 'Pixels', 'Pixels', 'Pixels', 'Pixels',
402 'Pixels', 'Number', 'Degrees', 'Pixels', 'Pixels', 'Pixels',
403 'Pixels', 'Pixels', 'Pixels', 'Pixels', 'Pixels', 'ADU',
404 'ADU', 'ADU', 'ADU', 'ADU', 'ADU', 'ADU', 'ADU', 'ADU',
405 'ADU', 'ADU', 'ADU', 'ADU', 'ADU', 'ADU', 'ADU', 'ADU',
406 'ADU', 'ADU', 'ADU', 'ADU', 'ADU', 'ADU', 'ADU', 'ADU',
407 'ADU', 'ADU', 'ADU', 'Pixels', 'Pixels', 'Pixels', 'ADU',
408 'ADU', 'ADU', 'ADU', 'ADU', 'ADU', 'Number', 'ADU', 'ADU',
409 'Number', 'RADIANS', 'RADIANS', 'Flag', 'N-sigma',
410 'Blank63', 'Blank64', 'Blank65', 'Blank66', 'Blank67',
411 'Blank68', 'Blank69', 'Blank70', 'Blank71', 'Blank72',
412 'Blank73', 'Blank74', 'Blank75', 'Blank76', 'Blank77',
413 'Blank78', 'Blank79', 'Blank80']
414 orderedListVSA = ['Sequence_number', 'Isophotal_flux', 'X_coordinate',
415 'X_coordinate_err', 'Y_coordinate', 'Y_coordinate_err',
416 'Gaussian_sigma', 'Ellipticity', 'Position_angle',
417 'Areal_1_profile', 'Areal_2_profile', 'Areal_3_profile',
418 'Areal_4_profile', 'Areal_5_profile', 'Areal_6_profile',
419 'Areal_7_profile', 'Areal_8_profile', 'Peak_height',
420 'Peak_height_err', 'Aper_flux_1', 'Aper_flux_1_err',
421 'Aper_flux_2', 'Aper_flux_2_err', 'Aper_flux_3',
422 'Aper_flux_3_err', 'Aper_flux_4', 'Aper_flux_4_err',
423 'Aper_flux_5', 'Aper_flux_5_err', 'Aper_flux_6',
424 'Aper_flux_6_err', 'Aper_flux_7', 'Aper_flux_7_err',
425 'Aper_flux_8', 'Aper_flux_8_err', 'Aper_flux_9',
426 'Aper_flux_9_err', 'Aper_flux_10', 'Aper_flux_10_err',
427 'Aper_flux_11', 'Aper_flux_11_err', 'Aper_flux_12',
428 'Aper_flux_12_err', 'Aper_flux_13', 'Aper_flux_13_err',
429 'Petr_radius', 'Kron_radius', 'Half_radius', 'Petr_flux',
430 'Petr_flux_err', 'Kron_flux', 'Kron_flux_err', 'Half_flux',
431 'Half_flux_err', 'Error_bit_flag', 'Sky_level', 'Sky_rms',
432 'Av_conf', 'RA', 'DEC', 'Classification', 'Statistic',
433 'Blank63', 'Blank64', 'Blank65', 'Blank66', 'Blank67',
434 'Blank68', 'Blank69', 'Blank70', 'Blank71', 'Blank72',
435 'Blank73', 'Blank74', 'Blank75', 'Blank76', 'Blank77',
436 'Blank78', 'Blank79', 'Blank80']
437 finalUnitsVSA = ['Number', 'ADU', 'Pixels', 'Pixels', 'Pixels', 'Pixels',
438 'Pixels', 'Number', 'Degrees', 'Pixels', 'Pixels', 'Pixels',
439 'Pixels', 'Pixels', 'Pixels', 'Pixels', 'Pixels', 'ADU',
440 'ADU', 'ADU', 'ADU', 'ADU', 'ADU', 'ADU', 'ADU', 'ADU',
441 'ADU', 'ADU', 'ADU', 'ADU', 'ADU', 'ADU', 'ADU', 'ADU',
442 'ADU', 'ADU', 'ADU', 'ADU', 'ADU', 'ADU', 'ADU', 'ADU',
443 'ADU', 'ADU', 'ADU', 'Pixels', 'Pixels', 'Pixels', 'ADU',
444 'ADU', 'ADU', 'ADU', 'ADU', 'ADU', 'Number', 'ADU', 'ADU',
445 'Number', 'RADIANS', 'RADIANS', 'Flag', 'N-sigma',
446 'Blank63', 'Blank64', 'Blank65', 'Blank66', 'Blank67',
447 'Blank68', 'Blank69', 'Blank70', 'Blank71', 'Blank72',
448 'Blank73', 'Blank74', 'Blank75', 'Blank76', 'Blank77',
449 'Blank78', 'Blank79', 'Blank80']
450
451
452
453 dictCASUCatNames = {
454 'Sequence_number':'NUMBER',
455 'Isophotal_flux':'FLUX_ISO',
456 'X_coordinate':'MathLookUp',
457 'X_coordinate_err':'MathLookUp',
458 'Y_coordinate':'MathLookUp',
459 'Y_coordinate_err':'MathLookUp',
460 'Gaussian_sigma':'MathLookUp',
461 'Ellipticity':'ELLIPTICITY',
462 'Position_angle':'MathLookUp',
463 'Areal_1_profile':'ISO0',
464 'Areal_2_profile':'ISO1',
465 'Areal_3_profile':'ISO2',
466 'Areal_4_profile':'ISO3',
467 'Areal_5_profile':'ISO4',
468 'Areal_6_profile':'ISO5',
469 'Areal_7_profile':'ISO6',
470 'Areal_8_profile':'ISO7',
471 'Peak_height':'FLUX_MAX',
472 'Peak_height_err':'MathLookUp',
473 'Aper_flux_1':'MathLookUp',
474 'Aper_flux_1_err':'MathLookUp',
475 'Aper_flux_2':'MathLookUp',
476 'Aper_flux_2_err':'MathLookUp',
477 'Aper_flux_3':'MathLookUp',
478 'Aper_flux_3_err':'MathLookUp',
479 'Aper_flux_4':'MathLookUp',
480 'Aper_flux_4_err':'MathLookUp',
481 'Aper_flux_5':'MathLookUp',
482 'Aper_flux_5_err':'MathLookUp',
483 'Aper_flux_6':'MathLookUp',
484 'Aper_flux_6_err':'MathLookUp',
485 'Aper_flux_7':'MathLookUp',
486 'Aper_flux_7_err':'MathLookUp',
487 'Aper_flux_8':'MathLookUp',
488 'Aper_flux_8_err':'MathLookUp',
489 'Aper_flux_9':'MathLookUp',
490 'Aper_flux_9_err':'MathLookUp',
491 'Aper_flux_10':'MathLookUp',
492 'Aper_flux_10_err':'MathLookUp',
493 'Aper_flux_11':'MathLookUp',
494 'Aper_flux_11_err':'MathLookUp',
495 'Aper_flux_12':'MathLookUp',
496 'Aper_flux_12_err':'MathLookUp',
497 'Aper_flux_13':'MathLookUp',
498 'Aper_flux_13_err':'MathLookUp',
499 'Petr_radius':'MathLookUp',
500 'Kron_radius':'MathLookUp',
501 'Hall_radius':'realDefault',
502 'Petr_flux':'FLUX_PETRO',
503 'Petr_flux_err':'FLUXERR_PETRO',
504 'Kron_flux':'FLUX_AUTO',
505 'Kron_flux_err':'FLUXERR_AUTO',
506 'Hall_flux':'realDefault',
507 'Hall_flux_err':'realDefault',
508 'Error_bit_flag':'FLAGS',
509 'Sky_level':'BACKGROUND',
510 'Sky_rms':'THRESHOLD',
511 'Parent_or_child':'intDefault',
512 'RA':'MathLookUp',
513 'DEC':'MathLookUp',
514 'Classification':'MathLookUp',
515 'Statistic':'CLASS_STAR',
516 'Blank63':'realDefault',
517 'Blank64':'realDefault',
518 'Blank65':'realDefault',
519 'Blank66':'realDefault',
520 'Blank67':'realDefault',
521 'Blank68':'realDefault',
522 'Blank69':'realDefault',
523 'Blank70':'realDefault',
524 'Blank71':'realDefault',
525 'Blank72':'realDefault',
526 'Blank73':'realDefault',
527 'Blank74':'realDefault',
528 'Blank75':'realDefault',
529 'Blank76':'realDefault',
530 'Blank77':'realDefault',
531 'Blank78':'realDefault',
532 'Blank79':'realDefault',
533 'Blank80':'realDefault',
534 'Half_radius':'MathLookUp',
535 'Half_flux':'realDefault',
536 'Half_flux_err':'realDefault',
537 'Av_conf':'intDefault',
538 }
539
540
541
542 asciiCatName = self.bestImage.getAsciiCatName()
543 orderedList, finalUnits = None, None
544 if self.bestImage.sysc.isWSA():
545 orderedList = orderedListWSA
546 finalUnits = finalUnitsWSA
547 if self.bestImage.sysc.isVSA():
548 orderedList = orderedListVSA
549 finalUnits = finalUnitsVSA
550 ascFile = open(asciiCatName, 'r')
551 lines = ascFile.readlines()
552 ascFile.close()
553 header = []
554 for line in lines:
555 if line[0] == '#':
556
557 colNo, colName, units, comment = splitHeadLine(line)
558 header.append((colNo, colName, units, comment))
559 nDataCols = len(lines[len(header)].split())
560 dArray = map(numpy.array,
561 utils.extractColumns(asciiCatName, range(nDataCols), dataType=float))
562
563 if minX != -999.:
564 dArray = selectData(dArray, header, minX, maxX, minY, maxY)
565 nDataRows = len(dArray[0])
566 indx = 0
567
568 colArray = []
569
570 limits = None
571 for casuAttr in orderedList:
572 sexOut = dictCASUCatNames[casuAttr]
573 unit = finalUnits[indx]
574 indx += 1
575 if sexOut == 'intDefault':
576 dVector = numpy.ones([nDataRows]) * dbc.intDefault()
577 elif sexOut == 'realDefault':
578 dVector = numpy.ones([nDataRows]).astype(numpy.float32) * dbc.realDefault()
579 elif sexOut == 'MathLookUp':
580 if casuAttr == 'X_coordinate':
581 for ii in range(len(header)):
582 if header[ii][1] == 'X_IMAGE':
583 ximg = dArray[header[ii][0] - 1]
584 if minX > 0.:
585 dVector = numpy.subtract(ximg, minX)
586 else:
587 dVector = ximg
588 elif casuAttr == 'Y_coordinate':
589 for ii in range(len(header)):
590 if header[ii][1] == 'Y_IMAGE':
591 yimg = dArray[header[ii][0] - 1]
592 if minY > 0.:
593 dVector = numpy.subtract(yimg, minY)
594 else:
595 dVector = yimg
596 if casuAttr == 'X_coordinate_err':
597 for ii in range(len(header)):
598 if header[ii][1] == 'ERRX2_IMAGE':
599 errx2 = dArray[header[ii][0] - 1]
600 dVector = numpy.add(numpy.sqrt(numpy.multiply(errx2, numpy.greater(errx2, 0.))), numpy.multiply(numpy.less_equal(errx2, 0.), dbc.realDefault()))
601 elif casuAttr == 'Y_coordinate_err':
602 for ii in range(len(header)):
603 if header[ii][1] == 'ERRY2_IMAGE':
604 erry2 = dArray[header[ii][0] - 1]
605 dVector = numpy.add(numpy.sqrt(numpy.multiply(erry2, numpy.greater(erry2, 0.))), numpy.multiply(numpy.less_equal(erry2, 0.), dbc.realDefault()))
606 elif casuAttr == 'Gaussian_sigma':
607 fwhmColL = [header[ii][0] for ii in range(len(header))
608 if header[ii][1] == 'FWHM_IMAGE']
609 if len(fwhmColL) == 1:
610 aColNo = fwhmColL[0] - 1
611 dVector = dArray[aColNo] / 2.355
612 else:
613 aImgColIdx = [ii for ii in range(len(header))
614 if header[ii][1] == 'A_IMAGE'][0]
615 aColNo = header[aImgColIdx][0] - 1
616 bImgColIdx = [ii for ii in range(len(header))
617 if header[ii][1] == 'B_IMAGE'][0]
618 bColNo = header[bImgColIdx][0] - 1
619 dVector = numpy.add(numpy.multiply(numpy.multiply(numpy.sqrt(numpy.add(numpy.power(dArray[aColNo], 2.), numpy.power(dArray[bColNo], 2.))), 0.5), numpy.greater(dArray[bColNo], 0.)), numpy.multiply(numpy.less_equal(dArray[bColNo], 0.), dbc.realDefault()))
620
621 elif casuAttr == 'Peak_height_err':
622 for ii in range(len(header)):
623 if header[ii][1] == 'FLUX_MAX':
624 aColNo = header[ii][0] - 1
625 if header[ii][1] == 'FLUX_APER':
626 bColNo = header[ii][0] - 1
627 if header[ii][1] == 'FLUXERR_APER':
628 cColNo = header[ii][0] - 1
629
630 msk = numpy.where(dArray[bColNo] > 0, dArray[bColNo], 1.)
631 dVector = numpy.where(
632 numpy.greater(dArray[bColNo], 0.),
633 numpy.multiply(dArray[aColNo],
634 numpy.divide(dArray[cColNo], msk)),
635 dbc.realDefault())
636
637 elif casuAttr == 'Position_angle':
638 for ii in range(len(header)):
639 if header[ii][1] == 'THETA_IMAGE':
640 dVector = numpy.add(numpy.multiply(numpy.greater(dArray[header[ii][0] - 1], 0.), dArray[header[ii][0] - 1]), numpy.multiply(numpy.less(dArray[header[ii][0] - 1], 0.), (dArray[header[ii][0] - 1] + 180.)))
641 elif casuAttr == 'Petr_radius':
642 for ii in range(len(header)):
643 if header[ii][1] == 'A_IMAGE':
644 aColNo = header[ii][0] - 1
645 if header[ii][1] == 'PETRO_RADIUS':
646 bColNo = header[ii][0] - 1
647 dVector = numpy.multiply(dArray[aColNo], dArray[bColNo])
648 elif casuAttr == 'Kron_radius':
649 for ii in range(len(header)):
650 if header[ii][1] == 'A_IMAGE':
651 aColNo = header[ii][0] - 1
652 if header[ii][1] == 'KRON_RADIUS':
653 bColNo = header[ii][0] - 1
654 dVector = numpy.multiply(dArray[aColNo], dArray[bColNo])
655 elif casuAttr == 'Half_radius':
656 for ii in range(len(header)):
657 if header[ii][1] == 'A_IMAGE':
658 aColNo = header[ii][0] - 1
659 if header[ii][1] == 'FLUX_RADIUS':
660 bColNo = header[ii][0] - 1
661 dVector = numpy.multiply(dArray[aColNo], dArray[bColNo])
662 elif casuAttr == 'Classification':
663 for ii in range(len(header)):
664 if header[ii][1] == 'CLASS_STAR':
665 aColNo = header[ii][0] - 1
666 if header[ii][1] == 'MAG_APER':
667 bColNo = header[ii][0] + 1
668 if header[ii][1] == 'MAGERR_APER':
669 cColNo = header[ii][0] + 1
670 dVector = self.classification(dArray[aColNo], dArray[bColNo], dArray[cColNo], errorLimit=0.04)
671 elif casuAttr == 'RA':
672 for ii in range(len(header)):
673 if header[ii][1] == 'ALPHA_J2000':
674 aColNo = header[ii][0] - 1
675 dVector = numpy.multiply(dArray[aColNo], 0.017453292)
676 elif casuAttr == 'DEC':
677 for ii in range(len(header)):
678 if header[ii][1] == 'DELTA_J2000':
679 aColNo = header[ii][0] - 1
680 dVector = numpy.multiply(dArray[aColNo], 0.017453292)
681 elif casuAttr == 'Av_conf':
682
683
684 dVector = self.calcAvConf()
685 elif len(casuAttr) > 10 and len(casuAttr) < 13:
686 if casuAttr[:10] == 'Aper_flux_':
687 for ii in range(len(header)):
688 if header[ii][1] == 'FLUX_APER':
689 aColNo = header[ii][0] - 1
690 aperNo = int(casuAttr[10:])
691 dVector = dArray[aColNo + aperNo - 1]
692 elif len(casuAttr) > 13:
693 if casuAttr[:10] == 'Aper_flux_':
694 for ii in range(len(header)):
695 if header[ii][1] == 'FLUXERR_APER':
696 aColNo = header[ii][0] - 1
697 aperNo = int(casuAttr[10:-4])
698 dVector = dArray[aColNo + aperNo - 1]
699 else:
700 for ii in range(len(header)):
701 if header[ii][1] == sexOut:
702 dVector = dArray[header[ii][0] - 1]
703 colM = pyfits.Column(name=casuAttr, format='E', unit=unit, array=dVector)
704 colArray.append(colM)
705 cols = pyfits.ColDefs(colArray)
706 hdu = pyfits.new_table(cols)
707 if limits:
708 hdu.header.update("MGSAT", limits[0], "Saturation limit, mags")
709 hdu.header.update("MGCLSS", limits[1], "Good class limit, mags")
710 return hdu
711
712
740
741
742
744 """Find best fit function.
745 """
746 minMag, maxMag, bSize = 9., int(max(magArray)), 0.5
747 if maxMag > 30: maxMag = 30.
748 nBins = int((maxMag - minMag) / bSize)
749 data = []
750 dataML = []
751 for bb in range(nBins):
752 mnmag = minMag + (bb + 0.5) * bSize
753 mgErrs = []
754
755 mgErrs = magErrArray[numpy.nonzero(
756 numpy.logical_and(numpy.logical_and(numpy.logical_and(
757 numpy.greater(magArray, mnmag - (0.5 * bSize)), numpy.less_equal(
758 magArray, mnmag + (0.5 * bSize))), numpy.greater(magErrArray, 0.)),
759 numpy.less(magErrArray, 99.)))]
760
761 if mgErrs.size > 5:
762 meanErr, sdErr = stats.clippedMADHistogram(mgErrs, minSd=0.0001)
763
764 meanErr = float(meanErr)
765 sdErr = float(sdErr)
766 else:
767 meanErr, sdErr = 0., 99.
768 if sdErr < 99 and sdErr > 0.and meanErr > 0. and meanErr < 0.5:
769 data.append((mnmag, meanErr, sdErr))
770 if sdErr > 0. and meanErr > 0.3 * rmsLimit and meanErr < 3.*rmsLimit:
771 dataML.append((mnmag, meanErr, sdErr))
772 params = [50., -5., 0.1]
773 fit = LSQ.polynomialLeastSquaresFit(params, dataML)
774 self.sLimit = min(mag for mag, meanErr, _sdErr in data if meanErr < 0.01) + 0.5
775
776 dMag = 0.01
777 test = 0
778 minMag = data[0][0]
779 while (test == 0):
780 if minMag > dataML[0][0] and minMag < dataML[len(dataML) - 1][0]:
781 expErr = LSQ._polynomialModel(fit[0], minMag)
782 if abs(expErr - rmsLimit) < 0.01:
783 self.mLimit = minMag
784 return
785 else:
786 minMag += dMag
787 else:
788 minMag += 10.*dMag
789
790
791
793 """
794 Copy the detectionCatalog.inpar and detectionCatalog.param, and substitute
795 all relevant strings.
796
797 """
798 pxScl = fits.getPixelScale(self.bestImage.fileName())
799 aperRatios = [0.5, math.sqrt(0.5), 1., math.sqrt(2.), 2., 2.*math.sqrt(2.), 4.,
800 5., 6., 7., 8., 10., 12.]
801
802
803
804 inparFile = self.bestImage.sexInparFile()
805 if not inparFile:
806 return
807
808 sexInputsDir = self.bestImage.sysc.sExtractorInputsPath()
809 weightImages = [self.bestImage.confName()]
810 if self.dualImage:
811 weightImages = [self.dualImage.confName()] + weightImages
812
813 paramFile = None
814 if self.bestImage.sysc.isVSA():
815 paramFile = os.path.join(sexInputsDir, 'detectionCatalogVSA.param')
816 if self.bestImage.sysc.isWSA():
817 paramFile = os.path.join(sexInputsDir, 'detectionCatalogWSA.param')
818 dCore = 2. / pxScl
819 catName = self.bestImage.getAsciiCatName()
820 apertures = ','.join(str(aper * dCore) for aper in aperRatios)
821 self.zeropoint = self.bestImage.getZeropoint()
822 replacementList = [('@catalogName', catName),
823 ('@catalogType', 'ASCII_HEAD'),
824 ('@paramFile', paramFile),
825 ('@sxInputDir', sexInputsDir),
826 ('@apertures', apertures),
827 ('@pixelScale', str(pxScl)),
828 ('@weightImage', ','.join(weightImages)),
829 ('@zeropoint', '%s' % self.zeropoint)]
830 if self.bestImage.sysc.isVSA():
831 replacementList += [('@checkImageType', 'SEGMENTATION'),
832 ('@checkImageName', self.bestImage.getSegmName())]
833 else:
834 replacementList += [('@checkImageType', 'NONE'),
835 ('@checkImageName', 'NONE')]
836
837 newLines = []
838 for line in file(inparFile):
839 for attr, value in replacementList:
840 line = line.replace(attr, value)
841
842 newLines.append(line)
843
844 self.inparFileName = self.bestImage.getInparName()
845 file(self.inparFileName, "w").writelines(newLines)
846
847
849 """ Run Source Extractor
850 """
851 command = 'sex %s%s -c %s'
852 command %= (self.dualImage.fileName() + ' ' if isDualImageMode else '',
853 self.bestImage.fileName(), self.inparFileName)
854
855 if not os.path.exists(self.inparFileName.replace('inpar', 'cat')):
856 Logger.addMessage("Running SExtractor: %s" % command)
857 extp.run(command, raiseOnError=False, parseStdOut=False)
858
859
860
862 """ Select stars
863 """
864
865 classObj = hdu.data.field('Classification')
866 flags = hdu.data.field('Error_bit_flag')
867 aperMag3 = self.zeropoint + hdu.data.field('Aper_flux_3')
868 apFluxMsk = numpy.where(hdu.data.field('Aper_flux_3') > 0, hdu.data.field('Aper_flux_3'), 1.)
869 errorFrac = numpy.where(numpy.greater(hdu.data.field('Aper_flux_3'), 0.),
870 (hdu.data.field('Aper_flux_3_err') /
871 apFluxMsk), 0.)
872
873 starMask = numpy.logical_and(numpy.logical_and(numpy.logical_and(numpy.logical_or(
874 numpy.equal(classObj, -1), numpy.equal(classObj, -2)),
875 numpy.equal(flags, 0)), numpy.greater(aperMag3, (self.sLimit + 4))),
876 numpy.less_equal(errorFrac, 0.05))
877 return starMask
878
879
880
882 """
883 """
884
885 phdu = pyfits.PrimaryHDU()
886
887 phdu.header.add_comment("FITS (Flexible Image Transport System) format is defined in 'Astronomy")
888 phdu.header.add_comment("and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H")
889 return phdu
890
891
892
893
894 -def selectData(dArray, header, minX, maxX, minY, maxY):
895 """Selects data in correct part of mosaic only
896
897 """
898
899
900 for colNo, colName, _units, _comment in header:
901 if colName == 'X_IMAGE':
902 xColNo = colNo - 1
903 elif colName == 'Y_IMAGE':
904 yColNo = colNo - 1
905
906 tDArray = numpy.transpose(numpy.array(dArray))
907 newArray = []
908 for ii in range(len(dArray[0])):
909 if tDArray[ii][xColNo] >= minX and tDArray[ii][xColNo] < maxX and \
910 tDArray[ii][yColNo] >= minY and tDArray[ii][yColNo] < maxY:
911 newArray.append(tDArray[ii])
912 return numpy.transpose(newArray)
913
914
915
917 """Routine to find which x and y coordinates go to with
918 which image.
919 """
920 xycoordDict = dict()
921 for image in ingestImages:
922 imgHdu = fits.open(image)
923 crpix1 = imgHdu[1].header['CRPIX1']
924 crpix2 = imgHdu[1].header['CRPIX2']
925 naxis1 = imgHdu[1].header['NAXIS1']
926 naxis2 = imgHdu[1].header['NAXIS2']
927 imgHdu.close()
928 minMaxXY = ((crpix1 + 0.5), (crpix1 + naxis1 - 0.5), (crpix2 + 0.5),
929 (crpix2 + naxis2 - 0.5))
930 xycoordDict[image] = minMaxXY
931
932 return xycoordDict
933
934
935
937 """Simple rms model
938 """
939 sky = params[0]
940 zp = params[1]
941 flux = 10 ** (0.4 * (zp - mag))
942 if (flux + sky) >= 0 and flux != 0:
943 return 1.0857 * math.sqrt(flux + sky) / flux
944 else:
945
946 return (params[0] * 0.) + 99.
947
948
949
950
951
952
953 -def breakOutExt(imageObj, extNum, type='IMAGE'):
954 """Function to break out extension.
955 """
956 imageName = imageObj.fileName() if not type == 'CONF' else imageObj.confName()
957 extractName = imageName.replace(imageObj.sysc.mefType, '_%d%s' % (
958 extNum, imageObj.sysc.mefType))
959
960
961 fitsFile = fits.open(imageName)
962
963 if imageObj.isMef():
964
965 fits.combHeaders(fitsFile[0].header, fitsFile[extNum - 1].header)
966
967 if not os.path.exists(extractName):
968 phdu = pyfits.PrimaryHDU(fitsFile[extNum - 1].data)
969 phdu.header = fitsFile[0].header
970 newHDUList = pyfits.HDUList([phdu])
971 newHDUList.writeto(extractName)
972 fitsFile.close()
973 return extractName
974
975
976
977
978
980 """This function splits a header line into component parts
981 """
982 prts = line.split()
983 colNo = int(prts[1])
984 colName = prts[2]
985 units = 'Number'
986 if prts[len(prts) - 1][:1] == '[' and prts[len(prts) - 1][-1:] == ']':
987 units = prts[len(prts) - 1][1:-1]
988 noRemPrts = len(prts) - 4
989 else:
990 noRemPrts = len(prts) - 3
991 comment = ''
992 for rp in range(noRemPrts):
993 comment += prts[3 + rp] + ' '
994 return colNo, colName, units, comment
995
996
997
999 prts = line.split()
1000 nPrts = len(prts)
1001 dArray = []
1002 for pt in range(nPrts):
1003 if len(types) > pt:
1004 if types[pt] == 'int':
1005 dArray.append(int(prts[pt]))
1006 elif types[pt] == 'real':
1007 dArray.append(float(prts[pt]))
1008 else:
1009 dArray.append(prts[pt])
1010 else:
1011 try:
1012 int(prts[pt])
1013 types.append('int')
1014 dArray.append(int(prts[pt]))
1015 except:
1016 try:
1017 float(prts[pt])
1018 types.append('real')
1019 dArray.append(float(prts[pt]))
1020 except:
1021 types.append('char')
1022 dArray.append(prts[pt])
1023 return dArray, types
1024
1025
1026
1027