1 """
2 (Mostly deprecated) code to handle coordinate systems and transform
3 between them.
4
5 Basically all of this should be taken over by stc and astropysics.
6 """
7
8
9
10
11
12
13
14 import math
15 import new
16 from math import sin, cos, pi
17 import re
18
19 import numpy
20
21 from gavo import utils
22 from gavo.utils import DEG
23 from gavo.utils import pgsphere
24 from gavo.utils import pyfits
25
26
27 import __builtin__
28 __builtin__._ASTROPY_SETUP_ = True
29 wcs = utils.DeferredImport("wcs", "from astropy import wcs")
30
31
32 fitsKwPat = re.compile("[A-Z0-9_-]{1,8}$")
33
35 """returns a pyfits header with the cards of d.items().
36
37 Only keys "looking like" FITS header keywords are used, i.e. all-uppercase
38 and shorter than 9 characters.
39 """
40 res = pyfits.Header()
41 res.update(dict((str(key), val)
42 for key, val in d.iteritems()
43 if fitsKwPat.match(key) and val is not None))
44 return res
45
46
47 _wcsTestDict = {
48 "CRVAL1": 0, "CRVAL2": 0, "CRPIX1": 50, "CRPIX2": 50,
49 "CD1_1": 0.01, "CD1_2": 0, "CD2_1": 0, "CD2_2": 0.01,
50 "NAXIS1": 100, "NAXIS2": 100, "CUNIT1": "deg", "CUNIT2": "deg",
51 "CTYPE1": 'RA---TAN-SIP', "CTYPE2": 'DEC--TAN-SIP', "LONPOLE": 180.,
52 }
53
54
56 """is a 2D box.
57
58 The can be constructed either with two tuples, giving two points
59 delimiting the box, or with four arguments x0, x1, y0, y1.
60
61 To access the thing, you can either access the x[01], y[01] attributes
62 or use getitem to retrieve the upper right and lower left corner.
63
64 The slightly silly ordering of the bounding points (larger values
65 first) is for consistency with Postgresql.
66 """
67 - def __init__(self, x0, x1, y0=None, y1=None):
68 if y0 is None:
69 x0, y0 = x0
70 x1, y1 = x1
71 lowerLeft = (min(x0, x1), min(y0, y1))
72 upperRight = (max(x0, x1), max(y0, y1))
73 self.x0, self.y0 = upperRight
74 self.x1, self.y1 = lowerLeft
75
77 if index==0 or index==-2:
78 return (self.x0, self.y0)
79 elif index==1 or index==-1:
80 return (self.x1, self.y1)
81 else:
82 raise IndexError("len(box) is always 2")
83
85 return "((%.4g,%.4g), (%.4g,%.4g))"%(self.x0, self.y0, self.x1, self.y1)
86
88 return "Box((%g,%g), (%g,%g))"%(self.x0, self.y0, self.x1, self.y1)
89
90
92 """returns a bounding box for the sequence of 2-sequences points.
93
94 The thing returned is a coords.Box.
95
96 >>> getBbox([(0.25, 1), (-3.75, 1), (-2, 4)])
97 Box((0.25,4), (-3.75,1))
98 """
99 xCoos, yCoos = [[p[i] for p in points] for i in range(2)]
100 return Box(min(xCoos), max(xCoos), min(yCoos), max(yCoos))
101
102
104 while alpha>360:
105 alpha -= 360
106 while alpha<0:
107 alpha += 360
108 return alpha
109
110
112 return max(-90, min(90, delta))
113
114
116 """returns true if something bordered by minRA and maxRA presumably straddles
117 the stitching line.
118
119 This assumes minRA<maxRA, and that "something" is less than 180 degrees in
120 longitude.
121
122 Angles are in degrees here.
123 """
124 return maxRA>270 and minRA<90
125
126
145
146
148 """monkeypatches wcs.WCS instances for DaCHS' purposes.
149 """
150 wcsObj._dachs_header = wcsFields
151 wcsObj.longAxis = naxis[0]
152 if len(naxis)>1:
153 wcsObj.latAxis = naxis[1]
154 wcsObj._monkey_naxis_lengths = [wcsFields.get("NAXIS%d"%i)
155 for i in naxis]
156 wcsObj.calcFootprint = new.instancemethod(_calcFootprintMonkeypatch,
157 wcsObj, wcsObj.__class__)
158
159
160 -def getWCS(wcsFields, naxis=(1,2), relax=True):
161 """returns a WCS instance from wcsFields
162
163 wcsFields can be either a dictionary or a pyfits header giving
164 some kind of WCS information, or an wcs.WCS instance that is
165 returned verbatim.
166
167 This will return None if no (usable) WCS information is found in the header.
168
169 We monkeypatch the resulting WCS structure quite a bit. Among
170 others:
171
172 * there's longAxis and latAxis attributes taken from naxis
173 * there's _dachs_header, containing the incoming k-v pairs
174 * there's _monkey_naxis_length, the lengths along the WCS axes.
175 """
176 if isinstance(wcsFields, wcs.WCS):
177 return wcsFields
178 if isinstance(wcsFields, dict):
179 wcsFields = makePyfitsFromDict(wcsFields)
180
181
182
183 if ("CD1_1" not in wcsFields
184 and "CDELT1" not in wcsFields
185 and "PC1_1" not in wcsFields
186 and "CNPIX1" not in wcsFields):
187 return None
188
189
190 for key in ["AP_ORDER", "BP_ORDER", "A_ORDER", "B_ORDER"]:
191 if wcsFields.get(key)==0:
192 del wcsFields[key]
193
194
195
196
197 for key in ["XPIXELSZ", "YPIXELSZ"]:
198 if key in wcsFields:
199 del wcsFields[key]
200
201 wcsObj = wcs.WCS(wcsFields, relax=relax, naxis=naxis)
202 _monkeypatchWCS(wcsObj, naxis, wcsFields)
203
204 return wcsObj
205
206
208 """returns the focal plane coordindates for the 2-sequence pixels.
209
210 (this is a thin wrapper intended to abstract for pix2world's funky
211 calling convention; also, we fix on the silly "0 pixel is 1 convention")
212 """
213 wcsObj = getWCS(wcsFields)
214 val = wcsObj.pix2foc((pixels[0],), (pixels[1],), 1)
215 return val[0][0], val[1][0]
216
217
219 """returns the sky coordindates for the 2-sequence pixels.
220
221 (this is a thin wrapper intended to abstract for pix2world's funky
222 calling convention; also, we fix on the silly "0 pixel is 1 convention")
223 """
224 wcsObj = getWCS(wcsFields)
225 val = wcsObj.all_pix2world((pixels[0],), (pixels[1],), 1)
226 return val[0][0], val[1][0]
227
228
230 """returns the pixel coordindates for the 2-sequence longLad.
231
232 (this is a thin wrapper intended to abstract for world2pix's funky
233 calling convention; also, we fix on the silly "0 pixel is 1 convention")
234 """
235 val = getWCS(wcsFields).wcs_world2pix((longLat[0],), (longLat[1],), 1)
236 return val[0][0], val[1][0]
237
238
240 """returns the sizes of a pixel at roughly the center of the image for
241 wcsFields.
242
243 Near the pole, this gets a bit weird; we do some limitation of the width
244 of RA pixels there.
245 """
246 wcs = getWCS(wcsFields)
247 width, height = wcs._dachs_header["NAXIS1"], wcs._dachs_header["NAXIS2"]
248 cosDelta = max(0.01, math.cos(pix2sky(wcs, (width/2, height/2))[1]*DEG))
249
250 p0 = pix2sky(wcs, (width/2, height/2))
251 p1 = pix2sky(wcs, (width/2+1, height/2))
252 p2 = pix2sky(wcs, (width/2, height/2+1))
253 return abs(p1[0]-p0[0])*cosDelta, abs(p2[1]-p0[1])
254
255
257 """returns a callable transforming pixel to physical coordinates.
258
259 wcsFields is passed to getWCS, see there for legal types.
260 """
261 wcs = getWCS(wcsFields)
262 return lambda x, y: pix2sky(wcs, (x, y))
263
264
266 """returns a callable transforming physical to pixel coordinates.
267
268 wcsFields is passed to getWCS, see there for legal types.
269 """
270 wcs = getWCS(wcsFields)
271 return lambda ra, dec: sky2pix(wcs, (ra,dec))
272
273
275 """returns a bbox and a field center for WCS FITS header fields.
276
277 wcsFields is passed to getWCS, see there for legal types.
278
279 Warning: this is different from wcs.calcFootprint in that
280 we keep negative angles if the stitching line is crossed; also,
281 no distortions or anything like that are taken into account.
282
283 This code is only used for bboxSIAP, and you must not use it
284 for anything else; it's going to disappear with it.
285 """
286 wcs = getWCS(wcsFields)
287 width, height = float(wcs._dachs_header["NAXIS1"]
288 ), float(wcs._dachs_header["NAXIS2"])
289 cA, cD = pix2sky(wcs, (width/2., height/2.))
290 wA, wD = getPixelSizeDeg(wcs)
291 wA *= width/2.
292 wD *= height/2.
293
294 bounds = [(cA+wA, cD+wD), (cA-wA, cD-wD), (cA+wA, cD-wD),
295 (cA-wA, cD+wD)]
296 bbox = getBbox(bounds)
297 if bbox[0][1]>89:
298 bbox = Box((0, clampDelta(bbox[0][1])), (360, clampDelta(bbox[1][1])))
299 if bbox[1][1]<-89:
300 bbox = Box((0, clampDelta(bbox[0][1])), (360, clampDelta(bbox[1][1])))
301 return bbox
302
303
305 """returns a pgsphere spoly corresponding to wcsFields
306
307 wcsFields is passed to getWCS, see there for legal types.
308
309 The polygon returned is computed by using the four corner points
310 assuming a rectangular image. This typically is only loosely related
311 to a proper spherical polygon describing the shape, as image boundaries
312 in the usual projects are not great circles.
313
314 Also, the spoly will be in the coordinate system of the WCS. If that
315 is not ICRS, you'll probably get something incompatible with most of the
316 VO.
317 """
318 wcs = getWCS(wcsFields)
319 return pgsphere.SPoly([pgsphere.SPoint.fromDegrees(*p)
320 for p in wcs.calcFootprint(wcs._dachs_header)])
321
322
324 """returns RA and Dec of the center of an image described by wcsFields.
325
326 This will use the 1-based axes given by spatialAxes to figure out
327 the pixel lengths of the axes.
328 """
329 wcs = getWCS(wcsFields)
330 center1 = wcs._dachs_header["NAXIS%s"%spatialAxes[0]]/2.
331 center2 = wcs._dachs_header["NAXIS%s"%spatialAxes[1]]/2.
332 return pix2sky(wcs, (center1, center2))
333
334
336 """returns a pgsphere.scircle large enough to cover the image
337 described by wcsFields.
338 """
339 wcs = getWCS(wcsFields)
340 center = getCenterFromWCSFields(wcs)
341
342 height, width = (wcs._dachs_header["NAXIS%s"%spatialAxes[0]],
343 wcs._dachs_header["NAXIS%s"%spatialAxes[1]])
344 radius = max(
345 getGCDist(center, pix2sky(wcs, corner))
346 for corner in [(0, 0), (height, 0), (0, width), (height, width)])
347
348 return pgsphere.SCircle(pgsphere.SPoint.fromDegrees(*center),
349 radius*DEG)
350
351
353 """returns a pair of a wcs.WCS instance and a sequence of
354 the spatial axes.
355
356 This will be None, () if no WCS could be discerned. There's some
357 heuristics involved in the identification of the spatial coordinates
358 that will probably fail for unconventional datasets.
359 """
360 wcsAxes = []
361
362
363 for ind in range(1, hdr["NAXIS"]+1):
364 if "-" in hdr.get("CTYPE%s"%ind, ""):
365 wcsAxes.append(ind)
366
367 if not wcsAxes:
368
369 return None, ()
370
371 if len(wcsAxes)!=2:
372 raise utils.ValidationError("This FITS has !=2"
373 " spatial WCS axes. Please contact the DaCHS authors and"
374 " make them support it.", "PUBDID")
375
376 return getWCS(hdr, naxis=wcsAxes), wcsAxes
377
378
380 """returns pixel cutout slices for covering cooPairs in an image with
381 wcsFields.
382
383 cooPairs is a sequence of (ra, dec) tuples. wcsFields is a DaCHS-enhanced
384 wcs.WCS instance.
385
386 Behaviour if cooPairs use a different coordinate system from wcsFields
387 is undefined at this point.
388
389 Each cutout slice is a tuple of (FITS axis number, lower limit, upper limit).
390
391 If cooPairs is off the wcsFields coverage, a null cutout on the longAxis
392 is returned.
393 """
394 latAxis = wcsFields.latAxis
395 longAxis = wcsFields.longAxis
396 latPixels = wcsFields._dachs_header["NAXIS%d"%latAxis]
397 longPixels = wcsFields._dachs_header["NAXIS%d"%longAxis]
398
399
400
401
402
403 cooPairs = [(
404 min(359.99999, max(-89.9999, ra)),
405 min(89.9999, max(-89.9999, dec)))
406 for ra, dec in cooPairs]
407
408 slices = []
409 pixelFootprint = numpy.asarray(
410 numpy.round(wcsFields.wcs_world2pix(cooPairs, 1)), numpy.int32)
411 pixelLimits = [
412 [min(pixelFootprint[:,0]), max(pixelFootprint[:,0])],
413 [min(pixelFootprint[:,1]), max(pixelFootprint[:,1])]]
414
415
416 if pixelLimits[0][1]<0 or pixelLimits[1][1]<0:
417 return [[longAxis, 0, 0]]
418 if pixelLimits[0][0]>longPixels or pixelLimits[1][0]>latPixels:
419 return [[longAxis, 0, 0]]
420
421
422 pixelLimits = [
423 [max(pixelLimits[0][0], 1), min(pixelLimits[0][1], longPixels)],
424 [max(pixelLimits[1][0], 1), min(pixelLimits[1][1], latPixels)]]
425
426 if pixelLimits[0]!=[1, longPixels]:
427 slices.append([longAxis]+pixelLimits[0])
428 if pixelLimits[1]!=[1, latPixels]:
429 slices.append([latAxis]+pixelLimits[1])
430 return slices
431
432
433
434
436 """is a 3d vector that responds to both .x... and [0]...
437
438 >>> x, y = Vector3(1,2,3), Vector3(2,3,4)
439 >>> x+y
440 Vector3(3.000000,5.000000,7.000000)
441 >>> 4*x
442 Vector3(4.000000,8.000000,12.000000)
443 >>> x*4
444 Vector3(4.000000,8.000000,12.000000)
445 >>> x*y
446 20
447 >>> "%.6f"%abs(x)
448 '3.741657'
449 >>> print(abs((x+y).normalized()))
450 1.0
451 """
453 if isinstance(x, tuple):
454 self.coos = x
455 else:
456 self.coos = (x, y, z)
457
459 return "Vector3(%f,%f,%f)"%tuple(self.coos)
460
462 def cutoff(c):
463 if abs(c)<1e-10:
464 return 0
465 else:
466 return c
467 rounded = [cutoff(c) for c in self.coos]
468 return "[%.2g,%.2g,%.2g]"%tuple(rounded)
469
471 return self.coos[index]
472
474 """does either scalar multiplication if other is not a Vector3, or
475 a scalar product.
476 """
477 if isinstance(other, Vector3):
478 return self.x*other.x+self.y*other.y+self.z*other.z
479 else:
480 return Vector3(self.x*other, self.y*other, self.z*other)
481
482 __rmul__ = __mul__
483
485 return Vector3(self.x/scalar, self.y/scalar, self.z/scalar)
486
488 return Vector3(self.x+other.x, self.y+other.y, self.z+other.z)
489
491 return Vector3(self.x-other.x, self.y-other.y, self.z-other.z)
492
494 return math.sqrt(self.x**2+self.y**2+self.z**2)
495
497 return Vector3(self.y*other.z-self.z*other.y,
498 self.z*other.x-self.x*other.z,
499 self.x*other.y-self.y*other.x)
500
502 return self/abs(self)
503
504 - def getx(self): return self.coos[0]
505 - def setx(self, x): self.coos[0] = x
506 x = property(getx, setx)
507 - def gety(self): return self.coos[1]
508 - def sety(self, y): self.coos[1] = y
509 y = property(gety, sety)
510 - def getz(self): return self.coos[2]
511 - def setz(self, z): self.coos[2] = z
512 z = property(getz, setz)
513
514
516 if a<0:
517 return -1
518 elif a>0:
519 return 1
520 else:
521 return 0
522
523
525
526 """returns the 3d coordinates of the intersection of the direction
527 vector given by the spherical coordinates alpha and delta with the
528 unit sphere.
529
530 alpha and delta are given in degrees.
531
532 >>> print(computeUnitSphereCoords(0,0))
533 [1,0,0]
534 >>> print(computeUnitSphereCoords(0, 90))
535 [0,0,1]
536 >>> print(computeUnitSphereCoords(90, 90))
537 [0,0,1]
538 >>> print(computeUnitSphereCoords(90, 0))
539 [0,1,0]
540 >>> print(computeUnitSphereCoords(180, -45))
541 [-0.71,0,-0.71]
542 """
543 return Vector3(*utils.spherToCart(alpha*DEG, delta*DEG))
544
545
547 """returns alpha, delta in degrees for the direction vector dirVec.
548
549 >>> dirVecToCelCoos(computeUnitSphereCoords(25.25, 12.125))
550 (25.25, 12.125)
551 >>> dirVecToCelCoos(computeUnitSphereCoords(25.25, 12.125)*16)
552 (25.25, 12.125)
553 >>> "%g,%g"%dirVecToCelCoos(computeUnitSphereCoords(25.25, 12.125)+
554 ... computeUnitSphereCoords(30.75, 20.0))
555 '27.9455,16.0801'
556 """
557 dirVec = dirVec.normalized()
558 alpha = math.atan2(dirVec.y, dirVec.x)
559 if alpha<0:
560 alpha += 2*math.pi
561 return alpha*180./math.pi, math.asin(dirVec.z)*180./math.pi
562
563
565 """returns the unit vectors for RA and Dec at the unit circle position cPos.
566
567 We compute them by solving u_1*p_1+u_2*p_2=0 (we already know that
568 u_3=0) simultaneously with u_1^2+u_2^2=1 for RA, and by computing the
569 cross product of the RA unit and the radius vector for dec.
570
571 This becomes degenerate at the poles. If we're exactly on a pole,
572 we *define* the unit vectors as (1,0,0) and (0,1,0).
573
574 Orientation is a pain -- the convention used here is that unit delta
575 always points to the pole.
576
577 >>> cPos = computeUnitSphereCoords(45, -45)
578 >>> ua, ud = getTangentialUnits(cPos)
579 >>> print abs(ua), abs(ud), cPos*ua, cPos*ud
580 1.0 1.0 0.0 0.0
581 >>> print ua, ud
582 [-0.71,0.71,0] [-0.5,-0.5,-0.71]
583 >>> ua, ud = getTangentialUnits(computeUnitSphereCoords(180, 60))
584 >>> print ua, ud
585 [0,-1,0] [0.87,0,0.5]
586 >>> ua, ud = getTangentialUnits(computeUnitSphereCoords(0, 60))
587 >>> print ua, ud
588 [0,1,0] [-0.87,0,0.5]
589 >>> ua, ud = getTangentialUnits(computeUnitSphereCoords(0, -60))
590 >>> print ua, ud
591 [0,1,0] [-0.87,0,-0.5]
592 """
593 try:
594 normalizer = 1/math.sqrt(cPos.x**2+cPos.y**2)
595 except ZeroDivisionError:
596 return Vector3(1,0,0), Vector3(0,1,0)
597 alphaUnit = normalizer*Vector3(cPos.y, -cPos.x, 0)
598 deltaUnit = normalizer*Vector3(cPos.x*cPos.z, cPos.y*cPos.z,
599 -cPos.x**2-cPos.y**2)
600
601 if sgn(cPos.z)!=sgn(deltaUnit.z):
602 deltaUnit = -1*deltaUnit
603
604 if cPos.z<0:
605 if deltaUnit.cross(alphaUnit)*cPos<0:
606 alphaUnit = -1*alphaUnit
607 else:
608 if deltaUnit.cross(alphaUnit)*cPos>0:
609 alphaUnit = -1*alphaUnit
610 return alphaUnit, deltaUnit
611
612
613 -def movePm(alphaDeg, deltaDeg, pmAlpha, pmDelta, timeDiff, foreshort=0):
614 """returns alpha and delta for an object with pm pmAlpha after timeDiff.
615
616 pmAlpha has to have cos(delta) applied, everything is supposed to be
617 in degrees, the time unit is yours to choose.
618 """
619 alpha, delta = alphaDeg/180.*math.pi, deltaDeg/180.*math.pi
620 pmAlpha, pmDelta = pmAlpha/180.*math.pi, pmDelta/180.*math.pi
621 sd, cd = math.sin(delta), math.cos(delta)
622 sa, ca = math.sin(alpha), math.cos(alpha)
623 muAbs = math.sqrt(pmAlpha**2+pmDelta**2);
624 muTot = muAbs+0.5*foreshort*timeDiff;
625
626 if muAbs<1e-20:
627 return alphaDeg, deltaDeg
628
629 dirA = pmAlpha/muAbs;
630 dirD = pmDelta/muAbs;
631 sinMot = sin(muTot*timeDiff);
632 cosMot = cos(muTot*timeDiff);
633
634 dirVec = Vector3(-sd*ca*dirD*sinMot - sa*dirA*sinMot + cd*ca*cosMot,
635 -sd*sa*dirD*sinMot + ca*dirA*sinMot + cd*sa*cosMot,
636 +cd*dirD*sinMot + sd*cosMot)
637 return dirVecToCelCoos(dirVec)
638
639
641 """returns the distance along a great circle between two points.
642
643 The distance is in degrees, the input positions are in degrees.
644 """
645 scalarprod = computeUnitSphereCoords(*pos1)*computeUnitSphereCoords(*pos2)
646
647 if scalarprod>=1:
648 return 0
649 return math.acos(scalarprod)/DEG
650
651
655
656
657 if __name__=="__main__":
658 _test()
659