1# Licensed under a 3-clause BSD style license - see PYFITS.rst 2 3import ctypes 4import gc 5import itertools 6import math 7import re 8import time 9import warnings 10from contextlib import suppress 11 12import numpy as np 13 14from .base import DELAYED, ExtensionHDU, BITPIX2DTYPE, DTYPE2BITPIX 15from .image import ImageHDU 16from .table import BinTableHDU 17from astropy.io.fits import conf 18from astropy.io.fits.card import Card 19from astropy.io.fits.column import Column, ColDefs, TDEF_RE 20from astropy.io.fits.column import KEYWORD_NAMES as TABLE_KEYWORD_NAMES 21from astropy.io.fits.fitsrec import FITS_rec 22from astropy.io.fits.header import Header 23from astropy.io.fits.util import (_is_pseudo_integer, _pseudo_zero, _is_int, 24 _get_array_mmap) 25 26from astropy.utils import lazyproperty 27from astropy.utils.exceptions import AstropyUserWarning 28 29try: 30 from astropy.io.fits import compression 31 COMPRESSION_SUPPORTED = COMPRESSION_ENABLED = True 32except ImportError: 33 COMPRESSION_SUPPORTED = COMPRESSION_ENABLED = False 34 35 36# Quantization dithering method constants; these are right out of fitsio.h 37NO_DITHER = -1 38SUBTRACTIVE_DITHER_1 = 1 39SUBTRACTIVE_DITHER_2 = 2 40QUANTIZE_METHOD_NAMES = { 41 NO_DITHER: 'NO_DITHER', 42 SUBTRACTIVE_DITHER_1: 'SUBTRACTIVE_DITHER_1', 43 SUBTRACTIVE_DITHER_2: 'SUBTRACTIVE_DITHER_2' 44} 45DITHER_SEED_CLOCK = 0 46DITHER_SEED_CHECKSUM = -1 47 48COMPRESSION_TYPES = ('RICE_1', 'GZIP_1', 'GZIP_2', 'PLIO_1', 'HCOMPRESS_1') 49 50# Default compression parameter values 51DEFAULT_COMPRESSION_TYPE = 'RICE_1' 52DEFAULT_QUANTIZE_LEVEL = 16. 53DEFAULT_QUANTIZE_METHOD = NO_DITHER 54DEFAULT_DITHER_SEED = DITHER_SEED_CLOCK 55DEFAULT_HCOMP_SCALE = 0 56DEFAULT_HCOMP_SMOOTH = 0 57DEFAULT_BLOCK_SIZE = 32 58DEFAULT_BYTE_PIX = 4 59 60CMTYPE_ALIASES = {'RICE_ONE': 'RICE_1'} 61 62COMPRESSION_KEYWORDS = {'ZIMAGE', 'ZCMPTYPE', 'ZBITPIX', 'ZNAXIS', 'ZMASKCMP', 63 'ZSIMPLE', 'ZTENSION', 'ZEXTEND'} 64 65 66class CompImageHeader(Header): 67 """ 68 Header object for compressed image HDUs designed to keep the compression 69 header and the underlying image header properly synchronized. 70 71 This essentially wraps the image header, so that all values are read from 72 and written to the image header. However, updates to the image header will 73 also update the table header where appropriate. 74 75 Note that if no image header is passed in, the code will instantiate a 76 regular `~astropy.io.fits.Header`. 77 """ 78 79 # TODO: The difficulty of implementing this screams a need to rewrite this 80 # module 81 82 _keyword_remaps = { 83 'SIMPLE': 'ZSIMPLE', 'XTENSION': 'ZTENSION', 'BITPIX': 'ZBITPIX', 84 'NAXIS': 'ZNAXIS', 'EXTEND': 'ZEXTEND', 'BLOCKED': 'ZBLOCKED', 85 'PCOUNT': 'ZPCOUNT', 'GCOUNT': 'ZGCOUNT', 'CHECKSUM': 'ZHECKSUM', 86 'DATASUM': 'ZDATASUM' 87 } 88 89 _zdef_re = re.compile(r'(?P<label>^[Zz][a-zA-Z]*)(?P<num>[1-9][0-9 ]*$)?') 90 _compression_keywords = set(_keyword_remaps.values()).union( 91 ['ZIMAGE', 'ZCMPTYPE', 'ZMASKCMP', 'ZQUANTIZ', 'ZDITHER0']) 92 _indexed_compression_keywords = {'ZNAXIS', 'ZTILE', 'ZNAME', 'ZVAL'} 93 # TODO: Once it place it should be possible to manage some of this through 94 # the schema system, but it's not quite ready for that yet. Also it still 95 # makes more sense to change CompImageHDU to subclass ImageHDU :/ 96 97 def __new__(cls, table_header, image_header=None): 98 # 2019-09-14 (MHvK): No point wrapping anything if no image_header is 99 # given. This happens if __getitem__ and copy are called - our super 100 # class will aim to initialize a new, possibly partially filled 101 # header, but we cannot usefully deal with that. 102 # TODO: the above suggests strongly we should *not* subclass from 103 # Header. See also comment above about the need for reorganization. 104 if image_header is None: 105 return Header(table_header) 106 else: 107 return super().__new__(cls) 108 109 def __init__(self, table_header, image_header): 110 self._cards = image_header._cards 111 self._keyword_indices = image_header._keyword_indices 112 self._rvkc_indices = image_header._rvkc_indices 113 self._modified = image_header._modified 114 self._table_header = table_header 115 116 # We need to override and Header methods that can modify the header, and 117 # ensure that they sync with the underlying _table_header 118 119 def __setitem__(self, key, value): 120 # This isn't pretty, but if the `key` is either an int or a tuple we 121 # need to figure out what keyword name that maps to before doing 122 # anything else; these checks will be repeated later in the 123 # super().__setitem__ call but I don't see another way around it 124 # without some major refactoring 125 if self._set_slice(key, value, self): 126 return 127 128 if isinstance(key, int): 129 keyword, index = self._keyword_from_index(key) 130 elif isinstance(key, tuple): 131 keyword, index = key 132 else: 133 # We don't want to specify and index otherwise, because that will 134 # break the behavior for new keywords and for commentary keywords 135 keyword, index = key, None 136 137 if self._is_reserved_keyword(keyword): 138 return 139 140 super().__setitem__(key, value) 141 142 if index is not None: 143 remapped_keyword = self._remap_keyword(keyword) 144 self._table_header[remapped_keyword, index] = value 145 # Else this will pass through to ._update 146 147 def __delitem__(self, key): 148 if isinstance(key, slice) or self._haswildcard(key): 149 # If given a slice pass that on to the superclass and bail out 150 # early; we only want to make updates to _table_header when given 151 # a key specifying a single keyword 152 return super().__delitem__(key) 153 154 if isinstance(key, int): 155 keyword, index = self._keyword_from_index(key) 156 elif isinstance(key, tuple): 157 keyword, index = key 158 else: 159 keyword, index = key, None 160 161 if key not in self: 162 raise KeyError(f"Keyword {key!r} not found.") 163 164 super().__delitem__(key) 165 166 remapped_keyword = self._remap_keyword(keyword) 167 168 if remapped_keyword in self._table_header: 169 if index is not None: 170 del self._table_header[(remapped_keyword, index)] 171 else: 172 del self._table_header[remapped_keyword] 173 174 def append(self, card=None, useblanks=True, bottom=False, end=False): 175 # This logic unfortunately needs to be duplicated from the base class 176 # in order to determine the keyword 177 if isinstance(card, str): 178 card = Card(card) 179 elif isinstance(card, tuple): 180 card = Card(*card) 181 elif card is None: 182 card = Card() 183 elif not isinstance(card, Card): 184 raise ValueError( 185 'The value appended to a Header must be either a keyword or ' 186 '(keyword, value, [comment]) tuple; got: {!r}'.format(card)) 187 188 if self._is_reserved_keyword(card.keyword): 189 return 190 191 super().append(card=card, useblanks=useblanks, bottom=bottom, end=end) 192 193 remapped_keyword = self._remap_keyword(card.keyword) 194 195 # card.keyword strips the HIERARCH if present so this must be added 196 # back to avoid a warning. 197 if str(card).startswith("HIERARCH ") and not remapped_keyword.startswith("HIERARCH "): 198 remapped_keyword = "HIERARCH " + remapped_keyword 199 200 card = Card(remapped_keyword, card.value, card.comment) 201 202 # Here we disable the use of blank cards, because the call above to 203 # Header.append may have already deleted a blank card in the table 204 # header, thanks to inheritance: Header.append calls 'del self[-1]' 205 # to delete a blank card, which calls CompImageHeader.__deltitem__, 206 # which deletes the blank card both in the image and the table headers! 207 self._table_header.append(card=card, useblanks=False, 208 bottom=bottom, end=end) 209 210 def insert(self, key, card, useblanks=True, after=False): 211 if isinstance(key, int): 212 # Determine condition to pass through to append 213 if after: 214 if key == -1: 215 key = len(self._cards) 216 else: 217 key += 1 218 219 if key >= len(self._cards): 220 self.append(card, end=True) 221 return 222 223 if isinstance(card, str): 224 card = Card(card) 225 elif isinstance(card, tuple): 226 card = Card(*card) 227 elif not isinstance(card, Card): 228 raise ValueError( 229 'The value inserted into a Header must be either a keyword or ' 230 '(keyword, value, [comment]) tuple; got: {!r}'.format(card)) 231 232 if self._is_reserved_keyword(card.keyword): 233 return 234 235 # Now the tricky part is to determine where to insert in the table 236 # header. If given a numerical index we need to map that to the 237 # corresponding index in the table header. Although rare, there may be 238 # cases where there is no mapping in which case we just try the same 239 # index 240 # NOTE: It is crucial that remapped_index in particular is figured out 241 # before the image header is modified 242 remapped_index = self._remap_index(key) 243 remapped_keyword = self._remap_keyword(card.keyword) 244 245 super().insert(key, card, useblanks=useblanks, after=after) 246 247 card = Card(remapped_keyword, card.value, card.comment) 248 249 # Here we disable the use of blank cards, because the call above to 250 # Header.insert may have already deleted a blank card in the table 251 # header, thanks to inheritance: Header.insert calls 'del self[-1]' 252 # to delete a blank card, which calls CompImageHeader.__delitem__, 253 # which deletes the blank card both in the image and the table headers! 254 self._table_header.insert(remapped_index, card, useblanks=False, 255 after=after) 256 257 def _update(self, card): 258 keyword = card[0] 259 260 if self._is_reserved_keyword(keyword): 261 return 262 263 super()._update(card) 264 265 if keyword in Card._commentary_keywords: 266 # Otherwise this will result in a duplicate insertion 267 return 268 269 remapped_keyword = self._remap_keyword(keyword) 270 self._table_header._update((remapped_keyword,) + card[1:]) 271 272 # Last piece needed (I think) for synchronizing with the real header 273 # This one is tricky since _relativeinsert calls insert 274 def _relativeinsert(self, card, before=None, after=None, replace=False): 275 keyword = card[0] 276 277 if self._is_reserved_keyword(keyword): 278 return 279 280 # Now we have to figure out how to remap 'before' and 'after' 281 if before is None: 282 if isinstance(after, int): 283 remapped_after = self._remap_index(after) 284 else: 285 remapped_after = self._remap_keyword(after) 286 remapped_before = None 287 else: 288 if isinstance(before, int): 289 remapped_before = self._remap_index(before) 290 else: 291 remapped_before = self._remap_keyword(before) 292 remapped_after = None 293 294 super()._relativeinsert(card, before=before, after=after, 295 replace=replace) 296 297 remapped_keyword = self._remap_keyword(keyword) 298 299 card = Card(remapped_keyword, card[1], card[2]) 300 self._table_header._relativeinsert(card, before=remapped_before, 301 after=remapped_after, 302 replace=replace) 303 304 @classmethod 305 def _is_reserved_keyword(cls, keyword, warn=True): 306 msg = ('Keyword {!r} is reserved for use by the FITS Tiled Image ' 307 'Convention and will not be stored in the header for the ' 308 'image being compressed.'.format(keyword)) 309 310 if keyword == 'TFIELDS': 311 if warn: 312 warnings.warn(msg) 313 return True 314 315 m = TDEF_RE.match(keyword) 316 317 if m and m.group('label').upper() in TABLE_KEYWORD_NAMES: 318 if warn: 319 warnings.warn(msg) 320 return True 321 322 m = cls._zdef_re.match(keyword) 323 324 if m: 325 label = m.group('label').upper() 326 num = m.group('num') 327 if num is not None and label in cls._indexed_compression_keywords: 328 if warn: 329 warnings.warn(msg) 330 return True 331 elif label in cls._compression_keywords: 332 if warn: 333 warnings.warn(msg) 334 return True 335 336 return False 337 338 @classmethod 339 def _remap_keyword(cls, keyword): 340 # Given a keyword that one might set on an image, remap that keyword to 341 # the name used for it in the COMPRESSED HDU header 342 # This is mostly just a lookup in _keyword_remaps, but needs handling 343 # for NAXISn keywords 344 345 is_naxisn = False 346 if keyword[:5] == 'NAXIS': 347 with suppress(ValueError): 348 index = int(keyword[5:]) 349 is_naxisn = index > 0 350 351 if is_naxisn: 352 return f'ZNAXIS{index}' 353 354 # If the keyword does not need to be remapped then just return the 355 # original keyword 356 return cls._keyword_remaps.get(keyword, keyword) 357 358 def _remap_index(self, idx): 359 # Given an integer index into this header, map that to the index in the 360 # table header for the same card. If the card doesn't exist in the 361 # table header (generally should *not* be the case) this will just 362 # return the same index 363 # This *does* also accept a keyword or (keyword, repeat) tuple and 364 # obtains the associated numerical index with self._cardindex 365 if not isinstance(idx, int): 366 idx = self._cardindex(idx) 367 368 keyword, repeat = self._keyword_from_index(idx) 369 remapped_insert_keyword = self._remap_keyword(keyword) 370 371 with suppress(IndexError, KeyError): 372 idx = self._table_header._cardindex((remapped_insert_keyword, 373 repeat)) 374 375 return idx 376 377 378# TODO: Fix this class so that it doesn't actually inherit from BinTableHDU, 379# but instead has an internal BinTableHDU reference 380class CompImageHDU(BinTableHDU): 381 """ 382 Compressed Image HDU class. 383 """ 384 385 _manages_own_heap = True 386 """ 387 The calls to CFITSIO lay out the heap data in memory, and we write it out 388 the same way CFITSIO organizes it. In principle this would break if a user 389 manually changes the underlying compressed data by hand, but there is no 390 reason they would want to do that (and if they do that's their 391 responsibility). 392 """ 393 394 _default_name = "COMPRESSED_IMAGE" 395 396 def __init__(self, data=None, header=None, name=None, 397 compression_type=DEFAULT_COMPRESSION_TYPE, 398 tile_size=None, 399 hcomp_scale=DEFAULT_HCOMP_SCALE, 400 hcomp_smooth=DEFAULT_HCOMP_SMOOTH, 401 quantize_level=DEFAULT_QUANTIZE_LEVEL, 402 quantize_method=DEFAULT_QUANTIZE_METHOD, 403 dither_seed=DEFAULT_DITHER_SEED, 404 do_not_scale_image_data=False, 405 uint=False, scale_back=False, **kwargs): 406 """ 407 Parameters 408 ---------- 409 data : array, optional 410 Uncompressed image data 411 412 header : `~astropy.io.fits.Header`, optional 413 Header to be associated with the image; when reading the HDU from a 414 file (data=DELAYED), the header read from the file 415 416 name : str, optional 417 The ``EXTNAME`` value; if this value is `None`, then the name from 418 the input image header will be used; if there is no name in the 419 input image header then the default name ``COMPRESSED_IMAGE`` is 420 used. 421 422 compression_type : str, optional 423 Compression algorithm: one of 424 ``'RICE_1'``, ``'RICE_ONE'``, ``'PLIO_1'``, ``'GZIP_1'``, 425 ``'GZIP_2'``, ``'HCOMPRESS_1'`` 426 427 tile_size : int, optional 428 Compression tile sizes. Default treats each row of image as a 429 tile. 430 431 hcomp_scale : float, optional 432 HCOMPRESS scale parameter 433 434 hcomp_smooth : float, optional 435 HCOMPRESS smooth parameter 436 437 quantize_level : float, optional 438 Floating point quantization level; see note below 439 440 quantize_method : int, optional 441 Floating point quantization dithering method; can be either 442 ``NO_DITHER`` (-1; default), ``SUBTRACTIVE_DITHER_1`` (1), or 443 ``SUBTRACTIVE_DITHER_2`` (2); see note below 444 445 dither_seed : int, optional 446 Random seed to use for dithering; can be either an integer in the 447 range 1 to 1000 (inclusive), ``DITHER_SEED_CLOCK`` (0; default), or 448 ``DITHER_SEED_CHECKSUM`` (-1); see note below 449 450 Notes 451 ----- 452 The astropy.io.fits package supports 2 methods of image compression: 453 454 1) The entire FITS file may be externally compressed with the gzip 455 or pkzip utility programs, producing a ``*.gz`` or ``*.zip`` 456 file, respectively. When reading compressed files of this type, 457 Astropy first uncompresses the entire file into a temporary file 458 before performing the requested read operations. The 459 astropy.io.fits package does not support writing to these types 460 of compressed files. This type of compression is supported in 461 the ``_File`` class, not in the `CompImageHDU` class. The file 462 compression type is recognized by the ``.gz`` or ``.zip`` file 463 name extension. 464 465 2) The `CompImageHDU` class supports the FITS tiled image 466 compression convention in which the image is subdivided into a 467 grid of rectangular tiles, and each tile of pixels is 468 individually compressed. The details of this FITS compression 469 convention are described at the `FITS Support Office web site 470 <https://fits.gsfc.nasa.gov/registry/tilecompression.html>`_. 471 Basically, the compressed image tiles are stored in rows of a 472 variable length array column in a FITS binary table. The 473 astropy.io.fits recognizes that this binary table extension 474 contains an image and treats it as if it were an image 475 extension. Under this tile-compression format, FITS header 476 keywords remain uncompressed. At this time, Astropy does not 477 support the ability to extract and uncompress sections of the 478 image without having to uncompress the entire image. 479 480 The astropy.io.fits package supports 3 general-purpose compression 481 algorithms plus one other special-purpose compression technique that is 482 designed for data masks with positive integer pixel values. The 3 483 general purpose algorithms are GZIP, Rice, and HCOMPRESS, and the 484 special-purpose technique is the IRAF pixel list compression technique 485 (PLIO). The ``compression_type`` parameter defines the compression 486 algorithm to be used. 487 488 The FITS image can be subdivided into any desired rectangular grid of 489 compression tiles. With the GZIP, Rice, and PLIO algorithms, the 490 default is to take each row of the image as a tile. The HCOMPRESS 491 algorithm is inherently 2-dimensional in nature, so the default in this 492 case is to take 16 rows of the image per tile. In most cases, it makes 493 little difference what tiling pattern is used, so the default tiles are 494 usually adequate. In the case of very small images, it could be more 495 efficient to compress the whole image as a single tile. Note that the 496 image dimensions are not required to be an integer multiple of the tile 497 dimensions; if not, then the tiles at the edges of the image will be 498 smaller than the other tiles. The ``tile_size`` parameter may be 499 provided as a list of tile sizes, one for each dimension in the image. 500 For example a ``tile_size`` value of ``[100,100]`` would divide a 300 X 501 300 image into 9 100 X 100 tiles. 502 503 The 4 supported image compression algorithms are all 'lossless' when 504 applied to integer FITS images; the pixel values are preserved exactly 505 with no loss of information during the compression and uncompression 506 process. In addition, the HCOMPRESS algorithm supports a 'lossy' 507 compression mode that will produce larger amount of image compression. 508 This is achieved by specifying a non-zero value for the ``hcomp_scale`` 509 parameter. Since the amount of compression that is achieved depends 510 directly on the RMS noise in the image, it is usually more convenient 511 to specify the ``hcomp_scale`` factor relative to the RMS noise. 512 Setting ``hcomp_scale = 2.5`` means use a scale factor that is 2.5 513 times the calculated RMS noise in the image tile. In some cases it may 514 be desirable to specify the exact scaling to be used, instead of 515 specifying it relative to the calculated noise value. This may be done 516 by specifying the negative of the desired scale value (typically in the 517 range -2 to -100). 518 519 Very high compression factors (of 100 or more) can be achieved by using 520 large ``hcomp_scale`` values, however, this can produce undesirable 521 'blocky' artifacts in the compressed image. A variation of the 522 HCOMPRESS algorithm (called HSCOMPRESS) can be used in this case to 523 apply a small amount of smoothing of the image when it is uncompressed 524 to help cover up these artifacts. This smoothing is purely cosmetic 525 and does not cause any significant change to the image pixel values. 526 Setting the ``hcomp_smooth`` parameter to 1 will engage the smoothing 527 algorithm. 528 529 Floating point FITS images (which have ``BITPIX`` = -32 or -64) usually 530 contain too much 'noise' in the least significant bits of the mantissa 531 of the pixel values to be effectively compressed with any lossless 532 algorithm. Consequently, floating point images are first quantized 533 into scaled integer pixel values (and thus throwing away much of the 534 noise) before being compressed with the specified algorithm (either 535 GZIP, RICE, or HCOMPRESS). This technique produces much higher 536 compression factors than simply using the GZIP utility to externally 537 compress the whole FITS file, but it also means that the original 538 floating point value pixel values are not exactly preserved. When done 539 properly, this integer scaling technique will only discard the 540 insignificant noise while still preserving all the real information in 541 the image. The amount of precision that is retained in the pixel 542 values is controlled by the ``quantize_level`` parameter. Larger 543 values will result in compressed images whose pixels more closely match 544 the floating point pixel values, but at the same time the amount of 545 compression that is achieved will be reduced. Users should experiment 546 with different values for this parameter to determine the optimal value 547 that preserves all the useful information in the image, without 548 needlessly preserving all the 'noise' which will hurt the compression 549 efficiency. 550 551 The default value for the ``quantize_level`` scale factor is 16, which 552 means that scaled integer pixel values will be quantized such that the 553 difference between adjacent integer values will be 1/16th of the noise 554 level in the image background. An optimized algorithm is used to 555 accurately estimate the noise in the image. As an example, if the RMS 556 noise in the background pixels of an image = 32.0, then the spacing 557 between adjacent scaled integer pixel values will equal 2.0 by default. 558 Note that the RMS noise is independently calculated for each tile of 559 the image, so the resulting integer scaling factor may fluctuate 560 slightly for each tile. In some cases, it may be desirable to specify 561 the exact quantization level to be used, instead of specifying it 562 relative to the calculated noise value. This may be done by specifying 563 the negative of desired quantization level for the value of 564 ``quantize_level``. In the previous example, one could specify 565 ``quantize_level = -2.0`` so that the quantized integer levels differ 566 by 2.0. Larger negative values for ``quantize_level`` means that the 567 levels are more coarsely-spaced, and will produce higher compression 568 factors. 569 570 The quantization algorithm can also apply one of two random dithering 571 methods in order to reduce bias in the measured intensity of background 572 regions. The default method, specified with the constant 573 ``SUBTRACTIVE_DITHER_1`` adds dithering to the zero-point of the 574 quantization array itself rather than adding noise to the actual image. 575 The random noise is added on a pixel-by-pixel basis, so in order 576 restore each pixel from its integer value to its floating point value 577 it is necessary to replay the same sequence of random numbers for each 578 pixel (see below). The other method, ``SUBTRACTIVE_DITHER_2``, is 579 exactly like the first except that before dithering any pixel with a 580 floating point value of ``0.0`` is replaced with the special integer 581 value ``-2147483647``. When the image is uncompressed, pixels with 582 this value are restored back to ``0.0`` exactly. Finally, a value of 583 ``NO_DITHER`` disables dithering entirely. 584 585 As mentioned above, when using the subtractive dithering algorithm it 586 is necessary to be able to generate a (pseudo-)random sequence of noise 587 for each pixel, and replay that same sequence upon decompressing. To 588 facilitate this, a random seed between 1 and 10000 (inclusive) is used 589 to seed a random number generator, and that seed is stored in the 590 ``ZDITHER0`` keyword in the header of the compressed HDU. In order to 591 use that seed to generate the same sequence of random numbers the same 592 random number generator must be used at compression and decompression 593 time; for that reason the tiled image convention provides an 594 implementation of a very simple pseudo-random number generator. The 595 seed itself can be provided in one of three ways, controllable by the 596 ``dither_seed`` argument: It may be specified manually, or it may be 597 generated arbitrarily based on the system's clock 598 (``DITHER_SEED_CLOCK``) or based on a checksum of the pixels in the 599 image's first tile (``DITHER_SEED_CHECKSUM``). The clock-based method 600 is the default, and is sufficient to ensure that the value is 601 reasonably "arbitrary" and that the same seed is unlikely to be 602 generated sequentially. The checksum method, on the other hand, 603 ensures that the same seed is used every time for a specific image. 604 This is particularly useful for software testing as it ensures that the 605 same image will always use the same seed. 606 """ 607 608 if not COMPRESSION_SUPPORTED: 609 # TODO: Raise a more specific Exception type 610 raise Exception('The astropy.io.fits.compression module is not ' 611 'available. Creation of compressed image HDUs is ' 612 'disabled.') 613 614 compression_type = CMTYPE_ALIASES.get(compression_type, compression_type) 615 616 if data is DELAYED: 617 # Reading the HDU from a file 618 super().__init__(data=data, header=header) 619 else: 620 # Create at least a skeleton HDU that matches the input 621 # header and data (if any were input) 622 super().__init__(data=None, header=header) 623 624 # Store the input image data 625 self.data = data 626 627 # Update the table header (_header) to the compressed 628 # image format and to match the input data (if any); 629 # Create the image header (_image_header) from the input 630 # image header (if any) and ensure it matches the input 631 # data; Create the initially empty table data array to 632 # hold the compressed data. 633 self._update_header_data(header, name, 634 compression_type=compression_type, 635 tile_size=tile_size, 636 hcomp_scale=hcomp_scale, 637 hcomp_smooth=hcomp_smooth, 638 quantize_level=quantize_level, 639 quantize_method=quantize_method, 640 dither_seed=dither_seed) 641 642 # TODO: A lot of this should be passed on to an internal image HDU o 643 # something like that, see ticket #88 644 self._do_not_scale_image_data = do_not_scale_image_data 645 self._uint = uint 646 self._scale_back = scale_back 647 648 self._axes = [self._header.get('ZNAXIS' + str(axis + 1), 0) 649 for axis in range(self._header.get('ZNAXIS', 0))] 650 651 # store any scale factors from the table header 652 if do_not_scale_image_data: 653 self._bzero = 0 654 self._bscale = 1 655 else: 656 self._bzero = self._header.get('BZERO', 0) 657 self._bscale = self._header.get('BSCALE', 1) 658 self._bitpix = self._header['ZBITPIX'] 659 660 self._orig_bzero = self._bzero 661 self._orig_bscale = self._bscale 662 self._orig_bitpix = self._bitpix 663 664 def _remove_unnecessary_default_extnames(self, header): 665 """Remove default EXTNAME values if they are unnecessary. 666 667 Some data files (eg from CFHT) can have the default EXTNAME and 668 an explicit value. This method removes the default if a more 669 specific header exists. It also removes any duplicate default 670 values. 671 """ 672 if 'EXTNAME' in header: 673 indices = header._keyword_indices['EXTNAME'] 674 # Only continue if there is more than one found 675 n_extname = len(indices) 676 if n_extname > 1: 677 extnames_to_remove = [index for index in indices 678 if header[index] == self._default_name] 679 if len(extnames_to_remove) == n_extname: 680 # Keep the first (they are all the same) 681 extnames_to_remove.pop(0) 682 # Remove them all in reverse order to keep the index unchanged. 683 for index in reversed(sorted(extnames_to_remove)): 684 del header[index] 685 686 @property 687 def name(self): 688 # Convert the value to a string to be flexible in some pathological 689 # cases (see ticket #96) 690 # Similar to base class but uses .header rather than ._header 691 return str(self.header.get('EXTNAME', self._default_name)) 692 693 @name.setter 694 def name(self, value): 695 # This is a copy of the base class but using .header instead 696 # of ._header to ensure that the name stays in sync. 697 if not isinstance(value, str): 698 raise TypeError("'name' attribute must be a string") 699 if not conf.extension_name_case_sensitive: 700 value = value.upper() 701 if 'EXTNAME' in self.header: 702 self.header['EXTNAME'] = value 703 else: 704 self.header['EXTNAME'] = (value, 'extension name') 705 706 @classmethod 707 def match_header(cls, header): 708 card = header.cards[0] 709 if card.keyword != 'XTENSION': 710 return False 711 712 xtension = card.value 713 if isinstance(xtension, str): 714 xtension = xtension.rstrip() 715 716 if xtension not in ('BINTABLE', 'A3DTABLE'): 717 return False 718 719 if 'ZIMAGE' not in header or not header['ZIMAGE']: 720 return False 721 722 if COMPRESSION_SUPPORTED and COMPRESSION_ENABLED: 723 return True 724 elif not COMPRESSION_SUPPORTED: 725 warnings.warn('Failure matching header to a compressed image ' 726 'HDU: The compression module is not available.\n' 727 'The HDU will be treated as a Binary Table HDU.', 728 AstropyUserWarning) 729 return False 730 else: 731 # Compression is supported but disabled; just pass silently (#92) 732 return False 733 734 def _update_header_data(self, image_header, 735 name=None, 736 compression_type=None, 737 tile_size=None, 738 hcomp_scale=None, 739 hcomp_smooth=None, 740 quantize_level=None, 741 quantize_method=None, 742 dither_seed=None): 743 """ 744 Update the table header (`_header`) to the compressed 745 image format and to match the input data (if any). Create 746 the image header (`_image_header`) from the input image 747 header (if any) and ensure it matches the input 748 data. Create the initially-empty table data array to hold 749 the compressed data. 750 751 This method is mainly called internally, but a user may wish to 752 call this method after assigning new data to the `CompImageHDU` 753 object that is of a different type. 754 755 Parameters 756 ---------- 757 image_header : `~astropy.io.fits.Header` 758 header to be associated with the image 759 760 name : str, optional 761 the ``EXTNAME`` value; if this value is `None`, then the name from 762 the input image header will be used; if there is no name in the 763 input image header then the default name 'COMPRESSED_IMAGE' is used 764 765 compression_type : str, optional 766 compression algorithm 'RICE_1', 'PLIO_1', 'GZIP_1', 'GZIP_2', 767 'HCOMPRESS_1'; if this value is `None`, use value already in the 768 header; if no value already in the header, use 'RICE_1' 769 770 tile_size : sequence of int, optional 771 compression tile sizes as a list; if this value is `None`, use 772 value already in the header; if no value already in the header, 773 treat each row of image as a tile 774 775 hcomp_scale : float, optional 776 HCOMPRESS scale parameter; if this value is `None`, use the value 777 already in the header; if no value already in the header, use 1 778 779 hcomp_smooth : float, optional 780 HCOMPRESS smooth parameter; if this value is `None`, use the value 781 already in the header; if no value already in the header, use 0 782 783 quantize_level : float, optional 784 floating point quantization level; if this value is `None`, use the 785 value already in the header; if no value already in header, use 16 786 787 quantize_method : int, optional 788 floating point quantization dithering method; can be either 789 NO_DITHER (-1), SUBTRACTIVE_DITHER_1 (1; default), or 790 SUBTRACTIVE_DITHER_2 (2) 791 792 dither_seed : int, optional 793 random seed to use for dithering; can be either an integer in the 794 range 1 to 1000 (inclusive), DITHER_SEED_CLOCK (0; default), or 795 DITHER_SEED_CHECKSUM (-1) 796 """ 797 798 # Clean up EXTNAME duplicates 799 self._remove_unnecessary_default_extnames(self._header) 800 801 image_hdu = ImageHDU(data=self.data, header=self._header) 802 self._image_header = CompImageHeader(self._header, image_hdu.header) 803 self._axes = image_hdu._axes 804 del image_hdu 805 806 # Determine based on the size of the input data whether to use the Q 807 # column format to store compressed data or the P format. 808 # The Q format is used only if the uncompressed data is larger than 809 # 4 GB. This is not a perfect heuristic, as one can contrive an input 810 # array which, when compressed, the entire binary table representing 811 # the compressed data is larger than 4GB. That said, this is the same 812 # heuristic used by CFITSIO, so this should give consistent results. 813 # And the cases where this heuristic is insufficient are extreme and 814 # almost entirely contrived corner cases, so it will do for now 815 if self._has_data: 816 huge_hdu = self.data.nbytes > 2 ** 32 817 else: 818 huge_hdu = False 819 820 # Update the extension name in the table header 821 if not name and 'EXTNAME' not in self._header: 822 # Do not sync this with the image header since the default 823 # name is specific to the table header. 824 self._header.set('EXTNAME', self._default_name, 825 'name of this binary table extension', 826 after='TFIELDS') 827 elif name: 828 # Force the name into table and image headers. 829 self.name = name 830 831 # Set the compression type in the table header. 832 if compression_type: 833 if compression_type not in COMPRESSION_TYPES: 834 warnings.warn( 835 'Unknown compression type provided (supported are {}). ' 836 'Default ({}) compression will be used.' 837 .format(', '.join(map(repr, COMPRESSION_TYPES)), 838 DEFAULT_COMPRESSION_TYPE), 839 AstropyUserWarning) 840 compression_type = DEFAULT_COMPRESSION_TYPE 841 842 self._header.set('ZCMPTYPE', compression_type, 843 'compression algorithm', after='TFIELDS') 844 else: 845 compression_type = self._header.get('ZCMPTYPE', 846 DEFAULT_COMPRESSION_TYPE) 847 compression_type = CMTYPE_ALIASES.get(compression_type, 848 compression_type) 849 850 # If the input image header had BSCALE/BZERO cards, then insert 851 # them in the table header. 852 853 if image_header: 854 bzero = image_header.get('BZERO', 0.0) 855 bscale = image_header.get('BSCALE', 1.0) 856 after_keyword = 'EXTNAME' 857 858 if bscale != 1.0: 859 self._header.set('BSCALE', bscale, after=after_keyword) 860 after_keyword = 'BSCALE' 861 862 if bzero != 0.0: 863 self._header.set('BZERO', bzero, after=after_keyword) 864 865 try: 866 bitpix_comment = image_header.comments['BITPIX'] 867 except (AttributeError, KeyError): 868 bitpix_comment = 'data type of original image' 869 870 try: 871 naxis_comment = image_header.comments['NAXIS'] 872 except (AttributeError, KeyError): 873 naxis_comment = 'dimension of original image' 874 875 # Set the label for the first column in the table 876 877 self._header.set('TTYPE1', 'COMPRESSED_DATA', 'label for field 1', 878 after='TFIELDS') 879 880 # Set the data format for the first column. It is dependent 881 # on the requested compression type. 882 883 if compression_type == 'PLIO_1': 884 tform1 = '1QI' if huge_hdu else '1PI' 885 else: 886 tform1 = '1QB' if huge_hdu else '1PB' 887 888 self._header.set('TFORM1', tform1, 889 'data format of field: variable length array', 890 after='TTYPE1') 891 892 # Create the first column for the table. This column holds the 893 # compressed data. 894 col1 = Column(name=self._header['TTYPE1'], format=tform1) 895 896 # Create the additional columns required for floating point 897 # data and calculate the width of the output table. 898 899 zbitpix = self._image_header['BITPIX'] 900 901 if zbitpix < 0 and quantize_level != 0.0: 902 # floating point image has 'COMPRESSED_DATA', 903 # 'UNCOMPRESSED_DATA', 'ZSCALE', and 'ZZERO' columns (unless using 904 # lossless compression, per CFITSIO) 905 ncols = 4 906 907 # CFITSIO 3.28 and up automatically use the GZIP_COMPRESSED_DATA 908 # store floating point data that couldn't be quantized, instead 909 # of the UNCOMPRESSED_DATA column. There's no way to control 910 # this behavior so the only way to determine which behavior will 911 # be employed is via the CFITSIO version 912 913 ttype2 = 'GZIP_COMPRESSED_DATA' 914 # The required format for the GZIP_COMPRESSED_DATA is actually 915 # missing from the standard docs, but CFITSIO suggests it 916 # should be 1PB, which is logical. 917 tform2 = '1QB' if huge_hdu else '1PB' 918 919 # Set up the second column for the table that will hold any 920 # uncompressable data. 921 self._header.set('TTYPE2', ttype2, 'label for field 2', 922 after='TFORM1') 923 924 self._header.set('TFORM2', tform2, 925 'data format of field: variable length array', 926 after='TTYPE2') 927 928 col2 = Column(name=ttype2, format=tform2) 929 930 # Set up the third column for the table that will hold 931 # the scale values for quantized data. 932 self._header.set('TTYPE3', 'ZSCALE', 'label for field 3', 933 after='TFORM2') 934 self._header.set('TFORM3', '1D', 935 'data format of field: 8-byte DOUBLE', 936 after='TTYPE3') 937 col3 = Column(name=self._header['TTYPE3'], 938 format=self._header['TFORM3']) 939 940 # Set up the fourth column for the table that will hold 941 # the zero values for the quantized data. 942 self._header.set('TTYPE4', 'ZZERO', 'label for field 4', 943 after='TFORM3') 944 self._header.set('TFORM4', '1D', 945 'data format of field: 8-byte DOUBLE', 946 after='TTYPE4') 947 after = 'TFORM4' 948 col4 = Column(name=self._header['TTYPE4'], 949 format=self._header['TFORM4']) 950 951 # Create the ColDefs object for the table 952 cols = ColDefs([col1, col2, col3, col4]) 953 else: 954 # default table has just one 'COMPRESSED_DATA' column 955 ncols = 1 956 after = 'TFORM1' 957 958 # remove any header cards for the additional columns that 959 # may be left over from the previous data 960 to_remove = ['TTYPE2', 'TFORM2', 'TTYPE3', 'TFORM3', 'TTYPE4', 961 'TFORM4'] 962 963 for k in to_remove: 964 try: 965 del self._header[k] 966 except KeyError: 967 pass 968 969 # Create the ColDefs object for the table 970 cols = ColDefs([col1]) 971 972 # Update the table header with the width of the table, the 973 # number of fields in the table, the indicator for a compressed 974 # image HDU, the data type of the image data and the number of 975 # dimensions in the image data array. 976 self._header.set('NAXIS1', cols.dtype.itemsize, 977 'width of table in bytes') 978 self._header.set('TFIELDS', ncols, 'number of fields in each row', 979 after='GCOUNT') 980 self._header.set('ZIMAGE', True, 'extension contains compressed image', 981 after=after) 982 self._header.set('ZBITPIX', zbitpix, 983 bitpix_comment, after='ZIMAGE') 984 self._header.set('ZNAXIS', self._image_header['NAXIS'], naxis_comment, 985 after='ZBITPIX') 986 987 # Strip the table header of all the ZNAZISn and ZTILEn keywords 988 # that may be left over from the previous data 989 990 for idx in itertools.count(1): 991 try: 992 del self._header['ZNAXIS' + str(idx)] 993 del self._header['ZTILE' + str(idx)] 994 except KeyError: 995 break 996 997 # Verify that any input tile size parameter is the appropriate 998 # size to match the HDU's data. 999 1000 naxis = self._image_header['NAXIS'] 1001 1002 if not tile_size: 1003 tile_size = [] 1004 elif len(tile_size) != naxis: 1005 warnings.warn('Provided tile size not appropriate for the data. ' 1006 'Default tile size will be used.', AstropyUserWarning) 1007 tile_size = [] 1008 1009 # Set default tile dimensions for HCOMPRESS_1 1010 1011 if compression_type == 'HCOMPRESS_1': 1012 if (self._image_header['NAXIS1'] < 4 or 1013 self._image_header['NAXIS2'] < 4): 1014 raise ValueError('Hcompress minimum image dimension is ' 1015 '4 pixels') 1016 elif tile_size: 1017 if tile_size[0] < 4 or tile_size[1] < 4: 1018 # user specified tile size is too small 1019 raise ValueError('Hcompress minimum tile dimension is ' 1020 '4 pixels') 1021 major_dims = len([ts for ts in tile_size if ts > 1]) 1022 if major_dims > 2: 1023 raise ValueError( 1024 'HCOMPRESS can only support 2-dimensional tile sizes.' 1025 'All but two of the tile_size dimensions must be set ' 1026 'to 1.') 1027 1028 if tile_size and (tile_size[0] == 0 and tile_size[1] == 0): 1029 # compress the whole image as a single tile 1030 tile_size[0] = self._image_header['NAXIS1'] 1031 tile_size[1] = self._image_header['NAXIS2'] 1032 1033 for i in range(2, naxis): 1034 # set all higher tile dimensions = 1 1035 tile_size[i] = 1 1036 elif not tile_size: 1037 # The Hcompress algorithm is inherently 2D in nature, so the 1038 # row by row tiling that is used for other compression 1039 # algorithms is not appropriate. If the image has less than 30 1040 # rows, then the entire image will be compressed as a single 1041 # tile. Otherwise the tiles will consist of 16 rows of the 1042 # image. This keeps the tiles to a reasonable size, and it 1043 # also includes enough rows to allow good compression 1044 # efficiency. It the last tile of the image happens to contain 1045 # less than 4 rows, then find another tile size with between 14 1046 # and 30 rows (preferably even), so that the last tile has at 1047 # least 4 rows. 1048 1049 # 1st tile dimension is the row length of the image 1050 tile_size.append(self._image_header['NAXIS1']) 1051 1052 if self._image_header['NAXIS2'] <= 30: 1053 tile_size.append(self._image_header['NAXIS1']) 1054 else: 1055 # look for another good tile dimension 1056 naxis2 = self._image_header['NAXIS2'] 1057 for dim in [16, 24, 20, 30, 28, 26, 22, 18, 14]: 1058 if naxis2 % dim == 0 or naxis2 % dim > 3: 1059 tile_size.append(dim) 1060 break 1061 else: 1062 tile_size.append(17) 1063 1064 for i in range(2, naxis): 1065 # set all higher tile dimensions = 1 1066 tile_size.append(1) 1067 1068 # check if requested tile size causes the last tile to have 1069 # less than 4 pixels 1070 1071 remain = self._image_header['NAXIS1'] % tile_size[0] # 1st dimen 1072 1073 if remain > 0 and remain < 4: 1074 tile_size[0] += 1 # try increasing tile size by 1 1075 1076 remain = self._image_header['NAXIS1'] % tile_size[0] 1077 1078 if remain > 0 and remain < 4: 1079 raise ValueError('Last tile along 1st dimension has ' 1080 'less than 4 pixels') 1081 1082 remain = self._image_header['NAXIS2'] % tile_size[1] # 2nd dimen 1083 1084 if remain > 0 and remain < 4: 1085 tile_size[1] += 1 # try increasing tile size by 1 1086 1087 remain = self._image_header['NAXIS2'] % tile_size[1] 1088 1089 if remain > 0 and remain < 4: 1090 raise ValueError('Last tile along 2nd dimension has ' 1091 'less than 4 pixels') 1092 1093 # Set up locations for writing the next cards in the header. 1094 last_znaxis = 'ZNAXIS' 1095 1096 if self._image_header['NAXIS'] > 0: 1097 after1 = 'ZNAXIS1' 1098 else: 1099 after1 = 'ZNAXIS' 1100 1101 # Calculate the number of rows in the output table and 1102 # write the ZNAXISn and ZTILEn cards to the table header. 1103 nrows = 0 1104 1105 for idx, axis in enumerate(self._axes): 1106 naxis = 'NAXIS' + str(idx + 1) 1107 znaxis = 'ZNAXIS' + str(idx + 1) 1108 ztile = 'ZTILE' + str(idx + 1) 1109 1110 if tile_size and len(tile_size) >= idx + 1: 1111 ts = tile_size[idx] 1112 else: 1113 if ztile not in self._header: 1114 # Default tile size 1115 if not idx: 1116 ts = self._image_header['NAXIS1'] 1117 else: 1118 ts = 1 1119 else: 1120 ts = self._header[ztile] 1121 tile_size.append(ts) 1122 1123 if not nrows: 1124 nrows = (axis - 1) // ts + 1 1125 else: 1126 nrows *= ((axis - 1) // ts + 1) 1127 1128 if image_header and naxis in image_header: 1129 self._header.set(znaxis, axis, image_header.comments[naxis], 1130 after=last_znaxis) 1131 else: 1132 self._header.set(znaxis, axis, 1133 'length of original image axis', 1134 after=last_znaxis) 1135 1136 self._header.set(ztile, ts, 'size of tiles to be compressed', 1137 after=after1) 1138 last_znaxis = znaxis 1139 after1 = ztile 1140 1141 # Set the NAXIS2 header card in the table hdu to the number of 1142 # rows in the table. 1143 self._header.set('NAXIS2', nrows, 'number of rows in table') 1144 1145 self.columns = cols 1146 1147 # Set the compression parameters in the table header. 1148 1149 # First, setup the values to be used for the compression parameters 1150 # in case none were passed in. This will be either the value 1151 # already in the table header for that parameter or the default 1152 # value. 1153 for idx in itertools.count(1): 1154 zname = 'ZNAME' + str(idx) 1155 if zname not in self._header: 1156 break 1157 zval = 'ZVAL' + str(idx) 1158 if self._header[zname] == 'NOISEBIT': 1159 if quantize_level is None: 1160 quantize_level = self._header[zval] 1161 if self._header[zname] == 'SCALE ': 1162 if hcomp_scale is None: 1163 hcomp_scale = self._header[zval] 1164 if self._header[zname] == 'SMOOTH ': 1165 if hcomp_smooth is None: 1166 hcomp_smooth = self._header[zval] 1167 1168 if quantize_level is None: 1169 quantize_level = DEFAULT_QUANTIZE_LEVEL 1170 1171 if hcomp_scale is None: 1172 hcomp_scale = DEFAULT_HCOMP_SCALE 1173 1174 if hcomp_smooth is None: 1175 hcomp_smooth = DEFAULT_HCOMP_SCALE 1176 1177 # Next, strip the table header of all the ZNAMEn and ZVALn keywords 1178 # that may be left over from the previous data 1179 for idx in itertools.count(1): 1180 zname = 'ZNAME' + str(idx) 1181 if zname not in self._header: 1182 break 1183 zval = 'ZVAL' + str(idx) 1184 del self._header[zname] 1185 del self._header[zval] 1186 1187 # Finally, put the appropriate keywords back based on the 1188 # compression type. 1189 1190 after_keyword = 'ZCMPTYPE' 1191 idx = 1 1192 1193 if compression_type == 'RICE_1': 1194 self._header.set('ZNAME1', 'BLOCKSIZE', 'compression block size', 1195 after=after_keyword) 1196 self._header.set('ZVAL1', DEFAULT_BLOCK_SIZE, 'pixels per block', 1197 after='ZNAME1') 1198 1199 self._header.set('ZNAME2', 'BYTEPIX', 1200 'bytes per pixel (1, 2, 4, or 8)', after='ZVAL1') 1201 1202 if self._header['ZBITPIX'] == 8: 1203 bytepix = 1 1204 elif self._header['ZBITPIX'] == 16: 1205 bytepix = 2 1206 else: 1207 bytepix = DEFAULT_BYTE_PIX 1208 1209 self._header.set('ZVAL2', bytepix, 1210 'bytes per pixel (1, 2, 4, or 8)', 1211 after='ZNAME2') 1212 after_keyword = 'ZVAL2' 1213 idx = 3 1214 elif compression_type == 'HCOMPRESS_1': 1215 self._header.set('ZNAME1', 'SCALE', 'HCOMPRESS scale factor', 1216 after=after_keyword) 1217 self._header.set('ZVAL1', hcomp_scale, 'HCOMPRESS scale factor', 1218 after='ZNAME1') 1219 self._header.set('ZNAME2', 'SMOOTH', 'HCOMPRESS smooth option', 1220 after='ZVAL1') 1221 self._header.set('ZVAL2', hcomp_smooth, 'HCOMPRESS smooth option', 1222 after='ZNAME2') 1223 after_keyword = 'ZVAL2' 1224 idx = 3 1225 1226 if self._image_header['BITPIX'] < 0: # floating point image 1227 self._header.set('ZNAME' + str(idx), 'NOISEBIT', 1228 'floating point quantization level', 1229 after=after_keyword) 1230 self._header.set('ZVAL' + str(idx), quantize_level, 1231 'floating point quantization level', 1232 after='ZNAME' + str(idx)) 1233 1234 # Add the dither method and seed 1235 if quantize_method: 1236 if quantize_method not in [NO_DITHER, SUBTRACTIVE_DITHER_1, 1237 SUBTRACTIVE_DITHER_2]: 1238 name = QUANTIZE_METHOD_NAMES[DEFAULT_QUANTIZE_METHOD] 1239 warnings.warn('Unknown quantization method provided. ' 1240 'Default method ({}) used.'.format(name)) 1241 quantize_method = DEFAULT_QUANTIZE_METHOD 1242 1243 if quantize_method == NO_DITHER: 1244 zquantiz_comment = 'No dithering during quantization' 1245 else: 1246 zquantiz_comment = 'Pixel Quantization Algorithm' 1247 1248 self._header.set('ZQUANTIZ', 1249 QUANTIZE_METHOD_NAMES[quantize_method], 1250 zquantiz_comment, 1251 after='ZVAL' + str(idx)) 1252 else: 1253 # If the ZQUANTIZ keyword is missing the default is to assume 1254 # no dithering, rather than whatever DEFAULT_QUANTIZE_METHOD 1255 # is set to 1256 quantize_method = self._header.get('ZQUANTIZ', NO_DITHER) 1257 1258 if isinstance(quantize_method, str): 1259 for k, v in QUANTIZE_METHOD_NAMES.items(): 1260 if v.upper() == quantize_method: 1261 quantize_method = k 1262 break 1263 else: 1264 quantize_method = NO_DITHER 1265 1266 if quantize_method == NO_DITHER: 1267 if 'ZDITHER0' in self._header: 1268 # If dithering isn't being used then there's no reason to 1269 # keep the ZDITHER0 keyword 1270 del self._header['ZDITHER0'] 1271 else: 1272 if dither_seed: 1273 dither_seed = self._generate_dither_seed(dither_seed) 1274 elif 'ZDITHER0' in self._header: 1275 dither_seed = self._header['ZDITHER0'] 1276 else: 1277 dither_seed = self._generate_dither_seed( 1278 DEFAULT_DITHER_SEED) 1279 1280 self._header.set('ZDITHER0', dither_seed, 1281 'dithering offset when quantizing floats', 1282 after='ZQUANTIZ') 1283 1284 if image_header: 1285 # Move SIMPLE card from the image header to the 1286 # table header as ZSIMPLE card. 1287 1288 if 'SIMPLE' in image_header: 1289 self._header.set('ZSIMPLE', image_header['SIMPLE'], 1290 image_header.comments['SIMPLE'], 1291 before='ZBITPIX') 1292 1293 # Move EXTEND card from the image header to the 1294 # table header as ZEXTEND card. 1295 1296 if 'EXTEND' in image_header: 1297 self._header.set('ZEXTEND', image_header['EXTEND'], 1298 image_header.comments['EXTEND']) 1299 1300 # Move BLOCKED card from the image header to the 1301 # table header as ZBLOCKED card. 1302 1303 if 'BLOCKED' in image_header: 1304 self._header.set('ZBLOCKED', image_header['BLOCKED'], 1305 image_header.comments['BLOCKED']) 1306 1307 # Move XTENSION card from the image header to the 1308 # table header as ZTENSION card. 1309 1310 # Since we only handle compressed IMAGEs, ZTENSION should 1311 # always be IMAGE, even if the caller has passed in a header 1312 # for some other type of extension. 1313 if 'XTENSION' in image_header: 1314 self._header.set('ZTENSION', 'IMAGE', 1315 image_header.comments['XTENSION'], 1316 before='ZBITPIX') 1317 1318 # Move PCOUNT and GCOUNT cards from image header to the table 1319 # header as ZPCOUNT and ZGCOUNT cards. 1320 1321 if 'PCOUNT' in image_header: 1322 self._header.set('ZPCOUNT', image_header['PCOUNT'], 1323 image_header.comments['PCOUNT'], 1324 after=last_znaxis) 1325 1326 if 'GCOUNT' in image_header: 1327 self._header.set('ZGCOUNT', image_header['GCOUNT'], 1328 image_header.comments['GCOUNT'], 1329 after='ZPCOUNT') 1330 1331 # Move CHECKSUM and DATASUM cards from the image header to the 1332 # table header as XHECKSUM and XDATASUM cards. 1333 1334 if 'CHECKSUM' in image_header: 1335 self._header.set('ZHECKSUM', image_header['CHECKSUM'], 1336 image_header.comments['CHECKSUM']) 1337 1338 if 'DATASUM' in image_header: 1339 self._header.set('ZDATASUM', image_header['DATASUM'], 1340 image_header.comments['DATASUM']) 1341 else: 1342 # Move XTENSION card from the image header to the 1343 # table header as ZTENSION card. 1344 1345 # Since we only handle compressed IMAGEs, ZTENSION should 1346 # always be IMAGE, even if the caller has passed in a header 1347 # for some other type of extension. 1348 if 'XTENSION' in self._image_header: 1349 self._header.set('ZTENSION', 'IMAGE', 1350 self._image_header.comments['XTENSION'], 1351 before='ZBITPIX') 1352 1353 # Move PCOUNT and GCOUNT cards from image header to the table 1354 # header as ZPCOUNT and ZGCOUNT cards. 1355 1356 if 'PCOUNT' in self._image_header: 1357 self._header.set('ZPCOUNT', self._image_header['PCOUNT'], 1358 self._image_header.comments['PCOUNT'], 1359 after=last_znaxis) 1360 1361 if 'GCOUNT' in self._image_header: 1362 self._header.set('ZGCOUNT', self._image_header['GCOUNT'], 1363 self._image_header.comments['GCOUNT'], 1364 after='ZPCOUNT') 1365 1366 # When we have an image checksum we need to ensure that the same 1367 # number of blank cards exist in the table header as there were in 1368 # the image header. This allows those blank cards to be carried 1369 # over to the image header when the hdu is uncompressed. 1370 1371 if 'ZHECKSUM' in self._header: 1372 required_blanks = image_header._countblanks() 1373 image_blanks = self._image_header._countblanks() 1374 table_blanks = self._header._countblanks() 1375 1376 for _ in range(required_blanks - image_blanks): 1377 self._image_header.append() 1378 table_blanks += 1 1379 1380 for _ in range(required_blanks - table_blanks): 1381 self._header.append() 1382 1383 @lazyproperty 1384 def data(self): 1385 # The data attribute is the image data (not the table data). 1386 data = compression.decompress_hdu(self) 1387 1388 if data is None: 1389 return data 1390 1391 # Scale the data if necessary 1392 if (self._orig_bzero != 0 or self._orig_bscale != 1): 1393 new_dtype = self._dtype_for_bitpix() 1394 data = np.array(data, dtype=new_dtype) 1395 1396 zblank = None 1397 1398 if 'ZBLANK' in self.compressed_data.columns.names: 1399 zblank = self.compressed_data['ZBLANK'] 1400 else: 1401 if 'ZBLANK' in self._header: 1402 zblank = np.array(self._header['ZBLANK'], dtype='int32') 1403 elif 'BLANK' in self._header: 1404 zblank = np.array(self._header['BLANK'], dtype='int32') 1405 1406 if zblank is not None: 1407 blanks = (data == zblank) 1408 1409 if self._bscale != 1: 1410 np.multiply(data, self._bscale, data) 1411 if self._bzero != 0: 1412 # We have to explicitly cast self._bzero to prevent numpy from 1413 # raising an error when doing self.data += self._bzero, and we 1414 # do this instead of self.data = self.data + self._bzero to 1415 # avoid doubling memory usage. 1416 np.add(data, self._bzero, out=data, casting='unsafe') 1417 1418 if zblank is not None: 1419 data = np.where(blanks, np.nan, data) 1420 1421 # Right out of _ImageBaseHDU.data 1422 self._update_header_scale_info(data.dtype) 1423 1424 return data 1425 1426 @data.setter 1427 def data(self, data): 1428 if (data is not None) and (not isinstance(data, np.ndarray) or 1429 data.dtype.fields is not None): 1430 raise TypeError('CompImageHDU data has incorrect type:{}; ' 1431 'dtype.fields = {}'.format( 1432 type(data), data.dtype.fields)) 1433 1434 @lazyproperty 1435 def compressed_data(self): 1436 # First we will get the table data (the compressed 1437 # data) from the file, if there is any. 1438 compressed_data = super().data 1439 if isinstance(compressed_data, np.rec.recarray): 1440 # Make sure not to use 'del self.data' so we don't accidentally 1441 # go through the self.data.fdel and close the mmap underlying 1442 # the compressed_data array 1443 del self.__dict__['data'] 1444 return compressed_data 1445 else: 1446 # This will actually set self.compressed_data with the 1447 # pre-allocated space for the compression data; this is something I 1448 # might do away with in the future 1449 self._update_compressed_data() 1450 1451 return self.compressed_data 1452 1453 @compressed_data.deleter 1454 def compressed_data(self): 1455 # Deleting the compressed_data attribute has to be handled 1456 # with a little care to prevent a reference leak 1457 # First delete the ._coldefs attributes under it to break a possible 1458 # reference cycle 1459 if 'compressed_data' in self.__dict__: 1460 del self.__dict__['compressed_data']._coldefs 1461 1462 # Now go ahead and delete from self.__dict__; normally 1463 # lazyproperty.__delete__ does this for us, but we can prempt it to 1464 # do some additional cleanup 1465 del self.__dict__['compressed_data'] 1466 1467 # If this file was mmap'd, numpy.memmap will hold open a file 1468 # handle until the underlying mmap object is garbage-collected; 1469 # since this reference leak can sometimes hang around longer than 1470 # welcome go ahead and force a garbage collection 1471 gc.collect() 1472 1473 @property 1474 def shape(self): 1475 """ 1476 Shape of the image array--should be equivalent to ``self.data.shape``. 1477 """ 1478 1479 # Determine from the values read from the header 1480 return tuple(reversed(self._axes)) 1481 1482 @lazyproperty 1483 def header(self): 1484 # The header attribute is the header for the image data. It 1485 # is not actually stored in the object dictionary. Instead, 1486 # the _image_header is stored. If the _image_header attribute 1487 # has already been defined we just return it. If not, we must 1488 # create it from the table header (the _header attribute). 1489 if hasattr(self, '_image_header'): 1490 return self._image_header 1491 1492 # Clean up any possible doubled EXTNAME keywords that use 1493 # the default. Do this on the original header to ensure 1494 # duplicates are removed cleanly. 1495 self._remove_unnecessary_default_extnames(self._header) 1496 1497 # Start with a copy of the table header. 1498 image_header = self._header.copy() 1499 1500 # Delete cards that are related to the table. And move 1501 # the values of those cards that relate to the image from 1502 # their corresponding table cards. These include 1503 # ZBITPIX -> BITPIX, ZNAXIS -> NAXIS, and ZNAXISn -> NAXISn. 1504 # (Note: Used set here instead of list in case there are any duplicate 1505 # keywords, which there may be in some pathological cases: 1506 # https://github.com/astropy/astropy/issues/2750 1507 for keyword in set(image_header): 1508 if CompImageHeader._is_reserved_keyword(keyword, warn=False): 1509 del image_header[keyword] 1510 1511 if 'ZSIMPLE' in self._header: 1512 image_header.set('SIMPLE', self._header['ZSIMPLE'], 1513 self._header.comments['ZSIMPLE'], before=0) 1514 elif 'ZTENSION' in self._header: 1515 if self._header['ZTENSION'] != 'IMAGE': 1516 warnings.warn("ZTENSION keyword in compressed " 1517 "extension != 'IMAGE'", AstropyUserWarning) 1518 image_header.set('XTENSION', 'IMAGE', 1519 self._header.comments['ZTENSION'], before=0) 1520 else: 1521 image_header.set('XTENSION', 'IMAGE', before=0) 1522 1523 image_header.set('BITPIX', self._header['ZBITPIX'], 1524 self._header.comments['ZBITPIX'], before=1) 1525 1526 image_header.set('NAXIS', self._header['ZNAXIS'], 1527 self._header.comments['ZNAXIS'], before=2) 1528 1529 last_naxis = 'NAXIS' 1530 for idx in range(image_header['NAXIS']): 1531 znaxis = 'ZNAXIS' + str(idx + 1) 1532 naxis = znaxis[1:] 1533 image_header.set(naxis, self._header[znaxis], 1534 self._header.comments[znaxis], 1535 after=last_naxis) 1536 last_naxis = naxis 1537 1538 # Delete any other spurious NAXISn keywords: 1539 naxis = image_header['NAXIS'] 1540 for keyword in list(image_header['NAXIS?*']): 1541 try: 1542 n = int(keyword[5:]) 1543 except Exception: 1544 continue 1545 1546 if n > naxis: 1547 del image_header[keyword] 1548 1549 # Although PCOUNT and GCOUNT are considered mandatory for IMAGE HDUs, 1550 # ZPCOUNT and ZGCOUNT are optional, probably because for IMAGE HDUs 1551 # their values are always 0 and 1 respectively 1552 if 'ZPCOUNT' in self._header: 1553 image_header.set('PCOUNT', self._header['ZPCOUNT'], 1554 self._header.comments['ZPCOUNT'], 1555 after=last_naxis) 1556 else: 1557 image_header.set('PCOUNT', 0, after=last_naxis) 1558 1559 if 'ZGCOUNT' in self._header: 1560 image_header.set('GCOUNT', self._header['ZGCOUNT'], 1561 self._header.comments['ZGCOUNT'], 1562 after='PCOUNT') 1563 else: 1564 image_header.set('GCOUNT', 1, after='PCOUNT') 1565 1566 if 'ZEXTEND' in self._header: 1567 image_header.set('EXTEND', self._header['ZEXTEND'], 1568 self._header.comments['ZEXTEND']) 1569 1570 if 'ZBLOCKED' in self._header: 1571 image_header.set('BLOCKED', self._header['ZBLOCKED'], 1572 self._header.comments['ZBLOCKED']) 1573 1574 # Move the ZHECKSUM and ZDATASUM cards to the image header 1575 # as CHECKSUM and DATASUM 1576 if 'ZHECKSUM' in self._header: 1577 image_header.set('CHECKSUM', self._header['ZHECKSUM'], 1578 self._header.comments['ZHECKSUM']) 1579 1580 if 'ZDATASUM' in self._header: 1581 image_header.set('DATASUM', self._header['ZDATASUM'], 1582 self._header.comments['ZDATASUM']) 1583 1584 # Remove the EXTNAME card if the value in the table header 1585 # is the default value of COMPRESSED_IMAGE. 1586 if ('EXTNAME' in image_header and 1587 image_header['EXTNAME'] == self._default_name): 1588 del image_header['EXTNAME'] 1589 1590 # Look to see if there are any blank cards in the table 1591 # header. If there are, there should be the same number 1592 # of blank cards in the image header. Add blank cards to 1593 # the image header to make it so. 1594 table_blanks = self._header._countblanks() 1595 image_blanks = image_header._countblanks() 1596 1597 for _ in range(table_blanks - image_blanks): 1598 image_header.append() 1599 1600 # Create the CompImageHeader that syncs with the table header, and save 1601 # it off to self._image_header so it can be referenced later 1602 # unambiguously 1603 self._image_header = CompImageHeader(self._header, image_header) 1604 1605 return self._image_header 1606 1607 def _summary(self): 1608 """ 1609 Summarize the HDU: name, dimensions, and formats. 1610 """ 1611 class_name = self.__class__.__name__ 1612 1613 # if data is touched, use data info. 1614 if self._data_loaded: 1615 if self.data is None: 1616 _shape, _format = (), '' 1617 else: 1618 1619 # the shape will be in the order of NAXIS's which is the 1620 # reverse of the numarray shape 1621 _shape = list(self.data.shape) 1622 _format = self.data.dtype.name 1623 _shape.reverse() 1624 _shape = tuple(_shape) 1625 _format = _format[_format.rfind('.') + 1:] 1626 1627 # if data is not touched yet, use header info. 1628 else: 1629 _shape = () 1630 1631 for idx in range(self.header['NAXIS']): 1632 _shape += (self.header['NAXIS' + str(idx + 1)],) 1633 1634 _format = BITPIX2DTYPE[self.header['BITPIX']] 1635 1636 return (self.name, self.ver, class_name, len(self.header), _shape, 1637 _format) 1638 1639 def _update_compressed_data(self): 1640 """ 1641 Compress the image data so that it may be written to a file. 1642 """ 1643 1644 # Check to see that the image_header matches the image data 1645 image_bitpix = DTYPE2BITPIX[self.data.dtype.name] 1646 1647 if image_bitpix != self._orig_bitpix or self.data.shape != self.shape: 1648 self._update_header_data(self.header) 1649 1650 # TODO: This is copied right out of _ImageBaseHDU._writedata_internal; 1651 # it would be cool if we could use an internal ImageHDU and use that to 1652 # write to a buffer for compression or something. See ticket #88 1653 # deal with unsigned integer 16, 32 and 64 data 1654 old_data = self.data 1655 if _is_pseudo_integer(self.data.dtype): 1656 # Convert the unsigned array to signed 1657 self.data = np.array( 1658 self.data - _pseudo_zero(self.data.dtype), 1659 dtype=f'=i{self.data.dtype.itemsize}') 1660 should_swap = False 1661 else: 1662 should_swap = not self.data.dtype.isnative 1663 1664 if should_swap: 1665 1666 if self.data.flags.writeable: 1667 self.data.byteswap(True) 1668 else: 1669 # For read-only arrays, there is no way around making 1670 # a byteswapped copy of the data. 1671 self.data = self.data.byteswap(False) 1672 1673 try: 1674 nrows = self._header['NAXIS2'] 1675 tbsize = self._header['NAXIS1'] * nrows 1676 1677 self._header['PCOUNT'] = 0 1678 if 'THEAP' in self._header: 1679 del self._header['THEAP'] 1680 self._theap = tbsize 1681 1682 # First delete the original compressed data, if it exists 1683 del self.compressed_data 1684 1685 # Make sure that the data is contiguous otherwise CFITSIO 1686 # will not write the expected data 1687 self.data = np.ascontiguousarray(self.data) 1688 1689 # Compress the data. 1690 # The current implementation of compress_hdu assumes the empty 1691 # compressed data table has already been initialized in 1692 # self.compressed_data, and writes directly to it 1693 # compress_hdu returns the size of the heap for the written 1694 # compressed image table 1695 heapsize, self.compressed_data = compression.compress_hdu(self) 1696 finally: 1697 # if data was byteswapped return it to its original order 1698 if should_swap: 1699 self.data.byteswap(True) 1700 self.data = old_data 1701 1702 # CFITSIO will write the compressed data in big-endian order 1703 dtype = self.columns.dtype.newbyteorder('>') 1704 buf = self.compressed_data 1705 compressed_data = buf[:self._theap].view(dtype=dtype, 1706 type=np.rec.recarray) 1707 self.compressed_data = compressed_data.view(FITS_rec) 1708 self.compressed_data._coldefs = self.columns 1709 self.compressed_data._heapoffset = self._theap 1710 self.compressed_data._heapsize = heapsize 1711 1712 def scale(self, type=None, option='old', bscale=1, bzero=0): 1713 """ 1714 Scale image data by using ``BSCALE`` and ``BZERO``. 1715 1716 Calling this method will scale ``self.data`` and update the keywords of 1717 ``BSCALE`` and ``BZERO`` in ``self._header`` and ``self._image_header``. 1718 This method should only be used right before writing to the output 1719 file, as the data will be scaled and is therefore not very usable after 1720 the call. 1721 1722 Parameters 1723 ---------- 1724 1725 type : str, optional 1726 destination data type, use a string representing a numpy dtype 1727 name, (e.g. ``'uint8'``, ``'int16'``, ``'float32'`` etc.). If is 1728 `None`, use the current data type. 1729 1730 option : str, optional 1731 how to scale the data: if ``"old"``, use the original ``BSCALE`` 1732 and ``BZERO`` values when the data was read/created. If 1733 ``"minmax"``, use the minimum and maximum of the data to scale. 1734 The option will be overwritten by any user-specified bscale/bzero 1735 values. 1736 1737 bscale, bzero : int, optional 1738 user specified ``BSCALE`` and ``BZERO`` values. 1739 """ 1740 1741 if self.data is None: 1742 return 1743 1744 # Determine the destination (numpy) data type 1745 if type is None: 1746 type = BITPIX2DTYPE[self._bitpix] 1747 _type = getattr(np, type) 1748 1749 # Determine how to scale the data 1750 # bscale and bzero takes priority 1751 if (bscale != 1 or bzero != 0): 1752 _scale = bscale 1753 _zero = bzero 1754 else: 1755 if option == 'old': 1756 _scale = self._orig_bscale 1757 _zero = self._orig_bzero 1758 elif option == 'minmax': 1759 if isinstance(_type, np.floating): 1760 _scale = 1 1761 _zero = 0 1762 else: 1763 _min = np.minimum.reduce(self.data.flat) 1764 _max = np.maximum.reduce(self.data.flat) 1765 1766 if _type == np.uint8: # uint8 case 1767 _zero = _min 1768 _scale = (_max - _min) / (2. ** 8 - 1) 1769 else: 1770 _zero = (_max + _min) / 2. 1771 1772 # throw away -2^N 1773 _scale = (_max - _min) / (2. ** (8 * _type.bytes) - 2) 1774 1775 # Do the scaling 1776 if _zero != 0: 1777 # We have to explicitly cast self._bzero to prevent numpy from 1778 # raising an error when doing self.data -= _zero, and we 1779 # do this instead of self.data = self.data - _zero to 1780 # avoid doubling memory usage. 1781 np.subtract(self.data, _zero, out=self.data, casting='unsafe') 1782 self.header['BZERO'] = _zero 1783 else: 1784 # Delete from both headers 1785 for header in (self.header, self._header): 1786 with suppress(KeyError): 1787 del header['BZERO'] 1788 1789 if _scale != 1: 1790 self.data /= _scale 1791 self.header['BSCALE'] = _scale 1792 else: 1793 for header in (self.header, self._header): 1794 with suppress(KeyError): 1795 del header['BSCALE'] 1796 1797 if self.data.dtype.type != _type: 1798 self.data = np.array(np.around(self.data), dtype=_type) # 0.7.7.1 1799 1800 # Update the BITPIX Card to match the data 1801 self._bitpix = DTYPE2BITPIX[self.data.dtype.name] 1802 self._bzero = self.header.get('BZERO', 0) 1803 self._bscale = self.header.get('BSCALE', 1) 1804 # Update BITPIX for the image header specifically 1805 # TODO: Make this more clear by using self._image_header, but only once 1806 # this has been fixed so that the _image_header attribute is guaranteed 1807 # to be valid 1808 self.header['BITPIX'] = self._bitpix 1809 1810 # Update the table header to match the scaled data 1811 self._update_header_data(self.header) 1812 1813 # Since the image has been manually scaled, the current 1814 # bitpix/bzero/bscale now serve as the 'original' scaling of the image, 1815 # as though the original image has been completely replaced 1816 self._orig_bitpix = self._bitpix 1817 self._orig_bzero = self._bzero 1818 self._orig_bscale = self._bscale 1819 1820 def _prewriteto(self, checksum=False, inplace=False): 1821 if self._scale_back: 1822 self.scale(BITPIX2DTYPE[self._orig_bitpix]) 1823 1824 if self._has_data: 1825 self._update_compressed_data() 1826 1827 # Use methods in the superclass to update the header with 1828 # scale/checksum keywords based on the data type of the image data 1829 self._update_pseudo_int_scale_keywords() 1830 1831 # Shove the image header and data into a new ImageHDU and use that 1832 # to compute the image checksum 1833 image_hdu = ImageHDU(data=self.data, header=self.header) 1834 image_hdu._update_checksum(checksum) 1835 if 'CHECKSUM' in image_hdu.header: 1836 # This will also pass through to the ZHECKSUM keyword and 1837 # ZDATASUM keyword 1838 self._image_header.set('CHECKSUM', 1839 image_hdu.header['CHECKSUM'], 1840 image_hdu.header.comments['CHECKSUM']) 1841 if 'DATASUM' in image_hdu.header: 1842 self._image_header.set('DATASUM', image_hdu.header['DATASUM'], 1843 image_hdu.header.comments['DATASUM']) 1844 # Store a temporary backup of self.data in a different attribute; 1845 # see below 1846 self._imagedata = self.data 1847 1848 # Now we need to perform an ugly hack to set the compressed data as 1849 # the .data attribute on the HDU so that the call to _writedata 1850 # handles it properly 1851 self.__dict__['data'] = self.compressed_data 1852 1853 return super()._prewriteto(checksum=checksum, inplace=inplace) 1854 1855 def _writeheader(self, fileobj): 1856 """ 1857 Bypasses `BinTableHDU._writeheader()` which updates the header with 1858 metadata about the data that is meaningless here; another reason 1859 why this class maybe shouldn't inherit directly from BinTableHDU... 1860 """ 1861 1862 return ExtensionHDU._writeheader(self, fileobj) 1863 1864 def _writedata(self, fileobj): 1865 """ 1866 Wrap the basic ``_writedata`` method to restore the ``.data`` 1867 attribute to the uncompressed image data in the case of an exception. 1868 """ 1869 1870 try: 1871 return super()._writedata(fileobj) 1872 finally: 1873 # Restore the .data attribute to its rightful value (if any) 1874 if hasattr(self, '_imagedata'): 1875 self.__dict__['data'] = self._imagedata 1876 del self._imagedata 1877 else: 1878 del self.data 1879 1880 def _close(self, closed=True): 1881 super()._close(closed=closed) 1882 1883 # Also make sure to close access to the compressed data mmaps 1884 if (closed and self._data_loaded and 1885 _get_array_mmap(self.compressed_data) is not None): 1886 del self.compressed_data 1887 1888 # TODO: This was copied right out of _ImageBaseHDU; get rid of it once we 1889 # find a way to rewrite this class as either a subclass or wrapper for an 1890 # ImageHDU 1891 def _dtype_for_bitpix(self): 1892 """ 1893 Determine the dtype that the data should be converted to depending on 1894 the BITPIX value in the header, and possibly on the BSCALE value as 1895 well. Returns None if there should not be any change. 1896 """ 1897 1898 bitpix = self._orig_bitpix 1899 # Handle possible conversion to uints if enabled 1900 if self._uint and self._orig_bscale == 1: 1901 for bits, dtype in ((16, np.dtype('uint16')), 1902 (32, np.dtype('uint32')), 1903 (64, np.dtype('uint64'))): 1904 if bitpix == bits and self._orig_bzero == 1 << (bits - 1): 1905 return dtype 1906 1907 if bitpix > 16: # scale integers to Float64 1908 return np.dtype('float64') 1909 elif bitpix > 0: # scale integers to Float32 1910 return np.dtype('float32') 1911 1912 def _update_header_scale_info(self, dtype=None): 1913 if (not self._do_not_scale_image_data and 1914 not (self._orig_bzero == 0 and self._orig_bscale == 1)): 1915 for keyword in ['BSCALE', 'BZERO']: 1916 # Make sure to delete from both the image header and the table 1917 # header; later this will be streamlined 1918 for header in (self.header, self._header): 1919 with suppress(KeyError): 1920 del header[keyword] 1921 # Since _update_header_scale_info can, currently, be 1922 # called *after* _prewriteto(), replace these with 1923 # blank cards so the header size doesn't change 1924 header.append() 1925 1926 if dtype is None: 1927 dtype = self._dtype_for_bitpix() 1928 if dtype is not None: 1929 self.header['BITPIX'] = DTYPE2BITPIX[dtype.name] 1930 1931 self._bzero = 0 1932 self._bscale = 1 1933 self._bitpix = self.header['BITPIX'] 1934 1935 def _generate_dither_seed(self, seed): 1936 if not _is_int(seed): 1937 raise TypeError("Seed must be an integer") 1938 1939 if not -1 <= seed <= 10000: 1940 raise ValueError( 1941 "Seed for random dithering must be either between 1 and " 1942 "10000 inclusive, 0 for autogeneration from the system " 1943 "clock, or -1 for autogeneration from a checksum of the first " 1944 "image tile (got {})".format(seed)) 1945 1946 if seed == DITHER_SEED_CHECKSUM: 1947 # Determine the tile dimensions from the ZTILEn keywords 1948 naxis = self._header['ZNAXIS'] 1949 tile_dims = [self._header[f'ZTILE{idx + 1}'] 1950 for idx in range(naxis)] 1951 tile_dims.reverse() 1952 1953 # Get the first tile by using the tile dimensions as the end 1954 # indices of slices (starting from 0) 1955 first_tile = self.data[tuple(slice(d) for d in tile_dims)] 1956 1957 # The checksum algorithm used is literally just the sum of the bytes 1958 # of the tile data (not its actual floating point values). Integer 1959 # overflow is irrelevant. 1960 csum = first_tile.view(dtype='uint8').sum() 1961 1962 # Since CFITSIO uses an unsigned long (which may be different on 1963 # different platforms) go ahead and truncate the sum to its 1964 # unsigned long value and take the result modulo 10000 1965 return (ctypes.c_ulong(csum).value % 10000) + 1 1966 elif seed == DITHER_SEED_CLOCK: 1967 # This isn't exactly the same algorithm as CFITSIO, but that's okay 1968 # since the result is meant to be arbitrary. The primary difference 1969 # is that CFITSIO incorporates the HDU number into the result in 1970 # the hopes of heading off the possibility of the same seed being 1971 # generated for two HDUs at the same time. Here instead we just 1972 # add in the HDU object's id 1973 return ((sum(int(x) for x in math.modf(time.time())) + id(self)) % 1974 10000) + 1 1975 else: 1976 return seed 1977