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