1
2
3 """
4 Astrometric calculations. Routines for positional calculations in the
5 celestial sphere derived from FORTRAN SLALIB of P.T. Wallace. Functions to
6 convert between pixel (X, Y) co-ordinates in the detector plane and sky
7 (RA, Dec) co-ordinates in the celestial sphere, use our PyWCS wrapper to the
8 WCS routines of Starlink's AST library.
9
10 Pairing module replacement functions::
11 import wsatools.pairing as pairing
12 pairing.multiPair(inputList, matchRadius, workPath)
13 pairing.pairOff(set0, set1, pFile0, pFile1, pairCriterion)
14 becomes::
15 import wsatools.Astrometry as astro
16 astro.pairObjects(inputList, matchRadius, workPath)
17 astro.pairObjects([set0, set1], pairCriterion, pFile0, oneWay=True)
18
19 @author: R.S. Collins, N.C. Hambly
20 @org: WFAU, IfA, University of Edinburgh
21
22 @requires: Starlink-AST, Numpy, PyWCS
23 """
24
25 from __future__ import division, print_function
26 from future_builtins import zip
27
28 from bisect import bisect
29 from collections import namedtuple
30 import ctypes
31 import math
32 import numpy
33 import os
34 import struct
35
36 import wsatools.DbConnect.DbConstants as dbc
37 from wsatools.DbConnect.DbSession import DbSession, Join
38 import wsatools.pywcs as pywcs
39 import wsatools.Statistics as stats
40
41
42
43
44 Pixel = namedtuple("Pixel", "x y")
45
46
47 SkyPos = namedtuple("SkyPos", "ra dec")
48
49
50
51 -def angDiff(angle1, angle2, minVal= -180, cycleLength=360):
52 """
53 Calculates the angular difference between the given two angles accounting
54 for zero crossings.
55
56 @param angle1: First angle (default is degrees).
57 @type angle1: float
58 @param angle2: Second angle (default is degrees).
59 @type angle2: float
60 @param minVal: Value of the angle at the "zero crossing".
61 @type minVal: float
62 @param cycleLength: Difference between minimum and maximum values, i.e.
63 separation of values at "zero crossing".
64 @type cycleLength: float
65
66 @return: Angular difference between two angles.
67 @rtype: float
68
69 """
70 delta = angle2 - angle1
71 if delta < minVal:
72 delta += cycleLength
73 elif delta > minVal + cycleLength:
74 delta -= cycleLength
75
76 return delta
77
78
79
81 """
82 Angle between two points on a sphere.
83
84 @note: Based on the SLALIB algorithms dsep() and dsepv() by P.T. Wallace.
85
86 @param pos1: Co-ordinate on celestial sphere in radians of the first point,
87 (RA, Dec), though can be any spherical co-ordinate.
88 @type pos1: tuple(float, float)
89 @param pos2: Co-ordinate on celestial sphere in radians of the second
90 point, (RA, Dec), though can be any spherical co-ordinate.
91 @type pos2: tuple(float, float)
92
93 @return: Angle in radians (always positive) between the two points.
94 @rtype: float
95
96 """
97 v1 = convSkyToCartesian(pos1)
98 v2 = convSkyToCartesian(pos2)
99
100
101
102
103
104
105
106
107
108
109
110
111 v1xv2 = numpy.cross(v1, v2)
112 modulus = numpy.sqrt(numpy.vdot(v1xv2, v1xv2))
113 return math.atan2(modulus, numpy.vdot(v1, v2)) if modulus else 0.0
114
115
116
118 """
119 Angle between two arrays of points on a sphere or plane.
120
121 @note: Based on the SLALIB algorithms dsep() and dsepv() by P.T. Wallace.
122
123 @param v1: Co-ordinates in radians of the first points, e.g. (RA, Dec).
124 @type v1: numpy.array
125 @param v2: Co-ordinates in radians of the second points, e.g. (RA, Dec).
126 @type v2: numpy.array
127 @param isEqu: True, if equatorial (spherical/celestial) co-ordinates.
128 @type isEqu: bool
129
130 @return: Angle in radians (always positive) between the two points.
131 @rtype: float
132
133 """
134 if isEqu:
135
136 v1 = convSkyToCartesian(v1)
137 v2 = convSkyToCartesian(v2)
138
139
140
141
142
143
144
145
146
147
148
149
150 v1xv2 = numpy.cross(v1, v2)
151 modulus = \
152 numpy.sqrt(numpy.diagonal(numpy.dot(v1xv2, numpy.transpose(v1xv2))))
153
154 return numpy.where(numpy.greater(modulus, 0),
155 numpy.arctan2(modulus,
156 numpy.diagonal(numpy.dot(v1, numpy.transpose(v2)))),
157 0)
158
159
160
162 """
163 Takes list of positions in radians and returns centre in radians.
164
165 @param cartPosList: Array of RA, Dec co-ordinate positions in radians.
166 @type cartPosList: numpy.array
167 @param isWeighted: If True, calculate using weights.
168 @type isWeighted: bool
169
170 @return: Centre of list of positions in RA, Dec co-ordinates with errors
171 (radians).
172 @rtype: tuple(float, float, float, float)
173
174 """
175 if isWeighted:
176 mX, errX = \
177 stats.clippedWeightedMean(cartPosList[:, 0], cartPosList[:, 3], 3.)
178
179 mY, errY = \
180 stats.clippedWeightedMean(cartPosList[:, 1], cartPosList[:, 3], 3.)
181
182 mZ, errZ = \
183 stats.clippedWeightedMean(cartPosList[:, 2], cartPosList[:, 3], 3.)
184
185 nDof = cartPosList.size - 3
186 chiSq = numpy.sum(
187 numpy.power(((cartPosList[:, 0] - mX) / cartPosList[:, 3]), 2)
188 + numpy.power(((cartPosList[:, 1] - mY) / cartPosList[:, 3]), 2)
189 + numpy.power(((cartPosList[:, 2] - mZ) / cartPosList[:, 3]), 2)) \
190 / nDof
191
192 expErrSize = 1 \
193 / numpy.sqrt(numpy.sum(1 / (cartPosList[:, 3] * cartPosList[:, 3])))
194
195 expAstErr = numpy.hypot(expErrSize, expErrSize)
196
197 else:
198 mX, errX = stats.clippedHistogram(cartPosList[:, 0])
199 mY, errY = stats.clippedHistogram(cartPosList[:, 1])
200 mZ, errZ = stats.clippedHistogram(cartPosList[:, 2])
201
202
203
204
205
206
207 isNormalCase = abs(mX) < 1 and abs(mY) <= 1 and abs(mZ) < 1
208 if abs(mZ) == 1:
209
210 mRa, mDec = 0, mZ * math.pi / 2
211 errXd, errYd = errX, errY
212
213 elif isNormalCase or mX == 0:
214 mRa = (math.atan2(mY, mX) if isNormalCase else mY)
215 if -2 * math.pi < mRa < 0:
216 mRa += 2 * math.pi
217
218 mDec = math.asin(mZ)
219 xyC = math.hypot(mX, mY)
220 errXd = math.hypot(mY * errX / xyC, mX * errY / xyC)
221 errY2d = math.hypot(mX * mZ * errX / xyC, mY * mZ * errY / xyC)
222 errYd = math.hypot(xyC * errZ, errY2d)
223
224 else:
225 mRa = mDec = errXd = errYd = dbc.realDefault()
226
227 if isWeighted:
228 return mRa, mDec, errXd, errYd, chiSq, expAstErr
229 else:
230 return mRa, mDec, errXd, errYd
231
232
233
235 """
236 Calculates hSpaces for a set of coordinates
237
238 @param coords: List of Cartesian coordinates on the sky
239 @type coords: list(tuple)
240 @param initPos: Initial positions of each new hspace in coords
241 @type initPos: list(int)
242 @param nCoordInHSDef: Number of coordinates used to define an H-space
243 @type nCoordInHSDef: int
244
245 @return: list of values, 4 for each half-space. The first 3 are the
246 Cartesian coordinates of the normal vector to the plane and the
247 4th is the offset of the plane from the origin.
248 @rtype: list
249
250 """
251 hSpaces = []
252 for index in initPos:
253 hsCoords = coords[index:index + nCoordInHSDef]
254 if nCoordInHSDef > 2:
255 planeCoords = hsCoords[1:3] - hsCoords[:2]
256 else:
257 planeCoords = hsCoords[:2]
258
259 norm = numpy.cross(planeCoords[0], planeCoords[1])
260 norm /= numpy.sqrt(numpy.sum(numpy.square(norm)))
261 hSpaces.extend(norm)
262 hSpaces.append(numpy.dot(hsCoords[0], norm))
263
264 return hSpaces
265
266
267
269 """
270 Calculate an HTM index for the given position at the given level.
271
272 >>> calcHtmIndex(ra=200.981, dec=-1.06456, level=20)
273 11030377854233L
274
275 @param ra: Right ascension in degrees.
276 @type ra: float
277 @param dec: Declination in degrees.
278 @type dec: float
279 @param level: HTM level to calculate down to.
280 @type level: int
281
282 @return: HTM index of given position.
283 @rtype: long
284
285 """
286 libhtm = ctypes.cdll.LoadLibrary("libhtm.so")
287 libhtm.calcHtmId.restype = ctypes.c_ulonglong
288 ra = ctypes.c_double(ra)
289 dec = ctypes.c_double(dec)
290 level = ctypes.c_int(level)
291
292 return libhtm.calcHtmId(ra, dec, level)
293
294
295
297 """
298 Takes a set of required mosaics and calculates a set of positions for the
299 corners.
300
301 @param requiredMosaics: List of namedtuples of RequiredMosaics attributes
302 @type requiredMosaics: list(namedtuple)
303
304 @return: Dictionary of 2D numpy arrays of Cartesian coordinates for the
305 corners of the mosaic of the given productID
306 @rtype: dict(int: numpy.array)
307
308 """
309
310
311
312 corners = [(0, -1), (0, 1), (1, -1), (1, 1)]
313 mosaicCornerPos = dict()
314 for mosaic in requiredMosaics:
315 angCorner = math.atan(mosaic.decExtent / mosaic.raExtent)
316 distCorner = math.radians(math.hypot(mosaic.raExtent, mosaic.decExtent))
317 distCorner /= 2
318 coords = []
319 for multFact, addFact in corners:
320 cornDirect = math.radians(mosaic.posAngle) + multFact * math.pi
321 cornDirect += addFact * angCorner
322 dDelta = math.asin(math.sin(distCorner) * math.cos(cornDirect))
323 dAlpha = math.asin(math.sin(distCorner) * math.sin(cornDirect))
324 alpha = math.radians(mosaic.ra) + dAlpha
325 delta = math.radians(mosaic.dec) + dDelta
326 coords.append((alpha, delta))
327
328 cartCoords = convSkyToCartesian(numpy.array(coords))
329 mosaicCornerPos[mosaic.productID] = cartCoords
330
331 return mosaicCornerPos
332
333
334
336 """
337 Converts a J2000 (FK5) Cartesian vector to a Galactic coordinate.
338
339 @param skyVector: 3-vector of Cartesian coordinates (FK5 J2000).
340 @type skyVector: numpy.array(float, float, float)
341
342 @return: Galactic coordinate (l, b) in degrees.
343 @rtype: tuple(float, float)
344
345 """
346 galRotMatrix = numpy.array([
347 [-0.054875539726, -0.873437108010, -0.483834985808],
348 [0.494109453312, -0.444829589425, 0.746982251810],
349 [-0.867666135858, -0.198076386122, 0.455983795705]])
350
351 return convCartesianToSky(numpy.dot(galRotMatrix, skyVector))
352
353
354
356 """
357 Converts Cartesian vector to a spherical coordinate. The spherical
358 coordinate is longitude (+ve anticlockwise looking from the positive
359 latitude pole) and latitude. The Cartesian coordinate is right handed,
360 with the x-axis at zero longitude and latitude, and the z-axis at the
361 positive latitude pole.
362
363 @param skyVector: 3-vector of right-handed Cartesian coordinates.
364 @type skyVector: tuple(float, float, float)
365
366 @return: Coordinate on celestial sphere in degrees (RA, Dec),
367 though can be any spherical coordinate.
368 @rtype: tuple(float, float)
369
370 """
371 cx, cy, cz = skyVector
372 isNormalCase = abs(cx) < 1 and abs(cy) <= 1 and abs(cz) < 1
373 if abs(cz) == 1:
374
375 return 0, math.degrees(cz * math.pi / 2)
376
377 elif isNormalCase or cx == 0:
378 ra = math.atan2(cy, cx) if isNormalCase else cy
379 if -2 * math.pi < ra < 0:
380 ra += 2 * math.pi
381
382 dec = math.asin(cz)
383
384 return math.degrees(ra), math.degrees(dec)
385
386 return dbc.realDefault(), dbc.realDefault()
387
388
389
390 -def convPixelToCartSky(pixels, multiframeID=None, extNum=None, db=None,
391 wcs=None, projection="ZPN"):
392 """
393 Converts from pixel (X,Y) co-ordinates in the detector plane to Cartesian
394 sky co-ordinates in the celestial sphere, for a given detector's
395 observation (e.g. multiframeID, extNum). If multiframeID and extNum
396 supplied it looks up WCS values in the database for this observation,
397 otherwise WCS must be supplied. Uses PyWCS to perform the conversion.
398
399 @param pixels: List of pixel (X,Y) co-ordinates.
400 @type pixels: list(tuple(float, float))
401 @param multiframeID: MultiframeID of the observation.
402 @type multiframeID: int
403 @param extNum: Detector extension number.
404 @type extNum: int
405 @param db: Optional database connection to obtain WCS
406 co-ordinates.
407 @type db: DbSession
408 @param wcs: World Co-ordinate System for detector.
409 @type wcs: tuple(float)
410 @param projection: Projection, either ZPN or TAN.
411 @type projection: str
412
413 @return: Array of 3-vectors of right-handed Cartesian sky coordinates
414 ordered by and corresponding to the list of pixel co-ordinates.
415 @rtype: numpy.array(numpy.array(float, float, float))
416
417 @todo: Would be more efficient if PyWCS performed multiple co-ordinate
418 conversions for a given WCS value set, using ctypes.
419 """
420 wcs = wcs or _queryWCS(db or DbSession(), multiframeID, extNum)
421 coords = [pywcs.convPixelToSky(x, y, projection, wcs) for x, y in pixels]
422
423 return convSkyToCartesian(numpy.array(coords))
424
425
426
427 -def convPixelToSky(pixels, multiframeID=None, extNum=None, db=None, wcs=None,
428 projection="ZPN"):
429 """
430 Converts from pixel (X,Y) co-ordinates in the detector plane to sky
431 (RA, Dec) co-ordinates in the celestial sphere, for a given detector's
432 observation (e.g. multiframeID, extNum). If multiframeID and extNum
433 supplied it looks up WCS values in the database for this observation,
434 otherwise WCS must be supplied. Uses PyWCS to perform the conversion.
435
436 @param pixels: List of pixel (X,Y) co-ordinates.
437 @type pixels: list(tuple(float, float))
438 @param multiframeID: MultiframeID of the observation.
439 @type multiframeID: int
440 @param extNum: Detector extension number.
441 @type extNum: int
442 @param db: Optional database connection to obtain WCS
443 co-ordinates.
444 @type db: DbSession
445 @param wcs: World Co-ordinate System for detector.
446 @type wcs: tuple(float)
447 @param projection: Projection, either ZPN or TAN.
448 @type projection: str
449
450 @return: Sky (RA, Dec) co-ordinates in degrees ordered by and corresponding
451 to the list of pixel co-ordinates.
452 @rtype: list(SkyPos)
453
454 @todo: Would be more efficient if PyWCS performed multiple co-ordinate
455 conversions for a given WCS value set, using ctypes.
456 """
457 wcs = wcs or _queryWCS(db or DbSession(), multiframeID, extNum)
458
459 coords = [pywcs.convPixelToSky(x, y, projection, wcs) for x, y in pixels]
460
461 return [SkyPos(math.degrees(normalAngle(ra)), math.degrees(dec))
462 for ra, dec in coords]
463
464
465
466 -def convSkyToPixel(positions, multiframeID=None, extNum=None, db=None,
467 wcs=None, projection="ZPN"):
468 """
469 Converts from sky (RA, Dec) co-ordinates in the celestial sphere to pixel
470 (X,Y) co-ordinates in the detector plane, for a given detector's
471 observation (e.g. multiframeID, extNum). If multiframeID and extNum
472 supplied it looks up WCS values in the database for this observation,
473 otherwise WCS must be supplied. Uses PyWCS to perform the conversion.
474
475 @param positions: List of sky (RA, Dec) co-ordinates in degrees.
476 @type positions: list(tuple(float, float))
477 @param multiframeID: MultiframeID of the observation.
478 @type multiframeID: int
479 @param extNum: Detector extension number.
480 @type extNum: int
481 @param db: Optional database connection to obtain WCS
482 co-ordinates.
483 @type db: DbSession
484 @param wcs: World Co-ordinate System for detector.
485 @type wcs: tuple(float)
486 @param projection: Projection, either ZPN or TAN.
487 @type projection: str
488
489 @return: Pixel (X, Y) co-ordinates ordered by and corresponding to the list
490 of sky position co-ordinates.
491 @rtype: list(Pixel)
492
493 @todo: Would be more efficient if PyWCS performed multiple co-ordinate
494 conversions for a given WCS value set, using ctypes.
495 """
496 wcs = wcs or _queryWCS(db or DbSession(), multiframeID, extNum)
497 return [Pixel(*pywcs.convSkyToPixel(math.radians(ra), math.radians(dec),
498 projection, wcs))
499 for ra, dec in positions]
500
501
502
504 """
505 Converts spherical coordinate(s) to Cartesian vector(s). The spherical
506 coordinates are longitude (+ve anticlockwise looking from the positive
507 latitude pole) and latitude. The Cartesian coordinates are right handed,
508 with the x-axis at zero longitude and latitude, and the z-axis at the
509 positive latitude pole.
510
511 @note: Based on the SLALIB algorithm dcs2c() by P.T. Wallace.
512
513 @param coords: Co-ordinate(s) on celestial sphere in radians (RA, Dec),
514 though can be any spherical co-ordinate.
515 @type coords: tuple(float, float) or numpy.array(numpy.array(float, float))
516
517 @return: (Array of) 3-vector(s) of right-handed Cartesian coordinates.
518 @rtype: numpy.array(float, float, float)
519 or numpy.array(numpy.array(float, float, float))
520
521 """
522 if not isinstance(coords, numpy.ndarray) or coords.ndim == 1:
523 ra, dec = coords
524 cosDec = math.cos(dec)
525
526 return numpy.array([math.cos(ra) * cosDec,
527 math.sin(ra) * cosDec,
528 math.sin(dec)])
529
530 cosDec = numpy.cos(coords[:, 1])
531 cartCoords = numpy.zeros([len(coords), 3])
532 cartCoords[:, 0] = numpy.cos(coords[:, 0]) * cosDec
533 cartCoords[:, 1] = numpy.sin(coords[:, 0]) * cosDec
534 cartCoords[:, 2] = numpy.sin(coords[:, 1])
535
536 return cartCoords
537
538
539
541 """
542 Converts B1950 (FK4) sky coordinate to Galactic coordinate.
543
544 @note: Based on Practical Astronomy With Your Calculator. Section 29 p43.
545
546 @param coords: Co-ordinate on celestial sphere in radians (RA, Dec).
547 @type coords: tuple(float, float)
548
549 @return: Galactic coordinate (l, b) in degrees.
550 @rtype: tuple(float, float)
551
552 """
553 gRa, gDec = math.radians(192.25), math.radians(27.4)
554 cosGDec = math.cos(gDec)
555 sinGDec = math.sin(gDec)
556 ra, dec = coords
557 ra -= gRa
558 cosDec = math.cos(dec)
559 sinDec = math.sin(dec)
560 sinB = cosDec * cosGDec * math.cos(ra) + sinDec * sinGDec
561 l = math.radians(33)
562 l += math.atan2(sinDec - sinB * sinGDec, cosDec * math.sin(ra) * cosGDec)
563 b = math.asin(sinB)
564
565 return math.degrees(l), math.degrees(b)
566
567
568
570 """
571 Projection of spherical coordinates onto the tangent plane, also
572 known as the gnomic projection. Function translated from the SLALIB FORTRAN
573 source.
574
575 @note: Based on the SLALIB algorithm ds2tp() by P.T. Wallace.
576
577 @param pos: Co-ordinate on celestial sphere in radians to be projected
578 (RA, Dec), though can be any spherical co-ordinate.
579 @type pos: tuple(float, float)
580 @param axis: Co-ordinate on celestial sphere in radians of the tangent
581 point (RA, Dec), though can be any spherical co-ordinate.
582 @type axis: tuple(float, float)
583
584 @return: Tangent plane co-ordinates (Xi, Eta) in radians.
585 @rtype: float, float
586
587 """
588 ra, dec = pos
589 raZ, decZ = axis
590
591
592 sDec, sDecZ = math.sin(dec), math.sin(decZ)
593 cDec, cDecZ = math.cos(dec), math.cos(decZ)
594 raDiff = ra - raZ
595 sRaDiff, cRaDiff = math.sin(raDiff), math.cos(raDiff)
596
597
598 denom = sDec * sDecZ + cDec * cDecZ * cRaDiff
599
600
601 tiny = 1.0e-6
602
603 if denom <= tiny:
604 if denom >= 0.0:
605 denom = tiny
606 elif denom > -tiny:
607 denom = -tiny
608
609
610
611
612 xi = cDec * sRaDiff / denom
613 eta = (sDec * cDecZ - cDec * sDecZ * cRaDiff) / denom
614 return xi, eta
615
616
617
619 """
620 Transform tangent plane coordinates back to spherical given the
621 gnomic projection (standard) and tangent point coordinates. Function
622 translated from the SLALIB Fortran source.
623
624 @note: Based on the SLALIB algorithm dtp2s() by P.T. Wallace.
625
626 @param pos: Co-ordinate on tangent plane in radians (Xi, Eta).
627 @type pos: tuple(float, float)
628 @param axis: Co-ordinate on celestial sphere in radians of the tangent
629 point (RA, Dec), though can be any spherical co-ordinate.
630 @type axis: tuple(float, float)
631
632 @return: Co-ordinate on celestial sphere in radians (RA, Dec), though can
633 be any spherical co-ordinate.
634 @rtype: float, float
635
636 """
637 xi, eta = pos
638 raZ, decZ = axis
639 sDecZ, cDecZ = math.sin(decZ), math.cos(decZ)
640 denom = cDecZ - eta * sDecZ
641 ra = normalAngle(math.atan2(xi, denom) + raZ)
642 dec = math.atan2(sDecZ + eta * cDecZ, math.sqrt(xi * xi + denom * denom))
643
644 return ra, dec
645
646
647
649 """
650 Calculates the sky (RA, Dec) co-ordinates of the given detector's corners.
651
652 @param multiframeID: MultiframeID of the observation.
653 @type multiframeID: int
654 @param extNum: Detector extension number.
655 @type extNum: int
656 @param db: Database connection to use for metadata queries.
657 @type db: DbSession
658
659 @return: Sky (RA, Dec) co-ordinates in degrees of each corner.
660 @rtype: list(SkyPos)
661
662 """
663 imageShape = db.query(
664 selectStr="axis1Length, axis2Length, tctyp3",
665 fromStr=Join([("MultiframeDetector", 'D'), ("CurrentAstrometry", 'A')],
666 ["multiframeID", "extNum"]),
667 whereStr="D.multiframeID=%s AND D.extNum=%s" % (multiframeID, extNum),
668 firstOnly=True)
669
670 pixels = [(1, 1), (1, imageShape.axis2Length),
671 (imageShape.axis1Length, imageShape.axis2Length),
672 (imageShape.axis1Length, 1)]
673
674 proj = imageShape.tctyp3[-3:]
675
676 return convPixelToSky(pixels, multiframeID, extNum, db, projection=proj)
677
678
679
681 """ Takes coordinates of pawprint, offsetPos and posAngle.
682 """
683 offsetX, offsetY = getOffSetsFromPos(offsetPos)
684
685
686 offRa = math.radians((offsetX * math.cos(math.radians(posAngle))
687 + offsetY * math.sin(math.radians(posAngle))) / 3600)
688
689 offDec = math.radians((-offsetX * math.sin(math.radians(posAngle))
690 + offsetY * math.cos(math.radians(posAngle))) / 3600)
691
692 return convTangentToSky((offRa, offDec), offsetCoords)
693
694
695
696 -def getMedAngle(angles, expRange, minVal= -180, cycleLength=360):
697 """ @return: Median angle, taking into account cyclical nature of
698 measurement.
699 @rtype: float
700 """
701 if len(angles) is 1:
702 return angles[0]
703
704
705 if max(angles) > 0.5 * cycleLength - 2 * expRange \
706 and min(angles) < -0.5 * cycleLength + 2 * expRange:
707
708 angles = [(angle + cycleLength if angle < 0 else angle)
709 for angle in angles]
710
711 median = float(stats.clippedMADHistogram(angles, retValues="med")[0])
712 if median < minVal:
713 median += cycleLength
714 if median > minVal + cycleLength:
715 median -= cycleLength
716
717 return median
718
719
720
722 """ Get position within tile.
723 """
724 if offPos == 6:
725 return dbc.realDefault(), dbc.realDefault()
726
727 if offPos == 7:
728 return 0, 0
729
730 tileOffSet = 329.779
731 xPos = offPos // 3 or -1
732 yPos = offPos % 3 - 1
733
734 return xPos * tileOffSet, yPos * tileOffSet
735
736
737
739 """ Get position within tile.
740 """
741 if offX < -9.e8 and offY < -9.e8:
742 return 6
743
744 if offX == 0 and offY == 0:
745 return 7
746
747 if offY < 0:
748 jj = 0
749 elif offY == 0:
750 jj = 1
751 else:
752 jj = 2
753
754 return 3 * (0 if offX < 0 else 1) + jj
755
756
757
758 -def matchFrames(frames, raIndex, centrePos, tolerance, extNum=None,
759 extIndex=None):
760 """
761 Returns a list of frames matching an input position. Takes a position, a
762 positional tolerance and a list of frames (or a frame-like list) and return
763 the subset of those matching that position within the specified tolerance.
764
765 Method: identifies the first frame that could be a match based on Dec in
766 the (ordered) input list, then cycles through the list checking all
767 possible matches up to the most distant Dec possible. Note that if the RA
768 and Dec come from a frame that is contained within the given list, that
769 frame will appear as a matched frame in the returned list in addition to
770 any other matches. The optional arguments are provided so that if frame
771 matching required that only common extension numbers are matched, and that
772 information is provided in those input arguments, then the function will
773 only match those.
774
775 @param frames: List of frames to match.
776 @type frames: list(tuple)
777 @param raIndex: Index in frame tuple of the RA co-ordinate (Dec is
778 assumed to be index+1).
779 @type raIndex: int
780 @param centrePos: Co-ordinate of central position (RA, Dec) in degrees.
781 @type centrePos: tuple(float, float)
782 @param tolerance: Positional matching tolerance for frames in a set
783 (degrees).
784 @type tolerance: float
785 @param extNum: Optionally specify a particular image extension to be
786 matched.
787 @type extNum: int
788 @param extIndex: If matching by image extension then provide the index of
789 the frame tuple that contains the extension number.
790 @type extIndex: int
791
792 @return: A subset of frames that match the given position.
793 @rtype: list(tuple)
794
795 """
796
797 if not frames:
798 return []
799
800 ra, dec = centrePos
801 matchCentre = (math.radians(ra), math.radians(dec))
802 matchTol = math.radians(tolerance)
803 matchingFrames = list()
804
805
806
807
808 i = bisect([frame[raIndex + 1] for frame in frames], dec - tolerance)
809 i = max(i - 1, 0)
810
811
812 while frames[i][raIndex + 1] <= dec + tolerance:
813
814 if not extNum or frames[i][extIndex] is extNum:
815 frameCentre = (math.radians(frames[i][raIndex]),
816 math.radians(frames[i][raIndex + 1]))
817 if angularSep(matchCentre, frameCentre) <= matchTol:
818 matchingFrames.append(frames[i])
819
820 i += 1
821 if i == len(frames):
822 break
823
824 return matchingFrames
825
826
827
828 -def _matchPairs(sourcesM, sourcesS, decsS, tolerance, cartSep):
829 """
830 Matches given lists of detections and returns a pointing array.
831
832 @param sourcesM: List of sorted master detection declinations, record
833 number and Cartesian co-ordinates.
834 @type sourcesM: list(tuple(float, (int, float, float, float)))
835 @param sourcesS: List of sorted secondary detection declinations, record
836 number and Cartesian co-ordinates.
837 @type sourcesS: list(tuple(float, (int, float, float, float)))
838 @param decsS: List of sorted secondary detection declinations.
839 @type decsS: list(float)
840 @param tolerance: Pairing criterion in arcseconds with additional 10%.
841 @type tolerance: float
842 @param cartSep: Pairing criterion in Cartesian distance.
843 @type cartSep: float
844
845 @return: Pointing array of detection indices.
846 @rtype: numpy.array(int)
847
848 @note: Using xrange here instead of range gives major performance benefits.
849 @todo: Need to try a fully vectorised calculation in NumPy for performance.
850 @todo: Should be able to combine all three sub for-loops in this function
851 into one - only doing Cartesian separation calculations when in the
852 declination tolerance range. Also try slicing the list, but need to
853 be certain perf. is good, also can possibly do this as a find type
854 func.
855 """
856 pp = numpy.zeros(len(sourcesM), numpy.int32) - 1
857 if not sourcesM or not sourcesS:
858 return pp
859
860 for decM, (idxM, cxM, cyM, czM) in sourcesM:
861
862 idx = bisect(decsS, decM)
863 idx = max(idx - 1, 0)
864
865
866 for i1 in reversed(xrange(1, idx + 1)):
867 if decM - decsS[i1] >= tolerance:
868 break
869 else:
870 i1 = 1
871 i1 -= 1
872 for i2 in xrange(idx, len(decsS) - 1):
873 if decsS[i2] - decM >= tolerance:
874 break
875 else:
876 i2 = len(decsS) - 2
877
878
879
880 seps = [tuple()] * (2 + i2 - i1)
881 for i in xrange(2 + i2 - i1):
882 _decS, (idxS, cxS, cyS, czS) = sourcesS[i1 + i]
883 dx = cxS - cxM
884 dy = cyS - cyM
885 dz = czS - czM
886 seps[i] = (dx * dx + dy * dy + dz * dz, idxS)
887
888
889 minSep, idxS = min(seps)
890 if math.sqrt(minSep) < cartSep:
891 pp[idxM] = idxS
892
893 return pp
894
895
896
898 """
899 Normalize an angle into the range 0 to 2xPI radians.
900
901 @param radians: The angle to normalize (in radians).
902 @type radians: float
903
904 @return: The angle normalized to the range 0 to 2xPI.
905 @rtype: float
906
907 """
908 angle = math.fmod(radians, 2.0 * math.pi)
909 return angle if angle >= 0.0 else angle + 2.0 * math.pi
910
911
912
913 -def pairObjects(inputList, matchRadius, outPathName, oneWay=False):
914 """
915 Pairs object positions from given binary files, saving the resulting
916 pointings in binary files.
917
918 @param inputList: List of detection data sets, one for each passband, as
919 a tuple of (filePath, (raIndex, raSize, raUnit),
920 (decIndex, decSize, decUnit), recordLength)
921 @param inputList: list(tuple)
922 @param matchRadius: Pairing separation criterion in arcseconds.
923 @type matchRadius: float
924 @param outPathName: Path to write the pointing files too. File name must
925 be included if oneWay is True.
926 @type outPathName: str
927 @param oneWay: If True, don't prepare a reverse pointing file.
928 @type oneWay: bool
929
930 @todo: Make completely schema driven for both data types and units,
931 basically the caller gets raIndex/raSize/raUnit from schema file and
932 also passes raType, e.g. 'd' in this case.
933 @todo: Can be made even more Pythonic and faster if we include the merge
934 table in with this function, avoiding the need for index look-ups
935 and reading/writing to binary files... assuming memory req. are OK.
936 @note: NumPy is useful and fast for the pp array, but slow for the
937 Cartesian vectors since they are stored as a list. If the
938 calculations were fully vectorised then this wouldn't be an issue.
939 """
940 convForUnit = {"DEG": math.radians(1), "HRS": math.pi / 12, "RAD": 1}
941
942 decs = []
943 sources = []
944 for filePath, (raIndex, raSize, raUnit), (decIndex, decSize, decUnit), \
945 recordLength in inputList:
946 decList = []
947 srcList = []
948 if filePath:
949 fileSize = os.path.getsize(filePath)
950 if fileSize % recordLength != 0:
951 raise Exception("File size is a non-integer multiple of row "
952 "length %s: %s" % (recordLength, filePath))
953
954 with open(filePath, 'rb') as detFile:
955 raEnd = raIndex + raSize
956 decEnd = decIndex + decSize
957 for recNum in xrange(fileSize // recordLength):
958 row = detFile.read(recordLength)
959
960 ra = struct.unpack('<d', row[raIndex:raEnd])[0]
961 dec = struct.unpack('<d', row[decIndex:decEnd])[0]
962
963
964 ra *= convForUnit[raUnit]
965 dec *= convForUnit[decUnit]
966 cosDec = math.cos(dec)
967
968 decList.append(3600 * math.degrees(dec))
969 srcList.append((recNum,
970 math.cos(ra) * cosDec, math.sin(ra) * cosDec, math.sin(dec)))
971
972 decs.append(sorted(decList))
973 sources.append(sorted(zip(decList, srcList)))
974
975
976 tol = 1.1 * matchRadius
977 matchRadius = math.radians(matchRadius / 3600)
978 d = numpy.array((1 - math.cos(matchRadius), math.sin(matchRadius)))
979 cartSep = math.sqrt(numpy.dot(d, d))
980
981
982 outPath = (outPathName if oneWay else
983 os.path.join(outPathName, "pointers_%s_%s.dat"))
984 for i1 in xrange(len(inputList)):
985 for i2 in xrange(i1 + 1, len(inputList)):
986 pp1 = _matchPairs(sources[i1], sources[i2], decs[i2], tol, cartSep)
987 pp2 = _matchPairs(sources[i2], sources[i1], decs[i1], tol, cartSep)
988
989
990 for index in xrange(len(pp1)):
991 if pp1[index] != -1 and pp2[pp1[index]] != index:
992 pp1[index] = -2
993 pp1.tofile(outPath if oneWay else outPath % (i1, i2))
994
995 if not oneWay:
996 for index in xrange(len(pp2)):
997 if pp2[index] != -1 and pp1[pp2[index]] != index:
998 pp2[index] = -2
999 pp2.tofile(outPath % (i2, i1))
1000
1001
1002
1004 """
1005 Looks up WCS values in the database for given observation.
1006
1007 @param db: Database connection to obtain WCS co-ordinates.
1008 @type db: DbSession
1009 @param multiframeID: MultiframeID of the observation.
1010 @type multiframeID: int
1011 @param extNum: Detector extension number.
1012 @type extNum: int
1013
1014 @return: WCS co-ordinates tuple.
1015 @rtype: tuple(float)
1016
1017 """
1018 if not multiframeID or not extNum:
1019 raise ValueError("Invalid multiframeID (%s) and/or extNum (%s)"
1020 % (multiframeID, extNum))
1021
1022 wcsCols = "crPixX, crPixY, crValX, crValY, cd11, cd12, cd21, cd22"
1023 wcsCols += ", pv21, pv22, pv23" + (", pv24, pv25" if db.sysc.isVSA() else
1024 ", 0., 0.")
1025
1026 return db.query(wcsCols, "CurrentAstrometry",
1027 whereStr="multiframeID=%s AND extNum=%s" % (multiframeID, extNum),
1028 firstOnly=True)
1029
1030
1032 """ Returns correct format of tileRa and tileDec
1033 """
1034 raHMS = raDeg / 15.
1035 raH = int(raHMS)
1036 raHMS = 60.*(raHMS - raH)
1037 raM = int(raHMS)
1038 raHMS = 60.*(raHMS - raM)
1039 raS = round(raHMS, 3)
1040 tileRa = raH * 10000 + raM * 100 + raS
1041 sign = -1 if decDeg < 0 else 1.
1042 decDMS = abs(decDeg)
1043 decD = int(decDMS)
1044 decDMS = 60.*(decDMS - decD)
1045 decM = int(decDMS)
1046 decDMS = 60.*(decDMS - decM)
1047 decS = round(decDMS, 3)
1048 tileDec = sign * (decD * 10000 + decM * 100 + decS)
1049 if outString:
1050 tileRa = str(tileRa)
1051 tileDec = str(tileDec)
1052 return (tileRa, tileDec)
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063