1 """
2 Some utility functions to deal with FITS files.
3
4 Note: pyfits is not thread-safe at least up to version 3.0.8. We therefore
5 provide the fitsLock context manager here that you should use to protect
6 places where you use pyfits in a core (or a similar spot).
7 """
8
9
10
11
12
13
14
15
16
17
18
19
20 import datetime
21 import gzip
22 import itertools
23 import os
24 import re
25 import struct
26 import tempfile
27 import threading
28 import warnings
29 from contextlib import contextmanager
30
31 from . import codetricks
32 from . import excs
33 from . import misctricks
34 from . import ostricks
35
36 try:
37 import numpy
38 import __builtin__
39
40 __builtin__._ASTROPY_SETUP_ = True
41
42 from astropy.io import fits as pyfits
43 from astropy.utils.exceptions import AstropyDeprecationWarning
44 warnings.filterwarnings('ignore', category=AstropyDeprecationWarning)
45 except ImportError:
46
47
48 pyfits = misctricks.NotInstalledModuleStub(
49 "astropy and/or numpy")
50
51 else:
55 raise excs.ReportableError("Incompatible pyfits version."
56 " Please complain to the maintainers.")
57
58 if hasattr(pyfits, "core") and hasattr(pyfits.core, "_TempHDU"):
59 _TempHDU = pyfits.core._TempHDU
60 elif hasattr(pyfits, "_TempHDU"):
61 _TempHDU = pyfits._TempHDU
62 elif hasattr(pyfits.Header, "fromstring"):
64 """a wrapper around modern pyfits to provide some ancient whacko
65 functionality."""
68
70 self.header = pyfits.Header.fromstring(self._raw)
71 return self
72
73
74
75
76 try:
77 from pyfits import core as pyfitscore
78 warnings.showwarning = pyfitscore.showwarning
79 warnings.formatwarning = pyfitscore.formatwarning
80 except (ImportError, AttributeError) as ex:
81
82 pass
83
84
85
86
87
88
89
90 if not hasattr(pyfits.Header, "ascardlist"):
93 pyfits.Header.ascardlist = _ascardlist
94 del _ascardlist
95
96 if not hasattr(pyfits.Header, "append"):
97 - def _append(self, card, end=False):
98 self.ascard.append(card, bottom=end)
99 pyfits.Header.append = _append
100 del _append
101
102 _FITS_TABLE_LOCK = threading.RLock()
111
112
113 CARD_SIZE = 80
114
115 END_CARD = 'END'+' '*(CARD_SIZE-3)
116
117 FITS_BLOCK_SIZE = CARD_SIZE*36
122
123
124
125
126 warnings.filterwarnings('ignore', category=DeprecationWarning)
127 try:
128 warnings.filterwarnings('ignore', category=pyfits.PyfitsDeprecationWarning)
129 except AttributeError:
130 pass
133 """pads input (a string) with blanks until len(result)%80=0
134
135 The length keyword argument lets you adjust the "card size". Use
136 this to pad headers with length=FITS_BLOCK_SIZE
137
138 >>> padCard("")
139 ''
140 >>> len(padCard("end"))
141 80
142 >>> len(padCard("whacko"*20))
143 160
144 >>> len(padCard("junkodumnk"*17, 27))%27
145 0
146 """
147
148
149 l = len(input)
150 if not l%length:
151 return input
152 return input+' '*(length-l%length)
153
156 """returns the bytes beloning to a FITS header starting at the current
157 position within the file f.
158
159 If the header is not complete after reading maxHeaderBlocks blocks,
160 a FITSError is raised.
161 """
162 parts = []
163
164 while True:
165 block = f.read(FITS_BLOCK_SIZE)
166 if not block:
167 raise EOFError('Premature end of file while reading header')
168
169 parts.append(block)
170 endCardPos = block.find(END_CARD)
171 if not endCardPos%CARD_SIZE:
172 break
173
174 if len(parts)>=maxHeaderBlocks:
175 raise FITSError("No end card found within %d blocks"%maxHeaderBlocks)
176 return "".join(parts)
177
180 """returns a pyfits header for the primary hdu of the opened file f.
181
182 This is mostly code lifted from pyfits._File._readHDU. The way
183 that class is made, it's hard to use it with stuff from a gzipped
184 source, and that's why this function is here. It is used in the quick
185 mode of fits grammars.
186
187 This function is adapted from pyfits.
188 """
189 hdu = _TempHDU()
190 hdu._raw = readHeaderBytes(f, maxHeaderBlocks)
191 hdu._extver = 1
192
193 hdu._new = 0
194 hdu = hdu.setupHDU()
195 return hdu.header
196
199 """returns a list of pyfits Cards parsed from aString.
200
201 This will raise a ValueError if aString's length is not divisible by 80.
202 It may also raise pyfits errors for malformed cards.
203
204 Empty (i.e., all-whitespace) cards are ignored. If an END card is
205 encountered processing is aborted.
206 """
207 cards = []
208 if len(aString)%CARD_SIZE:
209 raise ValueError("parseCards argument has impossible length %s"%(
210 len(aString)))
211 for offset in range(0, len(aString), CARD_SIZE):
212 rawCard = aString[offset:offset+CARD_SIZE]
213 if rawCard==END_CARD:
214 break
215 if not rawCard.strip():
216 continue
217 cards.append(pyfits.Card().fromstring(rawCard))
218 return cards
219
222 """returns the FITS serialization of a FITS header hdr.
223 """
224 parts = []
225 for card in hdr.ascardlist():
226 r = card.image
227 assert not len(r)%CARD_SIZE
228 parts.append(r)
229 serForm = "".join(parts)+padCard('END')
230 return padCard(serForm, length=FITS_BLOCK_SIZE)
231
234 """writes a FITS to targetFile having newHeader as the primary header,
235 where the rest of the data is taken from inputFile.
236
237 inputFile must be a file freshly opened for reading, targetFile one
238 freshly opened for writing.
239
240 This function is (among other things) a workaround for pyfits' misfeature of
241 unscaling scaled data in images when extending a header.
242 """
243 readPrimaryHeaderQuick(inputFile)
244 targetFile.write(serializeHeader(newHeader))
245 while True:
246 buf = inputFile.read(bufSize)
247 if not buf:
248 break
249 targetFile.write(buf)
250
253 """replaces the primary header of fitsName with newHeader.
254
255 Doing this, it tries to minimize the amount of writing necessary; if
256 fitsName has enough space for newHeader, just the header is written,
257 and newHeader is extended if necessary. Only if newHeader is longer than
258 the existing header is fitsName actually copied. We try to be safe in
259 this case, only overwriting the old entry when the new data is safely
260 on disk.
261
262 gzipped inputs used to be supported here, but they aren't any more.
263 """
264 if fitsName.endswith(".gz"):
265 raise NotImplementedError("replacePrimaryHeaderInPlace no longer"
266 " supports gzipped files.")
267
268 serializedNew = serializeHeader(newHeader)
269 with open(fitsName) as inputFile:
270 serializedOld = readHeaderBytes(inputFile)
271 inputFile.seek(0)
272
273 if len(serializedNew)<len(serializedOld):
274
275
276 serializedNew = serializedNew+(
277 len(serializedOld)-len(serializedNew))*" "
278 serializedNew = serializedNew.replace(END_CARD, " "*len(END_CARD))
279 serializedNew = serializedNew[:-len(END_CARD)]+END_CARD
280 assert len(serializedNew)==len(serializedOld)
281
282 if len(serializedNew)==len(serializedOld):
283
284
285 with open(fitsName, "r+") as targetFile:
286 targetFile.seek(0)
287 targetFile.write(serializedNew)
288
289 else:
290
291 with ostricks.safeReplaced(fitsName) as targetFile:
292 replacePrimaryHeader(inputFile, newHeader, targetFile)
293
294
295
296 STANDARD_CARD_SEQUENCE = [
297 ("SIMPLE", True),
298 ("BITPIX", True),
299 ("NAXIS", True),
300 ("NAXIS1", False),
301 ("NAXIS2", False),
302 ("NAXIS3", False),
303 ("NAXIS4", False),
304 ("EXTEND", False),
305 ("BZERO", False),
306 ("BSCALE", False),
307 ]
310 """helps _enforceHeaderConstraints.
311 """
312 if seq is None:
313 return
314 for item in seq:
315 if isinstance(item, tuple):
316 yield item
317 else:
318 yield (item, False)
319
322 """returns a pyfits header containing the cards in cardList with FITS
323 sequence constraints satisfied.
324
325 This may raise a FITSError if mandatory cards are missing.
326
327 This will only work correctly for image FITSes with less than five
328 dimensions.
329
330 On cardSequence, see sortHeaders (except that this function always
331 requries sequences of pairs).
332 """
333
334
335 cardsAdded, newCards = set(), []
336 cardDict = dict((card.keyword, card) for card in cardList)
337
338 for kw, mandatory in itertools.chain(
339 STANDARD_CARD_SEQUENCE, _iterForcedPairs(cardSequence)):
340 if isinstance(kw, pyfits.Card):
341 newCards.append(kw)
342 continue
343
344 if kw in cardsAdded:
345 continue
346
347 try:
348 newCards.append(cardDict[kw])
349 cardsAdded.add(kw)
350 except KeyError:
351 if mandatory:
352 raise FITSError("Mandatory card '%s' missing"%kw)
353
354 for card in cardList:
355
356 if card.keyword not in cardsAdded or card.keyword=='':
357 newCards.append(card)
358 return pyfits.Header(newCards)
359
363 """returns a pyfits header with "real" cards first, then history, then
364 comment cards.
365
366 Blanks in the input are discarded, and one blank each is added in
367 between the sections of real cards, history and comments.
368
369 Header can be an iterable yielding Cards or a pyfits header.
370
371 Duplicate history or comment entries will be swallowed.
372
373 cardSequence, if present, is a sequence of (item, mandatory) pairs.
374 There item can be a card name, in which case the corresponding
375 card will be inserted at this point in the sequence. If mandatory is
376 True, a missing card is an error. Keywords already in
377 fitstools.STANDARD_CARD_SEQUENCE are ignored.
378
379 Item can also a pyfits.Card instance; it will be put into the header
380 as-is.
381
382 As a shortcut, a sequence item may be something else then a tuple;
383 it will then be combined with a False to make one.
384
385 These days, if you think you need this, have a look at
386 fitstricks.makeHeaderFromTemplate first.
387 """
388 commentCs, historyCs, realCs = [], [], []
389 if hasattr(header, "ascardlist"):
390 iterable = header.ascardlist()
391 else:
392 iterable = header
393 for card in iterable:
394 if card.keyword=="COMMENT":
395 commentCs.append(card)
396 elif card.keyword=="HISTORY":
397 historyCs.append(card)
398 else:
399 realCs.append(card)
400
401 newCards = []
402 for card in realCs:
403 newCards.append(card)
404
405 historySeen = set()
406 if historyCs:
407 newCards.append(pyfits.Card(keyword=""))
408 for card in historyCs:
409 if historyFilter is None or historyFilter(card.value):
410 if card.value not in historySeen:
411 newCards.append(card)
412 historySeen.add(card.value)
413
414 commentsSeen = set()
415 if commentCs:
416 newCards.append(pyfits.Card(keyword=""))
417 for card in commentCs:
418 if commentFilter is None or commentFilter(card.value):
419 if card.value not in commentsSeen:
420 commentsSeen.add(card.value)
421 newCards.append(card)
422
423 return _enforceHeaderConstraints(newCards, cardSequence)
424
425
426 -def openGz(fitsName, tempDir=None):
427 """returns the hdus for the gzipped fits fitsName.
428
429 Scrap that as soon as we have gzipped fits support (i.e. newer pyfits)
430 in debian.
431 """
432 handle, pathname = tempfile.mkstemp(suffix="fits", dir=tempDir)
433 f = os.fdopen(handle, "w")
434 f.write(gzip.open(fitsName).read())
435 f.close()
436 hdus = pyfits.open(pathname)
437 hdus.readall()
438 os.unlink(pathname)
439 return hdus
440
441
442 -def writeGz(hdus, fitsName, compressLevel=5, mode=0664):
443 """writes and gzips hdus into fitsName. As a side effect, hdus will be
444 closed.
445
446 Appearently, not even recent versions of pyfits support writing of
447 zipped files (which is a bit tricky, admittedly). So, we'll probably
448 have to live with this kludge for a while.
449 """
450 handle, pathname = tempfile.mkstemp(suffix="fits")
451 with codetricks.silence():
452 hdus.writeto(pathname, clobber=True)
453 os.close(handle)
454 rawFitsData = open(pathname).read()
455 os.unlink(pathname)
456 handle, pathname = tempfile.mkstemp(suffix="tmp",
457 dir=os.path.dirname(fitsName))
458 os.close(handle)
459 dest = gzip.open(pathname, "w", compressLevel)
460 dest.write(rawFitsData)
461 dest.close()
462 os.rename(pathname, fitsName)
463 os.chmod(fitsName, mode)
464
467 """returns the hdus for fName.
468
469 (gzip detection is tacky, and we should look at the magic).
470
471 This is simulating legacy pyfits behaviour in that it will use
472 memmap I/O if it thinks that works. You'll probably want to use
473 astrpy directly these days if this gives you a headache.
474 """
475 if os.path.splitext(fitsName)[1].lower()==".gz":
476 return openGz(fitsName)
477 else:
478 hdus = pyfits.open(fitsName, memmap=True)
479
480
481 if (hdus[0].header.get("BZERO", 0)!=0
482 or hdus[0].header.get("BSCALE", 1)!=1
483 or "BLANK" in hdus[0].header):
484 hdus = pyfits.open(fitsName, memmap=False)
485 return hdus
486
489 """A class that allows header manipulation of fits files
490 without having to touch the data.
491
492 This class exists because pyfits insists on scaling scaled image data
493 on reading it. While we can scale back on writing, this is something
494 I'd rather not do. So, I have this base class to facilate the
495 HeaderManipulator that can handle gzipped fits files as well.
496 """
498 self.hdus = pyfits.open(fName, "update")
499 self.add_comment = self.hdus[0].header.add_comment
500 self.add_history = self.hdus[0].header.add_history
501 self.add_blank = self.hdus[0].header.add_blank
502 self.update = self.hdus[0].header.set
503 self.set = self.hdus[0].header.set
504
508
511
514 """is a class that allows header manipulation of fits files without
515 having to touch the data even for gzipped files.
516
517 See PlainHeaderManipulator. We only provide a decoration here that
518 transparently gzips and ungzips compressed fits files.
519 """
521 self.origFile = fName
522 handle, self.uncompressedName = tempfile.mkstemp(
523 suffix="fits")
524 destFile = os.fdopen(handle, "w")
525 destFile.write(gzip.open(fName).read())
526 destFile.close()
527 self.compressLevel = compressLevel
528 PlainHeaderManipulator.__init__(self, self.uncompressedName)
529
531 PlainHeaderManipulator.close(self)
532 destFile = gzip.open(self.origFile, "w", compresslevel=self.compressLevel)
533 destFile.write(open(self.uncompressedName).read())
534 destFile.close()
535 os.unlink(self.uncompressedName)
536
539 """returns a header manipulator for a FITS file.
540
541 (it's supposed to look like a class, hence the uppercase name)
542 It should automatically handle gzipped files.
543 """
544 if fName.lower().endswith(".gz"):
545 return GzHeaderManipulator(fName)
546 else:
547 return PlainHeaderManipulator(fName)
548
551 """returns x and y size a fits image.
552 """
553 hdr = readPrimaryHeaderQuick(open(fName))
554 return hdr["NAXIS1"], hdr["NAXIS2"]
555
558 """returns a sequence of the lengths of the axes of a FITS image
559 described by hdr.
560 """
561 return [hdr["NAXIS%d"%i] for i in range(1, hdr["NAXIS"]+1)]
562
565 """returns a array slices and a new pyfits header reflecting cuts.
566
567 (see cutoutFits for what cuts is).
568 """
569 cutDict = dict((c[0], c[1:]) for c in cuts)
570 slices = []
571 newHeader = header.copy()
572
573 for index, length in enumerate(getAxisLengths(header)):
574 firstPix, lastPix = cutDict.get(index+1, (None, None))
575
576 if firstPix is None:
577 firstPix = 1
578 if lastPix is None:
579 lastPix = length
580 firstPix = min(max(1, firstPix), length)
581 lastPix = min(length, max(1, lastPix))
582
583 if (firstPix, lastPix)==(1, length):
584 slices.append(slice(None, None, None))
585 else:
586 firstPix -= 1
587 newAxisLength = lastPix-firstPix
588 if newAxisLength==0:
589 newAxisLength = 1
590 lastPix = firstPix+1
591 slices.append(slice(int(firstPix), int(lastPix), 1))
592
593 newHeader["NAXIS%d"%(index+1)] = newAxisLength
594 refpixKey = "CRPIX%d"%(index+1)
595 newHeader.set(refpixKey, newHeader[refpixKey]-firstPix)
596
597 return slices, newHeader
598
601 """returns a cutout of hdu restricted to cuts.
602
603 hdu is a primary FITS hdu. cuts is a list of cut specs, each of which is
604 a triple (axis, lower, upper). axis is between 1 and naxis, lower and
605 upper a 1-based pixel coordinates of the limits, and "border" pixels
606 are included. Specifications outside of the image are legal and will
607 be cropped back. Open limits are supported via a specification of
608 None.
609
610 If an axis would vanish (i.e. length 0 or less), the function fudges
611 things such that the axis gets a length of 1.
612
613 axis is counted here in the FORTRAN/FITS sense, *not* in the C sense,
614 i.e., axis=1 cuts along NAXIS1, which is the *last* index in a numpy
615 array.
616
617 WCS CRPIXes in hdu's header will be updated. Axes and specified will
618 not be touched. It is an error to specifiy cuts for an axis twice
619 (behaviour is undefined).
620
621 Note that this will lose all extensions the orginal FITS file might have
622 had.
623 """
624 slices, newHeader = cutoutHeader(hdu.header, *cuts)
625 slices.reverse()
626 newHDU = hdu.__class__(data=hdu.data[tuple(slices)].copy(order='C'),
627 header=newHeader)
628
629 if "BZERO" in newHeader:
630 newHDU.header["BZERO"] = newHeader["BZERO"]
631 if "BSCALE" in newHeader:
632 newHDU.header["BSCALE"] = newHeader["BSCALE"]
633
634 return newHDU
635
638 """returns a FITS header suitable for a shrunken version of the image
639 described by oldHeader.
640
641 This only works for 2d images, scale must be an integer>1. The function
642 assumes no "fractional" pixels are handled, i.e, remainders in divisions
643 with factors are discarded. It is thus a companion for
644 iterScaledRows.
645
646 Note that oldHeader must be an actual pyfits header instance; a dictionary
647 will not do.
648
649 This is a pretty straight port of wcstools's imutil.c#ShrinkFITSHeader,
650 except we clear BZERO and BSCALE and set BITPIX to -32 (float array)
651 under the assumption that in the returned image, 32-bit floats are used.
652 """
653 assert oldHeader["NAXIS"]==2
654
655 factor = int(factor)
656 newHeader = oldHeader.copy()
657 newHeader.set("SIMPLE", True,"GAVO DaCHS, %s"%datetime.datetime.utcnow())
658 newHeader["NAXIS1"] = oldHeader["NAXIS1"]//factor
659 newHeader["NAXIS2"] = oldHeader["NAXIS2"]//factor
660 newHeader["BITPIX"] = -32
661
662 try:
663 ffac = float(factor)
664 newHeader["CRPIX1"] = oldHeader["CRPIX1"]/ffac+0.5
665 newHeader["CRPIX2"] = oldHeader["CRPIX2"]/ffac+0.5
666 for key in ("CDELT1", "CDELT2",
667 "CD1_1", "CD2_1", "CD1_2", "CD2_2"):
668 if key in oldHeader:
669 newHeader[key] = oldHeader[key]*ffac
670 except KeyError:
671 pass
672
673 newHeader.set("IMSHRINK", "Image scaled down %s-fold by DaCHS"%factor)
674
675 for hField in ["BZERO", "BSCALE"]:
676 if hField in newHeader:
677 del newHeader[hField]
678
679 return newHeader
680
681
682 NUM_CODE = {
683 8: 'uint8',
684 16: '>i2',
685 32: '>i4',
686 64: '>i8',
687 -32: '>f4',
688 -64: '>f8'}
691 """returns a decoder for the rows of FITS primary image data.
692
693 The decoder is called with an open file and returns the next row.
694 You need to keep track of the total number of rows yourself.
695 """
696 numType = NUM_CODE[hdr["BITPIX"]]
697 rowLength = hdr["NAXIS1"]
698
699 bzero, bscale = hdr.get("BZERO", 0), hdr.get("BSCALE", 1)
700 if bzero!=0 or bscale!=1:
701 def read(f):
702 return numpy.asarray(
703 numpy.fromfile(f, numType, rowLength), 'float32')*bscale+bzero
704 else:
705 def read(f):
706 return numpy.fromfile(f, numType, rowLength)
707
708 return read
709
712 """iterates over the rows of a FITS (primary) image.
713
714 You pass in a FITS header and a file positioned to the start of
715 the image data.
716
717 What's returned are 1d numpy arrays of the datatype implied by bitpix. The
718 function understands only very basic FITSes (BSCALE and BZERO are known,
719 though, and lead to floats arrays).
720
721 We do this ourselves since pyfits may pull in the whole thing or at least
722 mmaps it; both are not attractive when I want to stream-process large
723 images.
724 """
725 decoder = _makeDecoder(hdr)
726 for col in xrange(hdr["NAXIS2"]):
727 yield decoder(f)
728
731 """helps iterScaledRows for the case of a simple, real-file FITS.
732 """
733 if hdr is None:
734 hdr = readPrimaryHeaderQuick(inFile)
735 if hdr["NAXIS"]==0:
736
737 return _iterSetupCompatible(inFile, hdr)
738 return hdr, iterFITSRows(hdr, inFile)
739
742 """returns the index of a compressed image HDU if it immediately
743 follows extInd.
744
745 If hdus[extInd] is a CompImageHDU itself, this will return extInd
746 unchanged.
747 """
748 if isinstance(hdus[extInd], pyfits.CompImageHDU):
749 return extInd
750 try:
751 if isinstance(hdus[extInd+1], pyfits.CompImageHDU):
752 return extInd+1
753 except IndexError:
754
755 pass
756 return extInd
757
760 """helps iterScaledRows for when _iterSetupFast will not work.
761
762 Using extInd, you can select a different extension. extInd=0
763 will automatically select extension 1 if that's a compressed image
764 HDU.
765 """
766 hdus = pyfits.open(inFile)
767 extInd = fixImageExtind(hdus, extInd)
768
769 def iterRows():
770 for row in hdus[extInd].data[:,:]:
771 yield row
772
773 return hdus[extInd].header, iterRows()
774
775
776 -def iterScaledRows(inFile, factor=None, destSize=None, hdr=None, slow=False,
777 extInd=0):
778 """iterates over numpy arrays of pixel rows within the open FITS
779 stream inFile scaled by it integer in factor.
780
781 The arrays are always float32, regardless of the input. When the
782 image size is not a multiple of factor, border pixels are discarded.
783
784 A FITS header for this data can be obtained using shrinkWCSHeader.
785
786 You can pass in either a factor the file is to be scaled down, or
787 an approximate size desired. If both are given, factor takes precedence,
788 if none is given, it's an error.
789
790 If you pass in hdr, it will not be read from inFile; the function then
791 expects the file pointer to point to the start of the first data block.
792 Use this if you've already read the header of a non-seekable FITS.
793
794 extInd lets you select a different extension. extInd=0 means the first
795 image HDU, which may be in extension 1 for compressed images.
796
797 iterScaledRows will try to use a hand-hacked interface guaranteed to
798 stream. This only works for plain, 2D-FITSes from real files.
799 iterScaledRows normally notices when it should fall back to
800 pyfits code, but if it doesn't you can pass slow=True.
801 """
802 if isinstance(inFile, file) and not slow and extInd==0:
803 hdr, rows = _iterSetupFast(inFile, hdr)
804 else:
805 hdr, rows = _iterSetupCompatible(inFile, hdr, extInd)
806
807 if factor is None:
808 if destSize is None:
809 raise excs.DataError("iterScaledRows needs either factor or destSize.")
810 size = max(hdr["NAXIS1"], hdr["NAXIS2"])
811 factor = max(1, size//destSize+1)
812
813 factor = int(factor)
814 assert factor>0
815
816 rowLength = hdr["NAXIS1"]
817 destRowLength = rowLength//factor
818 summedInds = range(factor)
819
820 for index in xrange(hdr["NAXIS2"]//factor):
821 newRow = numpy.zeros((rowLength,), 'float32')
822 for i in summedInds:
823 try:
824 newRow += rows.next()
825 except StopIteration:
826 break
827 newRow /= factor
828
829
830
831 newRow = newRow[:destRowLength*factor]
832 yield sum(numpy.transpose(
833 (newRow/factor).reshape((destRowLength, factor))))
834
837 """iterates over the bytes for a simple FITS file generated by scaling
838 down the 2D image inFileName by factor.
839
840 factor must be an integer between 1 and something reasonable.
841
842 This function acquires the FITS lock itself.
843 """
844 with fitsLock():
845 with open(inFileName) as f:
846 oldHdr = readPrimaryHeaderQuick(f)
847 newHdr = shrinkWCSHeader(oldHdr, factor)
848 newHdr.update(extraCards)
849 yield serializeHeader(newHdr)
850
851 for row in iterScaledRows(f, factor, hdr=oldHdr):
852
853
854
855
856 yield struct.pack("!%df"%len(row), *row)
857
860 """returns a primary header sporting the key/value pairs given in the
861 dictionary d.
862
863 In all likelihood, this header will *not* work properly as a primary
864 header because, e.g., there are certain rules on the sequence of
865 header items. fitstricks.copyFields can make a valid header out
866 of what's returned here.
867
868 keys mapped to None are skipped, i.e., you have to do nullvalue handling
869 yourself.
870 """
871 hdr = pyfits.PrimaryHDU().header
872 for key, value in d.iteritems():
873 if value is not None:
874 hdr.set(key, value)
875 return hdr
876
879 """represents a single 1D WCS axis and allows easy metadata discovery.
880
881 You'll usually use the fromHeader constructor.
882
883 The idea of using this rather than astropy.wcs or similar is that this is
884 simple and robust. It doesn't know many of the finer points of WCS,
885 though, and in particular it's 1D only.
886
887 However, for the purposes of cutouts it probably should do for the
888 overwhelming majority of non-spatial FITS axes.
889
890 The default pixel coordinates are handled in the FITS sense here,
891 i.e., the first pixel has the index 1. Three are methods that have
892 pix0 in their names; these assume 0-based arrays. All the transforms
893 return Nones unchanged.
894
895 To retrieve the metadata shoved in, use the name, crval, crpix, cdelt,
896 ctype, cunit, and axisLength attributes.
897 """
898 - def __init__(self, name, crval, crpix, cdelt,
899 ctype="UNKNOWN", cunit="", axisLength=1):
900 assert cdelt!=0
901 self.crval, self.crpix, self.cdelt = crval, crpix, cdelt
902 self.ctype, self.cunit, self.axisLength = ctype, cunit, axisLength
903 self.name = name
904
906 """returns the physical value for a 1-based pixel coordinate.
907 """
908 if pixCoo is None:
909 return None
910 return self.crval+(pixCoo-self.crpix)*self.cdelt
911
913 """returns the physical value for a 0-based pixel coordinate.
914 """
915 if pix0Coo is None:
916 return None
917 return self.pixToPhys(pix0Coo+1)
918
920 """returns a 1-based pixel coordinate for a physical value.
921 """
922 if physCoo is None:
923 return None
924 return (physCoo-self.crval)/self.cdelt+self.crpix
925
927 """returns a 0-based pixel coordinate for a physical value.
928 """
929 if physCoo is None:
930 return None
931 return self.physToPix(physCoo)-1
932
934 """returns the minimal and maximal physical values this axis
935 takes within the image.
936 """
937 limits = self.pixToPhys(1), self.pixToPhys(self.axisLength)
938 return min(limits), max(limits)
939
940 @classmethod
942 """returns a WCSAxis for the specified axis in header.
943
944 If the axis is mentioned in a transformation matrix (CD or PC),
945 a ValueError is raised; this is strictly for 1D coordinates.
946
947 The axisIndex is 1-based; to get a transform for the axis described
948 by CTYPE1, pass 1 here.
949 """
950 if ("CD%d_%d"%(axisIndex, axisIndex) in header
951 or "PC%d_%d"%(axisIndex, axisIndex) in header):
952 if not forceSeparable:
953 raise ValueError("FITS axis %s appears not separable. WCSAxis"
954 " cannot handle this."%axisIndex)
955
956 def get(key, default):
957 return header.get("%s%d"%(key, axisIndex), default)
958
959 guessedName = get("CNAME", "").strip()
960 if not guessedName:
961 guessedName = get("CTYPE", "").split("-", 1)[0].strip()
962 if not guessedName or guessedName=="UNKNOWN":
963 guessedName = "COO"
964
965 guessedName = "%s_%d"%(guessedName, axisIndex)
966
967 return cls(guessedName,
968 get("CRVAL", 0), get("CRPIX", 0), get("CDELT", 1),
969 get("CTYPE", "UNKNOWN").strip(),
970 get("CUNIT", "").strip() or None,
971 get("NAXIS", 1))
972
977 """returns a WCSAxis instance from an axis index and a FITS header.
978
979 If the axis is mentioned in a transformation matrix (CD or PC),
980 a ``ValueError`` is raised (use ``forceSeparable`` to override).
981
982 The ``axisIndex`` is 1-based; to get a transform for the axis described
983 by CTYPE1, pass 1 here.
984
985 The object returned has methods like ``pixToPhys``, ``physToPix`` (and their
986 ``pix0`` brethren), and ``getLimits``.
987
988 Note that at this point WCSAxis only supports linear transforms (it's
989 a DaCHS-specific implementation). We'll extend it on request.
990 """
991 return WCSAxis.fromHeader(header, axisIndex, forceSeparable)
992
995 """is raised when something goes wrong while parsing ESO descriptors.
996 """
997
1000 return "(?P<%s>%s)"%(name, re)
1001
1004 """an ad-hoc parser for ESO's descriptors.
1005
1006 These are sometimes in FITS files produced by MIDAS. What I'm pretending
1007 to parse here is the contatenation of the cards' values without the
1008 boundaries.
1009
1010 The parse is happening at construction time, after which you fetch the result
1011 in the result attribute.
1012
1013 I'm adhoccing this. If someone digs up the docs on what the actual
1014 grammar is, I'll do it properly.
1015 """
1016
1017 stringPat = "'[^']*'"
1018 intPat = r"-?\d+"
1019 floatPat = r"-?\d+\.\d+E[+-]\d+"
1020 white = r"\s+"
1021 nextToken = r"\s*"
1022 headerSep = nextToken+','+nextToken
1023
1024 headerRE = re.compile(nextToken
1025 +_makeNamed("colName", stringPat)+headerSep
1026 +_makeNamed("typeCode", stringPat)+headerSep
1027 +_makeNamed("startInd", intPat)+headerSep
1028 +_makeNamed("endInd", intPat)+headerSep
1029 +stringPat+headerSep
1030 +stringPat+headerSep
1031 +stringPat)
1032 floatRE = re.compile(nextToken+floatPat)
1033 integerRE = re.compile(nextToken+intPat)
1034
1036 self.data, state = data, "header"
1037 self.result, self.curPos = {}, 0
1038 while state!="end":
1039 try:
1040 state = getattr(self, "_scan_"+state)()
1041 except Exception as msg:
1042 raise ESODescriptorsError(str(msg),
1043 location="character %s"%self.curPos,
1044 offending=repr(self.data[self.curPos:self.curPos+20]))
1045
1046 @classmethod
1048 """returns a parser from the data in the pyfits hdr.
1049 """
1050 descLines, collecting = [], False
1051
1052 for card in hdr.cards:
1053 if card.keyword=="HISTORY":
1054 if " ESO-DESCRIPTORS END" in card.value:
1055 collecting = False
1056 if collecting:
1057 descLines.append(card.value)
1058 if " ESO-DESCRIPTORS START" in card.value:
1059 collecting = True
1060
1061 return cls("\n".join(descLines))
1062
1064 """read the next descriptor header.
1065 """
1066 mat = self.headerRE.match(self.data, self.curPos)
1067 if not mat:
1068 if not self.data[self.curPos:].strip():
1069 return "end"
1070 else:
1071 raise ValueError("Could not find next header")
1072
1073 self.curPos = mat.end()
1074 self.curCol = []
1075 self.curColName = mat.group("colName")[1:-1]
1076 self.yetToRead = int(mat.group("endInd"))-int(mat.group("startInd"))+1
1077
1078 if mat.group("typeCode")[1:-1].startswith("R"):
1079 return "float"
1080 elif mat.group("typeCode")[1:-1].startswith("I"):
1081 return "integer"
1082 else:
1083 raise ValueError("Unknown type code %s"%mat.group("typeCode"))
1084
1086 """enter a parsed column into result and prepare for the next column.
1087 """
1088 self.result[self.curColName] = self.curCol
1089 del self.curCol
1090 del self.curColName
1091 del self.yetToRead
1092 return "header"
1093
1095 """returns a scanner function (invisible to the outside.
1096
1097 This is a helper method for building the parser class; typeName
1098 must be the same as the name in _scan_name. constructor is
1099 a function that turns string literals into objects of the desired
1100 type.
1101 """
1102 def scanner(self):
1103 mat = literalRE.match(self.data, self.curPos)
1104 if not mat:
1105 raise ValueError("Expected a %s here"%typeName)
1106 self.curPos = mat.end()
1107
1108 self.curCol.append(constructor(mat.group()))
1109 self.yetToRead -= 1
1110 if self.yetToRead==0:
1111 return "harvest"
1112 else:
1113 return typeName
1114
1115 return scanner
1116
1117 _scan_float = _makeTypeScanner("float", floatRE, float)
1118 _scan_integer = _makeTypeScanner("integer", integerRE, int)
1119
1120 del _makeTypeScanner
1121
1124 """returns parsed ESO descriptors from a pyfits header hdr.
1125
1126 ESO descriptors are data columns stuck into FITS history lines.
1127 They were produced by MIDAS. This implementation was made
1128 without actual documentation, is largely based on conjecture,
1129 and is certainly incomplete.
1130
1131 What's returned is a dictionary mapping column keywords to lists of
1132 column values.
1133 """
1134 return _ESODescriptorsParser.fromFITSHeader(hdr).result
1135
1140
1141 if __name__=="__main__":
1142 _test()
1143