1 """module for simple .fits image tasks (rotation, clipping out sections, making .pngs etc.)
2
3 (c) 2007-2014 Matt Hilton
4
5 U{http://astlib.sourceforge.net}
6
7 Some routines in this module will fail if, e.g., asked to clip a section from a
8 .fits image at a position not found within the image (as determined using the
9 WCS). Where this occurs, the function will return None. An error message will
10 be printed to the console when this happens if astImages.REPORT_ERRORS=True
11 (the default). Testing if an astImages function returns None can be used to
12 handle errors in scripts.
13
14 """
15
16 import os
17
18 import math
19 from astLib import astWCS
20 import numpy
21
22 REPORT_ERRORS = True
23
24
25
26 try:
27 from astropy.io import fits as pyfits
28 except:
29 try:
30 import pyfits
31 except:
32 raise Exception("couldn't import either pyfits or astropy.io.fits")
33
34 try:
35 from scipy import ndimage
36 from scipy import interpolate
37 except ImportError:
38 print("WARNING: astImages: failed to import scipy.ndimage - some "
39 "functions will not work.")
40 try:
41 import matplotlib
42 from matplotlib import pylab
43 matplotlib.interactive(False)
44 except ImportError:
45 print("WARNING: astImages: failed to import matplotlib - some functions "
46 "will not work.")
47
48
49
50 -def clipImageSectionWCS(imageData,
51 imageWCS,
52 RADeg,
53 decDeg,
54 clipSizeDeg,
55 returnWCS=True):
56 """Clips a square or rectangular section from an image array at the given
57 celestial coordinates. An updated WCS for the clipped section is optionally
58 returned, as well as the x, y pixel coordinates in the original image
59 corresponding to the clipped section.
60
61 Note that the clip size is specified in degrees on the sky. For projections
62 that have varying real pixel scale across the map (e.g. CEA), use
63 L{clipUsingRADecCoords} instead.
64
65 @type imageData: numpy array
66 @param imageData: image data array
67 @type imageWCS: astWCS.WCS
68 @param imageWCS: astWCS.WCS object
69 @type RADeg: float
70 @param RADeg: coordinate in decimal degrees
71 @type decDeg: float
72 @param decDeg: coordinate in decimal degrees
73 @type clipSizeDeg: float or list in format [widthDeg, heightDeg]
74 @param clipSizeDeg: if float, size of square clipped section in decimal
75 degrees; if list,
76 size of clipped section in degrees in x, y axes of image respectively
77 @type returnWCS: bool
78 @param returnWCS: if True, return an updated WCS for the clipped section
79 @rtype: dictionary
80 @return: clipped image section (numpy array), updated astWCS WCS object for
81 clipped image section, and coordinates of clipped section in imageData in
82 format {'data', 'wcs', 'clippedSection'}.
83
84 """
85
86 imHeight = imageData.shape[0]
87 imWidth = imageData.shape[1]
88 xImScale = imageWCS.getXPixelSizeDeg()
89 yImScale = imageWCS.getYPixelSizeDeg()
90
91 if type(clipSizeDeg) == float:
92 xHalfClipSizeDeg = clipSizeDeg / 2.0
93 yHalfClipSizeDeg = xHalfClipSizeDeg
94 elif type(clipSizeDeg) == list or type(clipSizeDeg) == tuple:
95 xHalfClipSizeDeg = clipSizeDeg[0] / 2.0
96 yHalfClipSizeDeg = clipSizeDeg[1] / 2.0
97 else:
98 raise Exception("did not understand clipSizeDeg: should be float, or "
99 "[widthDeg, heightDeg]")
100
101 xHalfSizePix = xHalfClipSizeDeg / xImScale
102 yHalfSizePix = yHalfClipSizeDeg / yImScale
103
104 cPixCoords = imageWCS.wcs2pix(RADeg, decDeg)
105
106 cTopLeft = [cPixCoords[0] + xHalfSizePix, cPixCoords[1] + yHalfSizePix]
107 cBottomRight = [cPixCoords[0] - xHalfSizePix, cPixCoords[1] - yHalfSizePix]
108
109 X = [int(round(cTopLeft[0])), int(round(cBottomRight[0]))]
110 Y = [int(round(cTopLeft[1])), int(round(cBottomRight[1]))]
111
112 X.sort()
113 Y.sort()
114
115 if X[0] < 0:
116 X[0] = 0
117 if X[1] > imWidth:
118 X[1] = imWidth
119 if Y[0] < 0:
120 Y[0] = 0
121 if Y[1] > imHeight:
122 Y[1] = imHeight
123
124 clippedData = imageData[Y[0]:Y[1], X[0]:X[1]]
125
126
127 if returnWCS:
128 try:
129 oldCRPIX1 = imageWCS.header['CRPIX1']
130 oldCRPIX2 = imageWCS.header['CRPIX2']
131 clippedWCS = imageWCS.copy()
132 clippedWCS.header['NAXIS1'] = clippedData.shape[1]
133 clippedWCS.header['NAXIS2'] = clippedData.shape[0]
134 clippedWCS.header['CRPIX1'] = oldCRPIX1 - X[0]
135 clippedWCS.header['CRPIX2'] = oldCRPIX2 - Y[0]
136 clippedWCS.updateFromHeader()
137
138 except KeyError:
139
140 if REPORT_ERRORS:
141
142 print("WARNING: astImages.clipImageSectionWCS() : no CRPIX1, "
143 "CRPIX2 keywords found - not updating clipped image "
144 "WCS.")
145
146 clippedData = imageData[Y[0]:Y[1], X[0]:X[1]]
147 clippedWCS = imageWCS.copy()
148 else:
149 clippedWCS = None
150
151 return {'data': clippedData,
152 'wcs': clippedWCS,
153 'clippedSection': [X[0], X[1], Y[0], Y[1]]}
154
155
156
158 """Clips a square or rectangular section from an image array at the given
159 pixel coordinates.
160
161 @type imageData: numpy array
162 @param imageData: image data array
163 @type XCoord: float
164 @param XCoord: coordinate in pixels
165 @type YCoord: float
166 @param YCoord: coordinate in pixels
167 @type clipSizePix: float or list in format [widthPix, heightPix]
168 @param clipSizePix: if float, size of square clipped section in pixels; if
169 list, size of clipped section in pixels in x, y axes of output image
170 respectively
171 @rtype: numpy array
172 @return: clipped image section
173
174 """
175
176 imHeight = imageData.shape[0]
177 imWidth = imageData.shape[1]
178
179 if type(clipSizePix) == float or type(clipSizePix) == int:
180 xHalfClipSizePix = int(round(clipSizePix / 2.0))
181 yHalfClipSizePix = xHalfClipSizePix
182 elif type(clipSizePix) == list or type(clipSizePix) == tuple:
183 xHalfClipSizePix = int(round(clipSizePix[0] / 2.0))
184 yHalfClipSizePix = int(round(clipSizePix[1] / 2.0))
185 else:
186 raise Exception("did not understand clipSizePix: should be float, or "
187 "[widthPix, heightPix]")
188
189 cTopLeft = [XCoord + xHalfClipSizePix, YCoord + yHalfClipSizePix]
190 cBottomRight = [XCoord - xHalfClipSizePix, YCoord - yHalfClipSizePix]
191
192 X = [int(round(cTopLeft[0])), int(round(cBottomRight[0]))]
193 Y = [int(round(cTopLeft[1])), int(round(cBottomRight[1]))]
194
195 X.sort()
196 Y.sort()
197
198 if X[0] < 0:
199 X[0] = 0
200 if X[1] > imWidth:
201 X[1] = imWidth
202 if Y[0] < 0:
203 Y[0] = 0
204 if Y[1] > imHeight:
205 Y[1] = imHeight
206
207 return imageData[Y[0]:Y[1], X[0]:X[1]]
208
209
216 """Clips a square or rectangular section from an image array at the given
217 celestial coordinates. The resulting clip is rotated and/or flipped such
218 that North is at the top, and East appears at the left. An updated WCS for
219 the clipped section is also returned. Note that the alignment of the
220 rotated WCS is currently not perfect - however, it is probably good enough
221 in most cases for use with L{ImagePlot} for plotting purposes.
222
223 Note that the clip size is specified in degrees on the sky. For projections
224 that have varying real pixel scale across the map (e.g. CEA), use
225 L{clipUsingRADecCoords} instead.
226
227 @type imageData: numpy array
228 @param imageData: image data array
229 @type imageWCS: astWCS.WCS
230 @param imageWCS: astWCS.WCS object
231 @type RADeg: float
232 @param RADeg: coordinate in decimal degrees
233 @type decDeg: float
234 @param decDeg: coordinate in decimal degrees
235 @type clipSizeDeg: float
236 @param clipSizeDeg: if float, size of square clipped section in decimal
237 degrees; if list, size of clipped section in degrees in RA, dec. axes of
238 output rotated image respectively
239 @type returnWCS: bool
240 @param returnWCS: if True, return an updated WCS for the clipped section
241 @rtype: dictionary
242 @return: clipped image section (numpy array), updated astWCS WCS object for
243 clipped image section, in format {'data', 'wcs'}.
244
245 @note: Returns 'None' if the requested position is not found within the
246 image. If the imaged WCS does not have keywords of the form CD1_1 etc., the
247 output WCS will not be rotated.
248
249 """
250
251 halfImageSize = imageWCS.getHalfSizeDeg()
252 imageCentre = imageWCS.getCentreWCSCoords()
253
254
255 if type(clipSizeDeg) == float:
256 xHalfClipSizeDeg = clipSizeDeg / 2.0
257 yHalfClipSizeDeg = xHalfClipSizeDeg
258 elif type(clipSizeDeg) == list or type(clipSizeDeg) == tuple:
259 xHalfClipSizeDeg = clipSizeDeg[0] / 2.0
260 yHalfClipSizeDeg = clipSizeDeg[1] / 2.0
261 else:
262 raise Exception("did not understand clipSizeDeg: should be float, or "
263 "[widthDeg, heightDeg]")
264
265 diagonalHalfSizeDeg = math.sqrt((xHalfClipSizeDeg * xHalfClipSizeDeg) +
266 (yHalfClipSizeDeg * yHalfClipSizeDeg))
267
268
269
270 if RADeg > imageCentre[0] - halfImageSize[0] and RADeg < imageCentre[0] + \
271 halfImageSize[0] and decDeg > imageCentre[1] - halfImageSize[1] and \
272 decDeg < imageCentre[1] + halfImageSize[1]:
273
274 imageDiagonalClip = clipImageSectionWCS(
275 imageData, imageWCS, RADeg, decDeg, diagonalHalfSizeDeg * 2.0)
276 diagonalClip = imageDiagonalClip['data']
277 diagonalWCS = imageDiagonalClip['wcs']
278
279 rotDeg = diagonalWCS.getRotationDeg()
280 imageRotated = ndimage.rotate(diagonalClip, rotDeg)
281 if diagonalWCS.isFlipped() == 1:
282 imageRotated = pylab.fliplr(imageRotated)
283
284
285 rotatedWCS = diagonalWCS.copy()
286 rotRadians = math.radians(rotDeg)
287
288 if returnWCS:
289 try:
290
291 CD11 = rotatedWCS.header['CD1_1']
292 CD21 = rotatedWCS.header['CD2_1']
293 CD12 = rotatedWCS.header['CD1_2']
294 CD22 = rotatedWCS.header['CD2_2']
295 if rotatedWCS.isFlipped() == 1:
296 CD11 = CD11 * -1
297 CD12 = CD12 * -1
298 CDMatrix = numpy.array([[CD11, CD12], [CD21, CD22]],
299 dtype=numpy.float64)
300
301 rotRadians = rotRadians
302 rot11 = math.cos(rotRadians)
303 rot12 = math.sin(rotRadians)
304 rot21 = -math.sin(rotRadians)
305 rot22 = math.cos(rotRadians)
306 rotMatrix = numpy.array([[rot11, rot12], [rot21, rot22]],
307 dtype=numpy.float64)
308 newCDMatrix = numpy.dot(rotMatrix, CDMatrix)
309
310 P1 = diagonalWCS.header['CRPIX1']
311 P2 = diagonalWCS.header['CRPIX2']
312 V1 = diagonalWCS.header['CRVAL1']
313 V2 = diagonalWCS.header['CRVAL2']
314
315 PMatrix = numpy.zeros((2, ), dtype=numpy.float64)
316 PMatrix[0] = P1
317 PMatrix[1] = P2
318
319
320 CMatrix = numpy.array(
321 [imageRotated.shape[1] / 2.0, imageRotated.shape[0] / 2.0])
322 centreCoords = diagonalWCS.getCentreWCSCoords()
323 alphaRad = math.radians(centreCoords[0])
324 deltaRad = math.radians(centreCoords[1])
325 thetaRad = math.asin(math.sin(deltaRad) *
326 math.sin(math.radians(V2)) +
327 math.cos(deltaRad) *
328 math.cos(math.radians(V2)) *
329 math.cos(alphaRad - math.radians(V1)))
330 phiRad = math.atan2(-math.cos(deltaRad) *
331 math.sin(alphaRad - math.radians(V1)),
332 math.sin(deltaRad) *
333 math.cos(math.radians(V2)) -
334 math.cos(deltaRad) *
335 math.sin(math.radians(V2)) *
336 math.cos(alphaRad - math.radians(V1))) + \
337 math.pi
338 RTheta = (180.0 / math.pi) * (1.0 / math.tan(thetaRad))
339
340 xy = numpy.zeros((2, ), dtype=numpy.float64)
341 xy[0] = RTheta * math.sin(phiRad)
342 xy[1] = -RTheta * math.cos(phiRad)
343 newPMatrix = CMatrix - numpy.dot(
344 numpy.linalg.inv(newCDMatrix), xy)
345
346
347
348
349
350
351
352
353 rotatedWCS.header['NAXIS1'] = imageRotated.shape[1]
354 rotatedWCS.header['NAXIS2'] = imageRotated.shape[0]
355 rotatedWCS.header['CRPIX1'] = newPMatrix[0]
356 rotatedWCS.header['CRPIX2'] = newPMatrix[1]
357 rotatedWCS.header['CRVAL1'] = V1
358 rotatedWCS.header['CRVAL2'] = V2
359 rotatedWCS.header['CD1_1'] = newCDMatrix[0][0]
360 rotatedWCS.header['CD2_1'] = newCDMatrix[1][0]
361 rotatedWCS.header['CD1_2'] = newCDMatrix[0][1]
362 rotatedWCS.header['CD2_2'] = newCDMatrix[1][1]
363 rotatedWCS.updateFromHeader()
364
365 except KeyError:
366
367 if REPORT_ERRORS:
368 print("WARNING: astImages.clipRotatedImageSectionWCS() : "
369 "no CDi_j keywords found - not rotating WCS.")
370
371 imageRotated = diagonalClip
372 rotatedWCS = diagonalWCS
373
374 imageRotatedClip = clipImageSectionWCS(imageRotated, rotatedWCS, RADeg,
375 decDeg, clipSizeDeg)
376
377 if returnWCS:
378 return {'data': imageRotatedClip['data'],
379 'wcs': imageRotatedClip['wcs']}
380 else:
381 return {'data': imageRotatedClip['data'], 'wcs': None}
382
383 else:
384
385 if REPORT_ERRORS:
386 print("""ERROR: astImages.clipRotatedImageSectionWCS() :
387 RADeg, decDeg are not within imageData.""")
388
389 return None
390
391
392
393 -def clipUsingRADecCoords(imageData,
394 imageWCS,
395 RAMin,
396 RAMax,
397 decMin,
398 decMax,
399 returnWCS=True):
400 """Clips a section from an image array at the pixel coordinates
401 corresponding to the given celestial coordinates.
402
403 @type imageData: numpy array
404 @param imageData: image data array
405 @type imageWCS: astWCS.WCS
406 @param imageWCS: astWCS.WCS object
407 @type RAMin: float
408 @param RAMin: minimum RA coordinate in decimal degrees
409 @type RAMax: float
410 @param RAMax: maximum RA coordinate in decimal degrees
411 @type decMin: float
412 @param decMin: minimum dec coordinate in decimal degrees
413 @type decMax: float
414 @param decMax: maximum dec coordinate in decimal degrees
415 @type returnWCS: bool
416 @param returnWCS: if True, return an updated WCS for the clipped section
417 @rtype: dictionary
418 @return: clipped image section (numpy array), updated astWCS WCS object for
419 clipped image section, and corresponding pixel coordinates in imageData in
420 format {'data', 'wcs', 'clippedSection'}.
421
422 @note: Returns 'None' if the requested position is not found within the
423 image.
424
425 """
426
427 imHeight = imageData.shape[0]
428 imWidth = imageData.shape[1]
429
430 xMin, yMin = imageWCS.wcs2pix(RAMin, decMin)
431 xMax, yMax = imageWCS.wcs2pix(RAMax, decMax)
432 xMin = int(round(xMin))
433 xMax = int(round(xMax))
434 yMin = int(round(yMin))
435 yMax = int(round(yMax))
436 X = [xMin, xMax]
437 X.sort()
438 Y = [yMin, yMax]
439 Y.sort()
440
441 if X[0] < 0:
442 X[0] = 0
443 if X[1] > imWidth:
444 X[1] = imWidth
445 if Y[0] < 0:
446 Y[0] = 0
447 if Y[1] > imHeight:
448 Y[1] = imHeight
449
450 clippedData = imageData[Y[0]:Y[1], X[0]:X[1]]
451
452
453 if returnWCS:
454 try:
455 oldCRPIX1 = imageWCS.header['CRPIX1']
456 oldCRPIX2 = imageWCS.header['CRPIX2']
457 clippedWCS = imageWCS.copy()
458 clippedWCS.header['NAXIS1'] = clippedData.shape[1]
459 clippedWCS.header['NAXIS2'] = clippedData.shape[0]
460 clippedWCS.header['CRPIX1'] = oldCRPIX1 - X[0]
461 clippedWCS.header['CRPIX2'] = oldCRPIX2 - Y[0]
462 clippedWCS.updateFromHeader()
463
464 except KeyError:
465
466 if REPORT_ERRORS:
467
468 print("WARNING: astImages.clipUsingRADecCoords() : no CRPIX1, "
469 "CRPIX2 keywords found - not updating clipped image"
470 "WCS.")
471
472 clippedData = imageData[Y[0]:Y[1], X[0]:X[1]]
473 clippedWCS = imageWCS.copy()
474 else:
475 clippedWCS = None
476
477 return {'data': clippedData,
478 'wcs': clippedWCS,
479 'clippedSection': [X[0], X[1], Y[0], Y[1]]}
480
481
482
484 """Scales image array and WCS by the given scale factor.
485
486 @type imageData: numpy array
487 @param imageData: image data array
488 @type imageWCS: astWCS.WCS
489 @param imageWCS: astWCS.WCS object
490 @type scaleFactor: float or list or tuple
491 @param scaleFactor: factor to resize image by - if tuple or list, in format
492 [x scale factor, y scale factor]
493 @rtype: dictionary
494 @return: image data (numpy array), updated astWCS WCS object for image, in
495 format {'data', 'wcs'}.
496
497 """
498
499 if type(scaleFactor) == int or type(scaleFactor) == float:
500 scaleFactor = [float(scaleFactor), float(scaleFactor)]
501 scaledData = ndimage.zoom(imageData, scaleFactor)
502
503
504
505 properDimensions = numpy.array(imageData.shape) * scaleFactor
506 offset = properDimensions - numpy.array(scaledData.shape)
507
508
509 try:
510 oldCRPIX1 = imageWCS.header['CRPIX1']
511 oldCRPIX2 = imageWCS.header['CRPIX2']
512 CD11 = imageWCS.header['CD1_1']
513 CD21 = imageWCS.header['CD2_1']
514 CD12 = imageWCS.header['CD1_2']
515 CD22 = imageWCS.header['CD2_2']
516 except KeyError:
517
518 try:
519 oldCRPIX1 = imageWCS.header['CRPIX1']
520 oldCRPIX2 = imageWCS.header['CRPIX2']
521 CD11 = imageWCS.header['CDELT1']
522 CD21 = 0
523 CD12 = 0
524 CD22 = imageWCS.header['CDELT2']
525 except KeyError:
526 if REPORT_ERRORS:
527 print("WARNING: astImages.rescaleImage() : no CDij or CDELT "
528 "keywords found - not updating WCS.")
529 scaledWCS = imageWCS.copy()
530 return {'data': scaledData, 'wcs': scaledWCS}
531
532 CDMatrix = numpy.array([[CD11, CD12], [CD21, CD22]], dtype=numpy.float64)
533 scaleFactorMatrix = numpy.array(
534 [[1.0 / scaleFactor[0], 0], [0, 1.0 / scaleFactor[1]]])
535 scaledCDMatrix = numpy.dot(scaleFactorMatrix, CDMatrix)
536
537 scaledWCS = imageWCS.copy()
538 scaledWCS.header['NAXIS1'] = scaledData.shape[1]
539 scaledWCS.header['NAXIS2'] = scaledData.shape[0]
540 scaledWCS.header['CRPIX1'] = oldCRPIX1 * scaleFactor[0] + offset[1]
541 scaledWCS.header['CRPIX2'] = oldCRPIX2 * scaleFactor[1] + offset[0]
542 scaledWCS.header['CD1_1'] = scaledCDMatrix[0][0]
543 scaledWCS.header['CD2_1'] = scaledCDMatrix[1][0]
544 scaledWCS.header['CD1_2'] = scaledCDMatrix[0][1]
545 scaledWCS.header['CD2_2'] = scaledCDMatrix[1][1]
546 scaledWCS.updateFromHeader()
547
548 return {'data': scaledData, 'wcs': scaledWCS}
549
550
552 """Creates a matplotlib.pylab plot of an image array with the specified
553 cuts in intensity applied. This routine is used by L{saveBitmap} and
554 L{saveContourOverlayBitmap}, which both produce output as .png, .jpg, etc.
555 images.
556
557 @type imageData: numpy array
558 @param imageData: image data array
559 @type cutLevels: list
560 @param cutLevels: sets the image scaling - available options:
561 - pixel values: cutLevels=[low value, high value].
562 - histogram equalisation: cutLevels=["histEq", number of bins ( e.g.
563 1024)]
564 - relative: cutLevels=["relative", cut per cent level (e.g. 99.5)]
565 - smart: cutLevels=["smart", cut per cent level (e.g. 99.5)]
566 ["smart", 99.5] seems to provide good scaling over a range of different
567 images.
568 @rtype: dictionary
569 @return: image section (numpy.array), matplotlib image normalisation
570 (matplotlib.colors.Normalize), in the format {'image', 'norm'}.
571
572 @note: If cutLevels[0] == "histEq", then only {'image'} is returned.
573
574 """
575
576
577 if cutLevels[0] == "histEq":
578
579 imageData = histEq(imageData, cutLevels[1])
580 anorm = pylab.normalize(imageData.min(), imageData.max())
581
582 elif cutLevels[0] == "relative":
583
584
585 sorted = numpy.sort(numpy.ravel(imageData))
586 maxValue = sorted.max()
587 minValue = sorted.min()
588
589
590 topCutIndex = len(sorted - 1) - \
591 int(math.floor(float((100.0 - cutLevels[1]) / 100.0) *
592 len(sorted - 1)))
593 bottomCutIndex = int(math.ceil(float((100.0 - cutLevels[1]) / 100.0) *
594 len(sorted - 1)))
595 topCut = sorted[topCutIndex]
596 bottomCut = sorted[bottomCutIndex]
597 anorm = pylab.normalize(bottomCut, topCut)
598
599 elif cutLevels[0] == "smart":
600
601
602 sorted = numpy.sort(numpy.ravel(imageData))
603 maxValue = sorted.max()
604 minValue = sorted.min()
605 numBins = 10000
606 binWidth = (maxValue - minValue) / float(numBins)
607 histogram = ndimage.histogram(sorted, minValue, maxValue, numBins)
608
609
610
611
612
613
614
615
616 backgroundValue = histogram.max()
617 foundBackgroundBin = False
618 foundTopBin = False
619 lastBin = -10000
620 for i in range(len(histogram)):
621
622 if histogram[i] >= lastBin and foundBackgroundBin:
623
624
625
626 if (minValue + (binWidth * i)) > bottomBinValue * 1.1:
627 topBinValue = minValue + (binWidth * i)
628 foundTopBin = True
629 break
630
631 if histogram[i] == backgroundValue and not foundBackgroundBin:
632 bottomBinValue = minValue + (binWidth * i)
633 foundBackgroundBin = True
634
635 lastBin = histogram[i]
636
637 if not foundTopBin:
638 topBinValue = maxValue
639
640
641 smartClipped = numpy.clip(sorted, bottomBinValue, topBinValue)
642 topCutIndex = len(smartClipped - 1) - \
643 int(math.floor(float((100.0 - cutLevels[1]) / 100.0) *
644 len(smartClipped - 1)))
645 bottomCutIndex = int(math.ceil(float((100.0 - cutLevels[1]) / 100.0) *
646 len(smartClipped - 1)))
647 topCut = smartClipped[topCutIndex]
648 bottomCut = smartClipped[bottomCutIndex]
649 anorm = pylab.normalize(bottomCut, topCut)
650 else:
651
652
653 anorm = pylab.normalize(cutLevels[0], cutLevels[1])
654
655 if cutLevels[0] == "histEq":
656 return {'image': imageData.copy()}
657 else:
658 return {'image': imageData.copy(), 'norm': anorm}
659
660
661
665 """Resamples an image and WCS to a tangent plane projection. Purely for
666 plotting purposes (e.g., ensuring RA, dec. coordinate axes perpendicular).
667
668 @type imageData: numpy array
669 @param imageData: image data array
670 @type imageWCS: astWCS.WCS
671 @param imageWCS: astWCS.WCS object
672 @type outputPixDimensions: list
673 @param outputPixDimensions: [width, height] of output image in pixels
674 @rtype: dictionary
675 @return: image data (numpy array), updated astWCS WCS object for image, in
676 format {'data', 'wcs'}.
677
678 """
679
680 RADeg, decDeg = imageWCS.getCentreWCSCoords()
681
682
683 xSizeDeg, ySizeDeg = imageWCS.getFullSizeSkyDeg()
684 xSizePix = int(round(outputPixDimensions[0]))
685 ySizePix = int(round(outputPixDimensions[1]))
686 xRefPix = xSizePix / 2.0
687 yRefPix = ySizePix / 2.0
688 xOutPixScale = xSizeDeg / xSizePix
689
690 newHead = pyfits.Header()
691 newHead['NAXIS'] = 2
692 newHead['NAXIS1'] = xSizePix
693 newHead['NAXIS2'] = ySizePix
694 newHead['CTYPE1'] = 'RA---TAN'
695 newHead['CTYPE2'] = 'DEC--TAN'
696 newHead['CRVAL1'] = RADeg
697 newHead['CRVAL2'] = decDeg
698 newHead['CRPIX1'] = xRefPix + 1
699 newHead['CRPIX2'] = yRefPix + 1
700 newHead['CDELT1'] = -xOutPixScale
701 newHead['CDELT2'] = xOutPixScale
702 newHead['CUNIT1'] = 'DEG'
703 newHead['CUNIT2'] = 'DEG'
704 newWCS = astWCS.WCS(newHead, mode='pyfits')
705 newImage = numpy.zeros([ySizePix, xSizePix])
706
707 tanImage = resampleToWCS(newImage,
708 newWCS,
709 imageData,
710 imageWCS,
711 highAccuracy=True,
712 onlyOverlapping=False)
713
714 return tanImage
715
716
717
718 -def resampleToWCS(im1Data,
719 im1WCS,
720 im2Data,
721 im2WCS,
722 highAccuracy=False,
723 onlyOverlapping=True):
724 """Resamples data corresponding to second image (with data im2Data, WCS
725 im2WCS) onto the WCS of the first image (im1Data, im1WCS). The output,
726 resampled image is of the pixel same dimensions of the first image. This
727 routine is for assisting in plotting - performing photometry on the output
728 is not recommended.
729
730 Set highAccuracy == True to sample every corresponding pixel in each image;
731 otherwise only every nth pixel (where n is the ratio of the image scales)
732 will be sampled, with values in between being set using a linear
733 interpolation (much faster).
734
735 Set onlyOverlapping == True to speed up resampling by only resampling the
736 overlapping area defined by both image WCSs.
737
738 @type im1Data: numpy array
739 @param im1Data: image data array for first image
740 @type im1WCS: astWCS.WCS
741 @param im1WCS: astWCS.WCS object corresponding to im1Data
742 @type im2Data: numpy array
743 @param im2Data: image data array for second image (to be resampled to match
744 first image)
745 @type im2WCS: astWCS.WCS
746 @param im2WCS: astWCS.WCS object corresponding to im2Data
747 @type highAccuracy: bool
748 @param highAccuracy: if True, sample every corresponding pixel in each
749 image; otherwise, sample
750 every nth pixel, where n = the ratio of the image scales.
751 @type onlyOverlapping: bool
752 @param onlyOverlapping: if True, only consider the overlapping area defined
753 by both image WCSs (speeds things up)
754 @rtype: dictionary
755 @return: numpy image data array and associated WCS in format {'data', 'wcs'}
756
757 """
758
759 resampledData = numpy.zeros(im1Data.shape)
760
761
762
763
764
765 xPixRatio = (im2WCS.getXPixelSizeDeg() / im1WCS.getXPixelSizeDeg()) / 2.0
766 yPixRatio = (im2WCS.getYPixelSizeDeg() / im1WCS.getYPixelSizeDeg()) / 2.0
767 xBorder = xPixRatio * 10.0
768 yBorder = yPixRatio * 10.0
769 if not highAccuracy:
770 if xPixRatio > 1:
771 xPixStep = int(math.ceil(xPixRatio))
772 else:
773 xPixStep = 1
774 if yPixRatio > 1:
775 yPixStep = int(math.ceil(yPixRatio))
776 else:
777 yPixStep = 1
778 else:
779 xPixStep = 1
780 yPixStep = 1
781
782 if onlyOverlapping:
783 overlap = astWCS.findWCSOverlap(im1WCS, im2WCS)
784 xOverlap = [overlap['wcs1Pix'][0], overlap['wcs1Pix'][1]]
785 yOverlap = [overlap['wcs1Pix'][2], overlap['wcs1Pix'][3]]
786 xOverlap.sort()
787 yOverlap.sort()
788 xMin = int(math.floor(xOverlap[0] - xBorder))
789 xMax = int(math.ceil(xOverlap[1] + xBorder))
790 yMin = int(math.floor(yOverlap[0] - yBorder))
791 yMax = int(math.ceil(yOverlap[1] + yBorder))
792 xRemainder = (xMax - xMin) % xPixStep
793 yRemainder = (yMax - yMin) % yPixStep
794 if xRemainder != 0:
795 xMax = xMax + xRemainder
796 if yRemainder != 0:
797 yMax = yMax + yRemainder
798
799
800 if xMin < 0:
801 xMin = 0
802 if xMax > im1Data.shape[1]:
803 xMax = im1Data.shape[1]
804 if yMin < 0:
805 yMin = 0
806 if yMax > im1Data.shape[0]:
807 yMax = im1Data.shape[0]
808 else:
809 xMin = 0
810 xMax = im1Data.shape[1]
811 yMin = 0
812 yMax = im1Data.shape[0]
813
814 for x in range(xMin, xMax, xPixStep):
815 for y in range(yMin, yMax, yPixStep):
816 RA, dec = im1WCS.pix2wcs(x, y)
817 x2, y2 = im2WCS.wcs2pix(RA, dec)
818 x2 = int(round(x2))
819 y2 = int(round(y2))
820 if x2 >= 0 and x2 < im2Data.shape[
821 1] and y2 >= 0 and y2 < im2Data.shape[0]:
822 resampledData[y][x] = im2Data[y2][x2]
823
824
825 if not highAccuracy:
826 for row in range(resampledData.shape[0]):
827 vals = resampledData[row, numpy.arange(xMin, xMax, xPixStep)]
828 index2data = interpolate.interp1d(
829 numpy.arange(0, vals.shape[0], 1), vals)
830 interpedVals = index2data(numpy.arange(0, vals.shape[0] - 1, 1.0 /
831 xPixStep))
832 resampledData[row, xMin:xMin + interpedVals.shape[
833 0]] = interpedVals
834 for col in range(resampledData.shape[1]):
835 vals = resampledData[numpy.arange(yMin, yMax, yPixStep), col]
836 index2data = interpolate.interp1d(
837 numpy.arange(0, vals.shape[0], 1), vals)
838 interpedVals = index2data(numpy.arange(0, vals.shape[0] - 1, 1.0 /
839 yPixStep))
840 resampledData[yMin:yMin + interpedVals.shape[0],
841 col] = interpedVals
842
843
844
845
846 return {'data': resampledData, 'wcs': im1WCS.copy()}
847
848
849 -def generateContourOverlay(backgroundImageData, backgroundImageWCS,
850 contourImageData, contourImageWCS, contourLevels,
851 contourSmoothFactor=0, highAccuracy=False):
852 """Rescales an image array to be used as a contour overlay to have the same
853 dimensions as the background image, and generates a set of contour levels.
854 The image array from which the contours are to be generated will be
855 resampled to the same dimensions as the background image data, and can be
856 optionally smoothed using a Gaussian filter. The sigma of the Gaussian
857 filter (contourSmoothFactor) is specified in arcsec.
858
859 @type backgroundImageData: numpy array
860 @param backgroundImageData: background image data array
861 @type backgroundImageWCS: astWCS.WCS
862 @param backgroundImageWCS: astWCS.WCS object of the background image data
863 array
864 @type contourImageData: numpy array
865 @param contourImageData: image data array from which contours are to be
866 generated
867 @type contourImageWCS: astWCS.WCS
868 @param contourImageWCS: astWCS.WCS object corresponding to contourImageData
869 @type contourLevels: list
870 @param contourLevels: sets the contour levels - available options:
871 - values: contourLevels=[list of values specifying each level]
872 - linear spacing: contourLevels=['linear', min level value, max level
873 value, number of levels] - can use "min", "max" to automatically set
874 min, max levels from image data
875 - log spacing: contourLevels=['log', min level value, max level value,
876 number of levels] - can use "min", "max" to automatically set min, max
877 levels from image data
878 @type contourSmoothFactor: float
879 @param contourSmoothFactor: standard deviation (in arcsec) of Gaussian
880 filter for pre-smoothing of contour image data (set to 0 for no smoothing)
881 @type highAccuracy: bool
882 @param highAccuracy: if True, sample every corresponding pixel in each
883 image; otherwise, sample every nth pixel, where n = the ratio of the image
884 scales.
885
886 """
887
888
889
890
891
892 if ("CD1_1" in backgroundImageWCS.header):
893 xScaleFactor = backgroundImageWCS.getXPixelSizeDeg() / (
894 contourImageWCS.getXPixelSizeDeg() / 5.0)
895 yScaleFactor = backgroundImageWCS.getYPixelSizeDeg() / (
896 contourImageWCS.getYPixelSizeDeg() / 5.0)
897 scaledBackground = scaleImage(backgroundImageData, backgroundImageWCS,
898 (xScaleFactor, yScaleFactor))
899 scaled = resampleToWCS(scaledBackground['data'],
900 scaledBackground['wcs'],
901 contourImageData,
902 contourImageWCS,
903 highAccuracy=highAccuracy)
904 scaledContourData = scaled['data']
905 scaledContourWCS = scaled['wcs']
906 scaledBackground = True
907 else:
908 scaled = resampleToWCS(backgroundImageData,
909 backgroundImageWCS,
910 contourImageData,
911 contourImageWCS,
912 highAccuracy=highAccuracy)
913 scaledContourData = scaled['data']
914 scaledContourWCS = scaled['wcs']
915 scaledBackground = False
916
917 if contourSmoothFactor > 0:
918 sigmaPix = (contourSmoothFactor /
919 3600.0) / scaledContourWCS.getPixelSizeDeg()
920 scaledContourData = ndimage.gaussian_filter(scaledContourData,
921 sigmaPix)
922
923
924
925 if contourLevels[0] == "linear":
926 if contourLevels[1] == "min":
927 xMin = contourImageData.flatten().min()
928 else:
929 xMin = float(contourLevels[1])
930 if contourLevels[2] == "max":
931 xMax = contourImageData.flatten().max()
932 else:
933 xMax = float(contourLevels[2])
934 nLevels = contourLevels[3]
935 xStep = (xMax - xMin) / (nLevels - 1)
936 cLevels = []
937 for j in range(nLevels + 1):
938 level = xMin + j * xStep
939 cLevels.append(level)
940
941 elif contourLevels[0] == "log":
942 if contourLevels[1] == "min":
943 xMin = contourImageData.flatten().min()
944 else:
945 xMin = float(contourLevels[1])
946 if contourLevels[2] == "max":
947 xMax = contourImageData.flatten().max()
948 else:
949 xMax = float(contourLevels[2])
950 if xMin <= 0.0:
951 raise Exception(
952 "minimum contour level set to <= 0 and log scaling chosen.")
953 xLogMin = math.log10(xMin)
954 xLogMax = math.log10(xMax)
955 nLevels = contourLevels[3]
956 xLogStep = (xLogMax - xLogMin) / (nLevels - 1)
957 cLevels = []
958 for j in range(nLevels + 1):
959 level = math.pow(10, xLogMin + j * xLogStep)
960 cLevels.append(level)
961
962 else:
963 cLevels = contourLevels
964
965
966 if scaledBackground:
967 scaledBack = scaleImage(scaledContourData, scaledContourWCS, (
968 1.0 / xScaleFactor, 1.0 / yScaleFactor))['data']
969 else:
970 scaledBack = scaledContourData
971
972 return {'scaledImage': scaledBack, 'contourLevels': cLevels}
973
974
975 -def saveBitmap(outputFileName, imageData, cutLevels, size, colorMapName):
976 """Makes a bitmap image from an image array; the image format is specified
977 by the filename extension. (e.g. ".jpg" =JPEG, ".png"=PNG).
978
979 @type outputFileName: string
980 @param outputFileName: filename of output bitmap image
981 @type imageData: numpy array
982 @param imageData: image data array
983 @type cutLevels: list
984 @param cutLevels: sets the image scaling - available options:
985 - pixel values: cutLevels=[low value, high value].
986 - histogram equalisation: cutLevels=["histEq", number of bins ( e.g.
987 1024)]
988 - relative: cutLevels=["relative", cut per cent level (e.g. 99.5)]
989 - smart: cutLevels=["smart", cut per cent level (e.g. 99.5)]
990 ["smart", 99.5] seems to provide good scaling over a range of different
991 images.
992 @type size: int
993 @param size: size of output image in pixels
994 @type colorMapName: string
995 @param colorMapName: name of a standard matplotlib colormap, e.g. "hot",
996 "cool", "gray" etc. (do "help(pylab.colormaps)" in the Python interpreter
997 to see available options)
998
999 """
1000
1001 cut = intensityCutImage(imageData, cutLevels)
1002
1003
1004 aspectR = float(cut['image'].shape[0]) / float(cut['image'].shape[1])
1005 pylab.figure(figsize=(10, 10 * aspectR))
1006 pylab.axes([0, 0, 1, 1])
1007
1008 try:
1009 colorMap = pylab.cm.get_cmap(colorMapName)
1010 except AssertionError:
1011 raise Exception(colorMapName +
1012 " is not a defined matplotlib colormap.")
1013
1014 if cutLevels[0] == "histEq":
1015 pylab.imshow(cut['image'],
1016 interpolation="bilinear",
1017 origin='lower',
1018 cmap=colorMap)
1019
1020 else:
1021 pylab.imshow(cut['image'],
1022 interpolation="bilinear",
1023 norm=cut['norm'],
1024 origin='lower',
1025 cmap=colorMap)
1026
1027 pylab.axis("off")
1028
1029 pylab.savefig("out_astImages.png")
1030 pylab.close("all")
1031
1032 try:
1033 from PIL import Image
1034 except:
1035 raise Exception("astImages.saveBitmap requires the Python Imaging "
1036 "Library to be installed.")
1037 im = Image.open("out_astImages.png")
1038 im.thumbnail((int(size), int(size)))
1039 im.save(outputFileName)
1040
1041 os.remove("out_astImages.png")
1042
1043
1044 -def saveContourOverlayBitmap(outputFileName, backgroundImageData,
1045 backgroundImageWCS, cutLevels,
1046 size, colorMapName, contourImageData,
1047 contourImageWCS,
1048 contourSmoothFactor, contourLevels, contourColor,
1049 contourWidth):
1050 """Makes a bitmap image from an image array, with a set of contours
1051 generated from a second image array overlaid. The image format is specified
1052 by the file extension (e.g. ".jpg"=JPEG, ".png"=PNG). The image array from
1053 which the contours are to be generated can optionally be pre-smoothed using
1054 a Gaussian filter.
1055
1056 @type outputFileName: string
1057 @param outputFileName: filename of output bitmap image
1058 @type backgroundImageData: numpy array
1059 @param backgroundImageData: background image data array
1060 @type backgroundImageWCS: astWCS.WCS
1061 @param backgroundImageWCS: astWCS.WCS object of the background image data
1062 array
1063 @type cutLevels: list
1064 @param cutLevels: sets the image scaling - available options:
1065 - pixel values: cutLevels=[low value, high value].
1066 - histogram equalisation: cutLevels=["histEq", number of bins ( e.g.
1067 1024)]
1068 - relative: cutLevels=["relative", cut per cent level (e.g. 99.5)]
1069 - smart: cutLevels=["smart", cut per cent level (e.g. 99.5)]
1070 ["smart", 99.5] seems to provide good scaling over a range of different
1071 images.
1072 @type size: int
1073 @param size: size of output image in pixels
1074 @type colorMapName: string
1075 @param colorMapName: name of a standard matplotlib colormap, e.g. "hot",
1076 "cool", "gray" etc. (do "help(pylab.colormaps)" in the Python interpreter
1077 to see available options)
1078 @type contourImageData: numpy array
1079 @param contourImageData: image data array from which contours are to be
1080 generated
1081 @type contourImageWCS: astWCS.WCS
1082 @param contourImageWCS: astWCS.WCS object corresponding to contourImageData
1083 @type contourSmoothFactor: float
1084 @param contourSmoothFactor: standard deviation (in pixels) of Gaussian
1085 filter for pre-smoothing of contour image data (set to 0 for no smoothing)
1086 @type contourLevels: list
1087 @param contourLevels: sets the contour levels - available options:
1088 - values: contourLevels=[list of values specifying each level]
1089 - linear spacing: contourLevels=['linear', min level value, max level
1090 value, number of levels] - can use "min", "max" to automatically set
1091 min, max levels from image datad
1092 - log spacing: contourLevels=['log', min level value, max level
1093 value,number of levels] - can use "min", "max" to automatically set
1094 min, max levels from image data
1095 @type contourColor: string
1096 @param contourColor: color of the overlaid contours, specified by the name
1097 of a standard matplotlib color, e.g., "black", "white", "cyan" etc. (do
1098 "help(pylab.colors)" in the Python interpreter to see available options)
1099 @type contourWidth: int
1100 @param contourWidth: width of the overlaid contours
1101
1102 """
1103
1104 cut = intensityCutImage(backgroundImageData, cutLevels)
1105
1106
1107 aspectR = float(cut['image'].shape[0]) / float(cut['image'].shape[1])
1108 pylab.figure(figsize=(10, 10 * aspectR))
1109 pylab.axes([0, 0, 1, 1])
1110
1111 try:
1112 colorMap = pylab.cm.get_cmap(colorMapName)
1113 except AssertionError:
1114 raise Exception(colorMapName +
1115 " is not a defined matplotlib colormap.")
1116
1117 if cutLevels[0] == "histEq":
1118 pylab.imshow(cut['image'],
1119 interpolation="bilinear",
1120 origin='lower',
1121 cmap=colorMap)
1122
1123 else:
1124 pylab.imshow(cut['image'],
1125 interpolation="bilinear",
1126 norm=cut['norm'],
1127 origin='lower',
1128 cmap=colorMap)
1129
1130 pylab.axis("off")
1131
1132
1133 contourData = generateContourOverlay(backgroundImageData,
1134 backgroundImageWCS, contourImageData,
1135 contourImageWCS, contourLevels,
1136 contourSmoothFactor)
1137
1138 pylab.contour(contourData['scaledImage'],
1139 contourData['contourLevels'],
1140 colors=contourColor,
1141 linewidths=contourWidth)
1142
1143 pylab.savefig("out_astImages.png")
1144 pylab.close("all")
1145
1146 try:
1147 from PIL import Image
1148 except ImportError:
1149 raise Exception("astImages.saveContourOverlayBitmap requires the "
1150 "Python Imaging Library to be installed")
1151
1152 im = Image.open("out_astImages.png")
1153 im.thumbnail((int(size), int(size)))
1154 im.save(outputFileName)
1155
1156 os.remove("out_astImages.png")
1157
1158
1159
1160 -def saveFITS(outputFileName, imageData, imageWCS=None):
1161 """Writes an image array to a new .fits file.
1162
1163 @type outputFileName: string
1164 @param outputFileName: filename of output FITS image
1165 @type imageData: numpy array
1166 @param imageData: image data array
1167 @type imageWCS: astWCS.WCS object
1168 @param imageWCS: image WCS object
1169
1170 @note: If imageWCS=None, the FITS image will be written with a rudimentary
1171 header containing no meta data.
1172
1173 """
1174
1175 if os.path.exists(outputFileName):
1176 os.remove(outputFileName)
1177
1178 newImg = pyfits.HDUList()
1179
1180 if imageWCS is not None:
1181 hdu = pyfits.PrimaryHDU(None, imageWCS.header)
1182 else:
1183 hdu = pyfits.PrimaryHDU(None, None)
1184
1185 hdu.data = imageData
1186 newImg.append(hdu)
1187 newImg.writeto(outputFileName)
1188 newImg.close()
1189
1190
1191
1192 -def histEq(inputArray, numBins):
1193 """Performs histogram equalisation of the input numpy array.
1194
1195 @type inputArray: numpy array
1196 @param inputArray: image data array
1197 @type numBins: int
1198 @param numBins: number of bins in which to perform the operation (e.g. 1024)
1199 @rtype: numpy array
1200 @return: image data array
1201
1202 """
1203
1204 imageData = inputArray
1205
1206
1207
1208 sortedDataIntensities = numpy.sort(numpy.ravel(imageData))
1209
1210
1211
1212 dataCumHist = numpy.zeros(numBins)
1213 minIntensity = sortedDataIntensities.min()
1214 maxIntensity = sortedDataIntensities.max()
1215 histRange = maxIntensity - minIntensity
1216 binWidth = histRange / float(numBins - 1)
1217 for i in range(len(sortedDataIntensities)):
1218 binNumber = int(math.ceil((sortedDataIntensities[i] - minIntensity) /
1219 binWidth))
1220 addArray = numpy.zeros(numBins)
1221 onesArray = numpy.ones(numBins - binNumber)
1222 onesRange = list(range(binNumber, numBins))
1223 numpy.put(addArray, onesRange, onesArray)
1224 dataCumHist = dataCumHist + addArray
1225
1226
1227 idealValue = dataCumHist.max() / float(numBins)
1228 idealCumHist = numpy.arange(idealValue, dataCumHist.max() + idealValue,
1229 idealValue)
1230
1231
1232 for y in range(imageData.shape[0]):
1233 for x in range(imageData.shape[1]):
1234
1235 intensityBin = int(math.ceil((imageData[y][x] - minIntensity) /
1236 binWidth))
1237
1238
1239 if intensityBin < 0:
1240 intensityBin = 0
1241 if intensityBin > len(dataCumHist) - 1:
1242 intensityBin = len(dataCumHist) - 1
1243
1244
1245 dataCumFreq = dataCumHist[intensityBin]
1246
1247
1248 idealBin = numpy.searchsorted(idealCumHist, dataCumFreq)
1249 idealIntensity = (idealBin * binWidth) + minIntensity
1250 imageData[y][x] = idealIntensity
1251
1252 return imageData
1253
1254
1255
1257 """Clips the inputArray in intensity and normalises the array such that
1258 minimum and maximum values are 0, 1. Clip in intensity is specified by
1259 clipMinMax, a list in the format [clipMin, clipMax]
1260
1261 Used for normalising image arrays so that they can be turned into RGB
1262 arrays that matplotlib can plot (see L{astPlots.ImagePlot}).
1263
1264 @type inputArray: numpy array
1265 @param inputArray: image data array
1266 @type clipMinMax: list
1267 @param clipMinMax: [minimum value of clipped array, maximum value of
1268 clipped array]
1269 @rtype: numpy array
1270 @return: normalised array with minimum value 0, maximum value 1
1271
1272 """
1273 clipped = inputArray.clip(clipMinMax[0], clipMinMax[1])
1274 slope = 1.0 / (clipMinMax[1] - clipMinMax[0])
1275 intercept = -clipMinMax[0] * slope
1276 clipped = clipped * slope + intercept
1277
1278 return clipped
1279