Package wsatools :: Module Calibration
[hide private]

Source Code for Module wsatools.Calibration

  1  #------------------------------------------------------------------------------ 
  2  # $Id: Calibration.py 10153 2013-12-06 11:14:03Z NicholasCross $ 
  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   
30 -class IllumTable(object):
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
44 - def __init__(self, illumFileName):
45 """ 46 """ 47 # @TODO: IllumFile may actually be csv file --> not all there yet? 48 # Does ETWS create csv files from illumFile 49 # +.csv 50 # first line - nXi, nEta, filterName 51 # then xi,eta,dmag 52 # @TODO: Is this the most efficient way of doing it? 53 if os.path.exists(illumFileName): 54 # xi repeated... 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 # filterName = dMag 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
82 - def setIntervalArray(self, valueArray, colNo):
83 """ Find interval in xi or eta 84 """ 85 # @TODO: change greater to greater_equal? 86 mask = numpy.greater_equal(valueArray, self.coordArrays[colNo][-1]) 87 # @TODO: output = numpy.array([valueArray.size(),4]) 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 # Set last value in greater_equal array to True (so at least one non_zero value 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
109 - def setInterval(self, value, colNo):
110 """ Find interval in xi or eta 111 """ 112 # @TODO: Work with arrays... 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 # arcsec 163 deltaZP = 0.02 # additional ZP error if poorer ZP 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 # Add in quadrature 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
246 -def _getMosaicIntData(intSourcesOfExt):
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
277 -def calcIllumCorrArray(illumTableData, xiArray, etaArray, frameType):
278 """ Calculate illumination correction 279 """ 280 281 illum = IllumTable(illumTableData) 282 # @TODO: convert to numpy arrays 283 if illum.isGoodTable: 284 illum.setIntervalArray(xiArray, 0) 285 illum.setIntervalArray(etaArray, 1) 286 # Illumination correction 287 illumCorr = 0.0 288 if not ('deep' in frameType or 'tile' in frameType or 'mosaic' in frameType): 289 # don't apply to deep images or mosaics 290 xiIntervall = illum.intervals[0] 291 etaIntervall = illum.intervals[1] 292 # Initialise 293 # f1, f2, f3, f4 = numpy.zeros([xiArray.size]), numpy.zeros([xiArray.size]), numpy.zeros([xiArray.size]), numpy.zeros([xiArray.size]) 294 illumCorr = numpy.zeros([xiArray.size]) 295 d1 = xiIntervall['val2'] - xiIntervall['val1'] 296 d2 = etaIntervall['val2'] - etaIntervall['val1'] 297 # create mask... 298 # case 1 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 # case 2 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 # case 3 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 # case 4 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 # First get average astrometry for input frames. 373 for j, intData in enumerate(allIntData): 374 # Find pairs, but also get int stack offsets. 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 # remove defaults 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
417 -def getIntermedPhotometry(compID, intPhotCounts, intStackPhotom, extNums, poorPhotList=None):
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