1
2
3 """
4 Calibration routines used by WFAU-product creating scripts and
5 re-calibration scripts.
6
7 @author: N.J.G. Cross
8 @org: WFAU, IfA, University of Edinburgh
9
10 @newfield contributors: Contributors, Contributors (Alphabetical Order)
11 @contributors:
12 """
13
14 from __future__ import division, print_function
15 from future_builtins import map, zip
16
17 from bisect import bisect
18 from collections import namedtuple
19 import numpy
20 from operator import itemgetter
21 import os
22 import math
23
24 import wsatools.CSV as csv
25 import wsatools.DbConnect.DbConstants as dbc
26 from wsatools.Logger import Logger
27 import wsatools.Statistics as stats
28
29
31 """ Reads in an illumination table and has functions to get dMag
32 """
33 xiArray = None
34 etaArray = None
35 xiInvDict = {}
36 etaInvDict = {}
37 dMagArray = None
38 isGoodTable = False
39
40 Interval = namedtuple("Interval", "x1 x2 val1 val2")
41
42
43
45 """
46 """
47
48
49
50
51
52
53 if os.path.exists(illumFileName):
54
55 self.isGoodTable = True
56 for ii, (xi, eta, dMag) in enumerate(csv.File(illumFileName)):
57 if ii == 0:
58
59 xiNo = int(xi)
60 etaNo = int(eta)
61
62 self.xiArray = numpy.zeros([xiNo], numpy.float64)
63 self.etaArray = numpy.zeros([etaNo], numpy.float64)
64 self.dMagArray = numpy.zeros([xiNo, etaNo], numpy.float64)
65 else:
66 jj = (ii - 1) // xiNo
67 kk = ii - 1 - (jj * xiNo)
68 self.dMagArray[jj][kk] = float(dMag)
69 if jj == 0:
70 self.etaArray[kk] = float(eta)
71 self.etaInvDict[float(eta)] = kk
72 if kk == 0:
73 self.xiArray[jj] = float(xi)
74 self.xiInvDict[float(xi)] = jj
75 self.coordArrays = [self.xiArray, self.etaArray]
76 self.coordInvDict = [self.xiInvDict, self.etaInvDict]
77 self.arraySize = [xiNo, etaNo]
78 self.intervals = [None, None]
79
80
81
83 """ Find interval in xi or eta
84 """
85
86 mask = numpy.greater_equal(valueArray, self.coordArrays[colNo][-1])
87
88 maskVA = numpy.where(mask == 0, valueArray, 0.)
89
90 coordArray2D = numpy.ones([self.coordArrays[colNo].size, valueArray.size])
91 coordArray2D = numpy.transpose(numpy.transpose(coordArray2D) * self.coordArrays[colNo])
92
93 lowBoundArray = numpy.nonzero(numpy.greater_equal(coordArray2D, maskVA))
94 lbArray2 = numpy.array([numpy.nonzero(numpy.equal(lowBoundArray[1], idx))[0] for idx in range(valueArray.size)])
95
96 interVal = numpy.zeros([valueArray.size],
97 dtype=[('x1', numpy.int32), ('x2', numpy.int32),
98 ('val1', numpy.float64), ('val2', numpy.float64)])
99
100 interVal['x2'] = numpy.where(mask == 0, lowBoundArray[0][numpy.array([lba[0] for lba in lbArray2])], self.arraySize[colNo] - 1)
101 x1 = numpy.where(interVal['x2'] > 0, interVal['x2'] - 1, 0)
102 interVal['x1'] = numpy.where(mask == 0, x1, self.arraySize[colNo] - 1)
103 interVal['val2'] = self.coordArrays[colNo][interVal['x2']]
104 interVal['val1'] = self.coordArrays[colNo][interVal['x1']]
105 self.intervals[colNo] = interVal
106
107
108
110 """ Find interval in xi or eta
111 """
112
113
114 if (value > self.coordArrays[colNo][self.arraySize[colNo] - 1]):
115 x1 = self.arraySize[colNo] - 1
116 x2 = self.arraySize[colNo] - 1
117 val1 = self.coordArrays[colNo][x1]
118 val2 = self.coordArrays[colNo][x2]
119 else:
120 x2 = max(numpy.nonzero(numpy.greater_equal(self.coordArrays[colNo], value))[0][0], 0)
121 val2 = self.coordArrays[colNo][x2]
122 if x2 > 0:
123 x1 = x2 - 1
124 else:
125 x1 = 0
126 val1 = self.coordArrays[colNo][x1]
127 self.intervals[colNo] = self.Interval(*[x1, x2, val1, val2])
128
129
130
131 -def _findAstromPhotomPairs(intData, deepData, isMosaic, compDict, isRecalMode,
132 isTrialRun=False, poorerZP=[]):
133 """
134 ?
135
136 @param intData: Pair of multiframe ID and dictionary of corresponding FITS
137 intermediate stack arrays (eNum, ra, dec, aMag3, aMag3Err)
138 referenced by extension number.
139 @type intData: tuple(int, dict(int: numpy.array))
140 @param deepData: List of pairs of extension numbers and corresponding
141 FITS deep stack array (ra, dec, aMag3, aMag3Err).
142 @type deepData: list(tuple(int, numpy.array))
143 @param isMosaic: If True the deep data represents a mosaic.
144 @type isMosaic: bool
145 @param isRecalMode: Gives additional logging if recalibrating
146 @type isRecalMode: bool
147 @param compDict: Dictionary of multiframeID: fileName
148 @type compDict: dict(int, str)
149 @param isTrialRun: Gives additional logging if trial run
150 @type isTrialRun: bool
151 @param poorerZP: List of lower quality calibrations
152 @type poorerZP: list(int, smallint)
153
154 @return: List of zero-point values as a tuple (bestZP, bestZPErr, medNumZP)
155 for each detector.
156 @rtype: list(tuple(real, real, int))
157
158 @todo: Maybe use pairing code.
159 @todo: Mosaic - want to use more data - many ext in int.
160
161 """
162 maxPosError = 0.5
163 deltaZP = 0.02
164 minSep = 2 / 3600
165 searchRadius = 3 / 3600
166 zeroPoints = []
167 multiframeID, intSourcesOfExt = intData
168 for extNum, deepSources in deepData:
169 if deepSources is None:
170 zeroPoints.append((dbc.realDefault(), dbc.realDefault(), 0))
171 continue
172
173 intSources = (intSourcesOfExt.get(extNum, numpy.array([[]]))
174 if not isMosaic else _getMosaicIntData(intSourcesOfExt))
175
176 intSources.dtype = [('ra', numpy.float64), ('dec', numpy.float64),
177 ('mag', numpy.float64), ('magErr', numpy.float64)]
178
179 intSources = numpy.sort(intSources, order=['dec'], axis=0)
180 photData = []
181 astromData = []
182 if intSources.size > 0:
183 decList = list(intSources['dec'])
184 for ra, dec, aMag3, _aMag3Err in deepSources:
185 minDec = dec - searchRadius
186 maxDec = dec + searchRadius
187 begIdx = max((bisect(decList, minDec) - 1), 0)
188 endIdx = min((bisect(decList, maxDec) + 1), len(decList))
189 srcSeps = numpy.hypot(intSources[begIdx:endIdx]['ra'] - ra,
190 intSources[begIdx:endIdx]['dec'] - dec)
191
192 minSrcSep = numpy.min(srcSeps)
193 nearMatch = numpy.where(minSrcSep < minSep,
194 numpy.nonzero(srcSeps == minSrcSep)[0][0], -1)
195
196 if nearMatch >= 0:
197 astromData.append(3600 * minSrcSep)
198 photData.append(
199 intSources[begIdx + nearMatch]['mag'] - aMag3)
200
201 medPos, psdPos = stats.clippedMADHistogram(astromData,
202 minSd=0.00001, retValues="med,psd")
203
204 else:
205 medPos, psdPos = dbc.realDefault(), dbc.realDefault()
206 if isTrialRun:
207 Logger.addMessage("Astrometric match: %s objects, median %s sigma %s"
208 % (len(astromData), medPos, psdPos))
209 Logger.addMessage("Photometry data: %s objects"
210 % (len(photData)))
211 if medPos > maxPosError or psdPos > maxPosError:
212 if not isRecalMode:
213 Logger.addMessage(
214 "<Info> Detector frame %s:%s has dubious astrometry "
215 "(Median = %f arcsec, stdev = %f arcsec) and will be"
216 " deprecated upon recalibration if not already."
217 % (compDict[multiframeID], extNum, medPos, psdPos),
218 alwaysLog=False)
219
220 zeroPoints.append(
221 (dbc.realDefault(), dbc.realDefault(), dbc.intDefault()))
222
223 elif medPos < 0 or psdPos < 0:
224 Logger.addMessage("<Info> Detector frame %s:%s has dubious "
225 "astrometry (no matches) and will be deprecated upon "
226 "recalibration if not already." % (compDict[multiframeID],
227 extNum))
228
229 zeroPoints.append((dbc.realDefault(), dbc.realDefault(), 0))
230
231 else:
232 medZPArr = stats.clippedMADHistogram(numpy.array(photData),
233 minSd=0.0001, retValues="med,psd")
234 zp = medZPArr[0][0]
235 zpErr = medZPArr[1][0]
236 if (multiframeID, extNum) in poorerZP:
237
238 zpErr = math.sqrt(zpErr * zpErr + deltaZP * deltaZP)
239 zeroPoints.append((zp, zpErr, len(photData)))
240 if isTrialRun:
241 Logger.addMessage("Photometry ZP: %s +/- %s" % (zeroPoints[-1][0], zeroPoints[-1][1]))
242 return zeroPoints
243
244
245
247 """
248 Produces array of intermediate data to compare to deep mosaic data.
249
250 @param intSourcesOfExt: Dictionary of FITS intermediate stack arrays
251 (eNum, ra, dec, aMag3, aMag3Err) referenced by
252 extension number.
253 @type intSourcesOfExt: dict(int: numpy.array)
254
255 @return: Array of intermediate data to compare to deep mosaic data.
256 @rtype: numpy.array
257
258 """
259 extNums = sorted(set(intSourcesOfExt))
260 xSizes = [len(intSourcesOfExt[extNum]) for extNum in extNums]
261 ySizes = [len(intSourcesOfExt[extNum][0])
262 for extNum in extNums if len(intSourcesOfExt[extNum]) > 0]
263
264 ySize = ySizes[0]
265 intSources = numpy.zeros([sum(xSizes), ySize])
266 xMax = 0
267 for xSize, extNum in zip(xSizes, extNums):
268 xMin = xMax
269 xMax += xSize
270 intSources[xMin : xMax, : ySize] = \
271 intSourcesOfExt.get(extNum, numpy.array([[]]))
272
273 return intSources
274
275
276
278 """ Calculate illumination correction
279 """
280
281 illum = IllumTable(illumTableData)
282
283 if illum.isGoodTable:
284 illum.setIntervalArray(xiArray, 0)
285 illum.setIntervalArray(etaArray, 1)
286
287 illumCorr = 0.0
288 if not ('deep' in frameType or 'tile' in frameType or 'mosaic' in frameType):
289
290 xiIntervall = illum.intervals[0]
291 etaIntervall = illum.intervals[1]
292
293
294 illumCorr = numpy.zeros([xiArray.size])
295 d1 = xiIntervall['val2'] - xiIntervall['val1']
296 d2 = etaIntervall['val2'] - etaIntervall['val1']
297
298
299 mask = numpy.logical_and(numpy.equal(d1, 0.), numpy.not_equal(d2, 0.))
300 d2msk = numpy.where(mask, d2, 1.)
301 f1 = illum.dMagArray[xiIntervall['x1'], etaIntervall['x1']]
302 f2 = illum.dMagArray[xiIntervall['x1'], etaIntervall['x2']] - f1
303 f2 *= (etaArray - etaIntervall['val1']) / d2msk
304 illumCorr = numpy.where(mask, f1 + f2, illumCorr)
305
306 mask = numpy.logical_and(numpy.equal(d2, 0.), numpy.not_equal(d1, 0.))
307 d1msk = numpy.where(mask, d1, 1.)
308 f1 = illum.dMagArray[xiIntervall['x1'], etaIntervall['x1']]
309 f2 = illum.dMagArray[xiIntervall['x2'], etaIntervall['x1']] - f1
310 f2 *= (xiArray - xiIntervall['val1']) / d1msk
311 illumCorr = numpy.where(mask, f1 + f2, illumCorr)
312
313 mask = numpy.logical_and(numpy.equal(d1, 0.), numpy.equal(d2, 0.))
314 illumCorr = numpy.where(mask, illum.dMagArray[xiIntervall['x1'], etaIntervall['x1']], illumCorr)
315
316 mask = numpy.logical_and(numpy.not_equal(d1, 0.), numpy.not_equal(d2, 0.))
317 d1msk = numpy.where(mask, d1, 1.)
318 d2msk = numpy.where(mask, d2, 1.)
319 f1 = illum.dMagArray[xiIntervall['x1'], etaIntervall['x1']]
320 f1 *= (xiIntervall['val2'] - xiArray)
321 f2 = illum.dMagArray[xiIntervall['x2'], etaIntervall['x1']]
322 f2 *= (xiArray - xiIntervall['val1'])
323 ic = (f1 + f2) * (etaIntervall['val2'] - etaArray)
324 f3 = illum.dMagArray[xiIntervall['x1'], etaIntervall['x2']]
325 f3 *= (xiIntervall['val2'] - xiArray)
326 f4 = illum.dMagArray[xiIntervall['x2'], etaIntervall['x2']]
327 f4 *= (xiArray - xiIntervall['val1'])
328 ic += (f3 + f4) * (etaArray - etaIntervall['val1'])
329 ic /= (d1msk * d2msk)
330 illumCorr = numpy.where(mask, ic, illumCorr)
331 else:
332 raise SystemExit("Missing illumination table")
333 return illumCorr
334
335
336
337 -def checkAstrometryPhotometry(allIntData, deepData, isMosaic, compDict,
338 isRecalMode=False, isTrialRun=False, poorerZP=[]):
339 """
340 Tests whether astrometry and photometry is good enough.
341
342 @param allIntData: List of pairs of multiframe IDs and dictionaries of
343 corresponding intermediate stack source measurements
344 (ra, dec, aMag3, aMag3Err) referenced by extension
345 number.
346 @type allIntData: list(tuple(int, dict(int: numpy.array)))
347 @param deepData: List of pairs of extension numbers and corresponding
348 deep stack source measurements
349 (ra, dec, aMag3, aMag3Err) ordered by extension number.
350 @type deepData: list(tuple(int, numpy.array))
351 @param isMosaic: If True, the deep data represents a mosaic.
352 @type isMosaic: bool
353 @param compDict: Dictionary of multiframeID: fileName
354 @type compDict: dict(int, str)
355 @param isRecalMode: If True, recalibrating intermediates, else finding ZP
356 for deep.
357 @type isRecalMode: bool
358 @param isTrialRun: If True, log additional test info,.
359 @type istrialRun: bool
360 @param poorerZP: List of lower quality calibrations
361 @type poorerZP: list(int, smallint)
362
363 @return: Dictionary of zero-points (bestZP, bestZPErr, medNumZP) referenced
364 by extension number.
365 @rtype: dict(int: tuple(real, real, int))
366
367 """
368 dataZP = numpy.zeros([len(deepData), len(allIntData)], numpy.float32)
369 dataZPErr = numpy.zeros([len(deepData), len(allIntData)], numpy.float32)
370 numZP = numpy.zeros([len(deepData), len(allIntData)], numpy.int32)
371
372
373 for j, intData in enumerate(allIntData):
374
375 zeroPoints = _findAstromPhotomPairs(intData, deepData, isMosaic,
376 compDict, isRecalMode, isTrialRun,
377 poorerZP)
378 for i, zp in enumerate(zeroPoints):
379 dataZP[i][j], dataZPErr[i][j], numZP[i][j] = zp
380 if isTrialRun:
381 Logger.addMessage("ZP data array shape: %s x %s" % dataZP.shape)
382 newZPDict = {}
383 for i, extNum in enumerate(map(itemgetter(0), deepData)):
384
385 goodIndices = [ind for ind, zp in enumerate(dataZP[i]) if zp > 0]
386 gdNumZP = [zp for ind, zp in enumerate(numZP[i]) if ind in goodIndices]
387 gdDataZP = \
388 [zp for ind, zp in enumerate(dataZP[i]) if ind in goodIndices]
389
390 gdDataZPErr = \
391 [zp for ind, zp in enumerate(dataZPErr[i]) if ind in goodIndices]
392 if isTrialRun:
393 Logger.addMessage("There are %s good indices: %s" %
394 (len(gdDataZP), ','.join(map(str, goodIndices))))
395
396 if len(gdDataZP) > 1:
397 newZPDict[extNum] = \
398 stats.clippedWeightedMean(gdDataZP, gdDataZPErr, 3.) + (int(
399 stats.clippedMADHistogram(gdNumZP, retValues="median")[0]),)
400
401 elif len(gdDataZP) == 1:
402 newZPDict[extNum] = (gdDataZP[0], gdDataZPErr[0], gdNumZP[0])
403 else:
404 newZPDict[extNum] = (dbc.realDefault(), dbc.realDefault(),
405 dbc.intDefault())
406 if isTrialRun:
407 Logger.addMessage("ExtNum %s, ZP: %s+/-%s, number of calibrators %s"
408 % (extNum, newZPDict[extNum][0], newZPDict[extNum][1], newZPDict[extNum][2]))
409 if newZPDict[extNum][0] < 0:
410 newZPDict[extNum] = (dbc.realDefault(), dbc.realDefault(),
411 dbc.intDefault())
412
413 return newZPDict
414
415
416
418 """
419 Calculate the photometry for the multiframe.
420
421 @param intStackPhotom: Database query results as a numpy array.
422 @type intStackPhotom: numpy.array(tuple)
423
424 """
425 photometryOfExt = {}
426 if poorPhotList:
427 poorStackPhotom, poorStackCounts = poorPhotList
428 else:
429 poorStackPhotom, poorStackCounts = [], []
430
431 for extNum in extNums:
432 photIdx = [idx for idx, (cpID, eNum, _numb) in enumerate(intPhotCounts)
433 if cpID == compID and eNum == extNum][0]
434 number = intPhotCounts[photIdx][2]
435 index = sum(intPhotCounts[idx][2] for idx in range(photIdx))
436 if number >= 20:
437 photometryOfExt[extNum] = numpy.array(
438 intStackPhotom[index:index + number, 2:], dtype=numpy.float64)
439 else:
440 pPhotIdx = [idx for idx, (cpID, eNum, _numb) in enumerate(poorStackCounts)
441 if cpID == compID and eNum == extNum][0]
442 begIdx = sum(poorStackCounts[idx][2] for idx in range(pPhotIdx))
443 pnumber = poorStackCounts[pPhotIdx][2]
444 endIdx = begIdx + pnumber
445 if pnumber < 20:
446 Logger.addMessage("Not enough matches for a decent ZP in component %s, extNum %s" % (compID, extNum))
447 photometryOfExt[extNum] = numpy.array(
448 poorStackPhotom[begIdx:endIdx, 2:], dtype=numpy.float64)
449 return photometryOfExt
450
451
452