[Follow up PR] Make HCOMPRESS Support tiles not along the first two dimensions

.gitattributes

@@ -1 +1,2 @@
 *.fits    -text
+astropy/io/fits/tiled_compression/tests/data/*fits binary

.github/workflows/ci_workflows.yml

@@ -45,7 +45,7 @@ jobs:
           }
 
   tests:
-    needs: [initial_checks]
+    # needs: [initial_checks]
     uses: OpenAstronomy/github-actions-workflows/.github/workflows/tox.yml@v1
     with:
       submodules: false
@@ -65,7 +65,12 @@ jobs:
         # leave open any files. This has a performance impact on running the
         # tests, hence why it is not enabled by default.
         - name: Python 3.9 with all optional dependencies
-          linux: py39-test-alldeps
+          linux: py39-test-alldeps-fitsio
+          libraries:
+            apt:
+              - language-pack-fr
+              - tzdata
+              - libcfitsio-dev
           toxargs: -v --develop
           posargs: --open-files --run-slow
 
@@ -82,7 +87,7 @@ jobs:
           posargs: --durations=50 --run-slow
 
   allowed_failures:
-    needs: [initial_checks]
+    # needs: [initial_checks]
     uses: OpenAstronomy/github-actions-workflows/.github/workflows/tox.yml@v1
     with:
       submodules: false
@@ -101,7 +106,7 @@ jobs:
   parallel_and_32bit:
     name: 32-bit and parallel
     runs-on: ubuntu-latest
-    needs: [initial_checks]
+    # needs: [initial_checks]
     container:
       image: quay.io/pypa/manylinux2014_i686:2022-04-03-da6ecb3
     steps:

astropy/io/fits/hdu/compressed.py

@@ -17,6 +17,7 @@
 from astropy.io.fits.column import TDEF_RE, ColDefs, Column
 from astropy.io.fits.fitsrec import FITS_rec
 from astropy.io.fits.header import Header
+from astropy.io.fits.tiled_compression import compress_hdu, decompress_hdu
 from astropy.io.fits.util import (
     _get_array_mmap,
     _is_int,
@@ -30,13 +31,7 @@
 from .image import ImageHDU
 from .table import BinTableHDU
 
-try:
-    from astropy.io.fits import compression
-
-    COMPRESSION_SUPPORTED = COMPRESSION_ENABLED = True
-except ImportError:
-    COMPRESSION_SUPPORTED = COMPRESSION_ENABLED = False
-
+COMPRESSION_ENABLED = True
 
 # Quantization dithering method constants; these are right out of fitsio.h
 NO_DITHER = -1
@@ -642,14 +637,6 @@ def __init__(
         same image will always use the same seed.
         """
 
-        if not COMPRESSION_SUPPORTED:
-            # TODO: Raise a more specific Exception type
-            raise Exception(
-                "The astropy.io.fits.compression module is not "
-                "available.  Creation of compressed image HDUs is "
-                "disabled."
-            )
-
         compression_type = CMTYPE_ALIASES.get(compression_type, compression_type)
 
         if data is DELAYED:
@@ -764,16 +751,8 @@ def match_header(cls, header):
         if "ZIMAGE" not in header or not header["ZIMAGE"]:
             return False
 
-        if COMPRESSION_SUPPORTED and COMPRESSION_ENABLED:
+        if COMPRESSION_ENABLED:
             return True
-        elif not COMPRESSION_SUPPORTED:
-            warnings.warn(
-                "Failure matching header to a compressed image "
-                "HDU: The compression module is not available.\n"
-                "The HDU will be treated as a Binary Table HDU.",
-                AstropyUserWarning,
-            )
-            return False
         else:
             # Compression is supported but disabled; just pass silently (#92)
             return False
@@ -1068,10 +1047,26 @@ def _update_header_data(
         # Set default tile dimensions for HCOMPRESS_1
 
         if compression_type == "HCOMPRESS_1":
-            if self._image_header["NAXIS1"] < 4 or self._image_header["NAXIS2"] < 4:
+
+            # HCOMPRESS_1 allows >2D tiles provided only two dimensions have
+            # non-1 values. We store the indices of those dimensions to use
+            # below.
+            keep_dims = []
+            for i in range(len(tile_size)):
+                if tile_size[i] > 1:
+                    keep_dims.append(i)
+
+            if (
+                self._image_header[f"NAXIS{keep_dims[0] + 1}"] < 4
+                or self._image_header[f"NAXIS{keep_dims[1] + 1}"] < 4
+            ):
                 raise ValueError("Hcompress minimum image dimension is 4 pixels")
             elif tile_size:
-                if tile_size[0] < 4 or tile_size[1] < 4:
+                if len(keep_dims) != 2:
+                    raise ValueError(
+                        "Hcompress tile size should contain exactly two non-1 values"
+                    )
+                if tile_size[keep_dims[0]] < 4 or tile_size[keep_dims[1]] < 4:
                     # user specified tile size is too small
                     raise ValueError("Hcompress minimum tile dimension is 4 pixels")
                 major_dims = len([ts for ts in tile_size if ts > 1])
@@ -1082,6 +1077,10 @@ def _update_header_data(
                         "to 1."
                     )
 
+            # TODO: is '0' a magic value that means all pixels along a certain
+            # dimension? If so we should perhaps check not the first two but
+            # the two dimensions that are not 1.
+
             if tile_size and (tile_size[0] == 0 and tile_size[1] == 0):
                 # compress the whole image as a single tile
                 tile_size[0] = self._image_header["NAXIS1"]
@@ -1125,24 +1124,33 @@ def _update_header_data(
             # check if requested tile size causes the last tile to have
             # less than 4 pixels
 
-            remain = self._image_header["NAXIS1"] % tile_size[0]  # 1st dimen
+            remain = (
+                self._image_header[f"NAXIS{keep_dims[0] + 1}"] % tile_size[keep_dims[0]]
+            )  # 1st dimen
 
             if remain > 0 and remain < 4:
-                tile_size[0] += 1  # try increasing tile size by 1
-
-                remain = self._image_header["NAXIS1"] % tile_size[0]
+                tile_size[keep_dims[0]] += 1  # try increasing tile size by 1
+                remain = (
+                    self._image_header[f"NAXIS{keep_dims[0] + 1}"]
+                    % tile_size[keep_dims[0]]
+                )
 
                 if remain > 0 and remain < 4:
                     raise ValueError(
                         "Last tile along 1st dimension has less than 4 pixels"
                     )
 
-            remain = self._image_header["NAXIS2"] % tile_size[1]  # 2nd dimen
+            remain = (
+                self._image_header[f"NAXIS{keep_dims[1] + 1}"] % tile_size[keep_dims[1]]
+            )  # 2nd dimen
 
             if remain > 0 and remain < 4:
-                tile_size[1] += 1  # try increasing tile size by 1
+                tile_size[keep_dims[1]] += 1  # try increasing tile size by 1
 
-                remain = self._image_header["NAXIS2"] % tile_size[1]
+                remain = (
+                    self._image_header[f"NAXIS{keep_dims[1] + 1}"]
+                    % tile_size[keep_dims[1]]
+                )
 
                 if remain > 0 and remain < 4:
                     raise ValueError(
@@ -1495,7 +1503,7 @@ def _update_header_data(
     @lazyproperty
     def data(self):
         # The data attribute is the image data (not the table data).
-        data = compression.decompress_hdu(self)
+        data = decompress_hdu(self)
 
         if data is None:
             return data
@@ -1505,6 +1513,11 @@ def data(self):
             new_dtype = self._dtype_for_bitpix()
             data = np.array(data, dtype=new_dtype)
 
+            # TODO: ZBLANK should not be dealt with here as it should apply
+            # to data before being unquantized, which has already happened in
+            # decompress_hdu, so we should remove references to ZBLANK here.
+            # However, it might be that we still need to handle BLANK.
+
             zblank = None
 
             if "ZBLANK" in self.compressed_data.columns.names:
@@ -1802,17 +1815,6 @@ def _update_compressed_data(self):
                 self.data - _pseudo_zero(self.data.dtype),
                 dtype=f"=i{self.data.dtype.itemsize}",
             )
-            should_swap = False
-        else:
-            should_swap = not self.data.dtype.isnative
-
-        if should_swap:
-            if self.data.flags.writeable:
-                self.data.byteswap(True)
-            else:
-                # For read-only arrays, there is no way around making
-                # a byteswapped copy of the data.
-                self.data = self.data.byteswap(False)
 
         try:
             nrows = self._header["NAXIS2"]
@@ -1836,11 +1838,8 @@ def _update_compressed_data(self):
             # self.compressed_data, and writes directly to it
             # compress_hdu returns the size of the heap for the written
             # compressed image table
-            heapsize, self.compressed_data = compression.compress_hdu(self)
+            heapsize, self.compressed_data = compress_hdu(self)
         finally:
-            # if data was byteswapped return it to its original order
-            if should_swap:
-                self.data.byteswap(True)
             self.data = old_data
 
         # CFITSIO will write the compressed data in big-endian order