Merge pull request #169 from elvinpoole/hdf5

Add support for hdf5 files

Author: erykoff

.github/workflows/python-package-pr.yml

@@ -24,7 +24,7 @@ jobs:
     - name: Install dependencies
       run: |
         python -m pip install --upgrade pip
-        pip install .[parquet,test,test_with_healpy]
+        pip install .[parquet,test,test_with_healpy,hdf5]
     - name: Lint with flake8
       run: |
         flake8

docs/basic_interface.rst

@@ -208,12 +208,17 @@ Writing a :code:`HealSparseMap` is easy.  To write a map in the default FITS for
 
     map3.write('output_file.hs', clobber=False)
 
-And to write a map in the Parquet format with ``pyarrow``:
+to write a map in the Parquet format with ``pyarrow``:
 
 .. code-block :: python
 
     map3.write('output_file.hsparquet', clobber=False, format='parquet')
 
+and to write a map in the HDF5 format with ``h5py``:
+
+.. code-block :: python
+
+    map3.write('output_file.h5', clobber=False, format='hdf5')
 
 Metadata
 --------

docs/filespec.rst

@@ -276,3 +276,123 @@ The other columns in the overflow coverage pixel are filled with the default sen
 
 If the sparse map is a bit-packed mask, the schema is the same as for a regular sparse map image.
 In this case, as with the FITS serialization, the sparse map is stored as an array of unsigned 8-bit integers which is the in-memory backing of the bit-packed array.
+
+.. _hdf5_format:
+
+HealSparseMap HDF5 Serialization
+================================
+
+A :code:`HealSparseMap` may also be serialized to an HDF5 file. 
+Multiple :code:`HealSparseMap` objects can be stored in the same HDF5 file.
+Each :code:`HealSparseMap` is stored in a different HDF5 group.
+They are not required to have the same mask, :code:`nside_coverage`, or :code:`nside_sparse`.
+
+All datasets are written with gzip compression and chunked such that each chunk corresponds to a single coverage pixel.
+
+HDF5 Group Layout
+-----------------
+
+A serialized :code:`HealSparseMap` is stored within a single HDF5 group (default name :code:`"map"`), which contains:
+
+* A dataset encoding the coverage map
+* One or more datasets encoding the sparse map
+* Attributes storing map metadata
+
+Coverage Map
+------------
+
+The coverage map is stored as a one-dimensional dataset:
+
+* **Dataset name:** :code:`cov_index_map`
+* **Datatype:** :code:`numpy.int64`
+* **Shape:** :code:`(12 * nside_coverage * nside_coverage,)`
+
+This dataset is a direct serialization of the in-memory coverage index map,
+following the same conventions described in :ref:`coverage_map`.
+
+Sparse Map
+----------
+
+The sparse map is stored as a two-dimensional dataset, where each row corresponds to a single coverage pixel.
+The first row always represents the overflow (sentinel) coverage pixel.
+
+Regular Sparse Map
+^^^^^^^^^^^^^^^^^^
+
+For regular (non-record, non-wide-mask, non-bit-packed) sparse maps:
+
+* **Dataset name:** :code:`sparse_map`
+* **Shape:** :code:`(ncov_in_sparse, nfine_per_cov)`
+* **Datatype:** Sparse map datatype
+
+where:
+
+* :code:`ncov_in_sparse`: The number of coverage pixels containing non-sentinel data, plus one for the overflow pixel
+ 
+Each row contains the :code:`nfine_per_cov` sparse pixel values associated with a single coverage pixel.
+The dataset is chunked as :code:`(1, nfine_per_cov)` to allow efficient access to individual coverage pixels.
+
+Sparse Map Record Array
+^^^^^^^^^^^^^^^^^^^^^^^
+
+If the sparse map is a numpy record array, each field is stored in a separate subgroup:
+
+* **Group name:** Name of the record field
+* **Dataset name:** :code:`sparse_map`
+* **Shape:** :code:`(ncov_in_sparse, nfine_per_cov)`
+* **Datatype:** Field datatype
+
+All fields share the same coverage map indexing.
+
+Sparse Map Wide Mask
+^^^^^^^^^^^^^^^^^^^^
+
+For wide mask maps, the sparse map is stored as a three-dimensional dataset:
+
+* **Dataset name:** :code:`sparse_map`
+* **Shape:** :code:`(ncov_in_sparse, nfine_per_cov, wide_mask_width)`
+* **Datatype:** :code:`numpy.uint8`
+
+Each sparse pixel is represented by :code:`wide_mask_width` bytes.
+The dataset is chunked as :code:`(1, nfine_per_cov, wide_mask_width)`.
+
+Sparse Map Bit-Packed Mask
+^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+For bit-packed mask maps, the sparse map is stored in its packed representation:
+
+* **Dataset name:** :code:`sparse_map`
+* **Shape:** :code:`(ncov_in_sparse, nfine_per_cov // 8)`
+* **Datatype:** :code:`numpy.uint8`
+
+Each byte encodes eight sparse pixels.
+
+Metadata and Attributes
+-----------------------
+
+All remaining map metadata is stored as HDF5 group attributes.
+The following attributes are required:
+
+* :code:`nside_sparse`
+* :code:`nside_coverage`
+* :code:`sentinel`
+* :code:`primary`
+* :code:`nest` (always :code:`True`)
+
+The following flags describe the sparse map type:
+
+* :code:`is_rec_array`
+* :code:`is_bit_packed`
+* :code:`is_wide_mask`
+* :code:`wide_mask_width`
+
+Any additional metadata attached to the :code:`HealSparseMap` object is also written as group attributes.
+
+Partial Reads
+-------------
+
+The HDF5 format supports reading a subset of coverage pixels.
+When a subset is requested, only the corresponding rows of the sparse map datasets are read.
+The overflow (sentinel) coverage pixel is always included as the first row.
+
+On read, the sparse map rows are reshaped and reordered to reconstruct the in-memory sparse map layout described in :ref:`sparse_map`.

healsparse/healSparseMap.py

@@ -119,7 +119,7 @@ def __init__(self, cov_map=None, cov_index_map=None, sparse_map=None, nside_spar
     @classmethod
     def read(cls, filename, nside_coverage=None, pixels=None, header=False,
              degrade_nside=None, weightfile=None, reduction='mean',
-             use_threads=False):
+             use_threads=False, hdf5_group='map'):
         """
         Read in a HealSparseMap.
 
@@ -149,6 +149,8 @@ def read(cls, filename, nside_coverage=None, pixels=None, header=False,
            Not yet implemented for parquet files.
         use_threads : `bool`, optional
            Use multithreaded reading for parquet files.
+        hdf5_group : `str`, optional
+           If file format is ``hdf5``, read from this group, otherwise unused.
 
         Returns
         -------
@@ -159,7 +161,8 @@ def read(cls, filename, nside_coverage=None, pixels=None, header=False,
         """
         return _read_map(cls, filename, nside_coverage=nside_coverage, pixels=pixels,
                          header=header, degrade_nside=degrade_nside,
-                         weightfile=weightfile, reduction=reduction, use_threads=use_threads)
+                         weightfile=weightfile, reduction=reduction, use_threads=use_threads,
+                         hdf5_group=hdf5_group)
 
     @classmethod
     def make_empty(cls, nside_coverage, nside_sparse, dtype, primary=None, sentinel=None,
@@ -359,7 +362,7 @@ def convert_healpix_map(healpix_map, nside_coverage, nest=True, sentinel=hpg.UNS
 
         return cov_map, sparse_map
 
-    def write(self, filename, clobber=False, nocompress=False, format='fits', nside_io=4):
+    def write(self, filename, clobber=False, nocompress=False, format='fits', nside_io=4, hdf5_group='map'):
         """
         Write a HealSparseMap to a file.  Use the `metadata` property from
         the map to persist additional information in the fits header.
@@ -379,17 +382,19 @@ def write(self, filename, clobber=False, nocompress=False, format='fits', nside_
             Must be less than or equal to nside_coverage, and not greater than 16.
             This option only applies if format=``parquet``.
         format : `str`, optional
-            File format.  May be ``fits``, ``parquet``, or ``healpix``. Note that
+            File format.  May be ``fits``, ``parquet``, ``hdf5`` or ``healpix``. Note that
             the ``healpix`` EXPLICIT format does not maintain all metadata and
             coverage information.
+        hdf5_group: `str`, optional
+            If file format is ``hdf5``, save to this group, otherwise unused.
 
         Raises
         ------
         NotImplementedError if file format is not supported.
         ValueError if nside_io is out of range.
         """
         _write_map(self, filename, clobber=clobber, nocompress=nocompress, format=format,
-                   nside_io=nside_io)
+                   nside_io=nside_io, hdf5_group=hdf5_group)
 
     def write_moc(self, filename, clobber=False):
         """

healsparse/io_map.py

@@ -5,10 +5,11 @@
 from .parquet_shim import check_parquet_dataset
 from .io_map_parquet import _read_map_parquet, _write_map_parquet
 from .io_map_healpix import _write_map_healpix
+from .io_map_hdf5 import _read_map_hdf5, _write_map_hdf5, check_hdf5_file
 
 
 def _read_map(healsparse_class, filename, nside_coverage=None, pixels=None, header=False,
-              degrade_nside=None, weightfile=None, reduction='mean', use_threads=False):
+              degrade_nside=None, weightfile=None, reduction='mean', use_threads=False, hdf5_group='map'):
     """
     Internal function to check the map filetype and read in a HealSparseMap.
 
@@ -37,6 +38,8 @@ def _read_map(healsparse_class, filename, nside_coverage=None, pixels=None, head
         (mean, median, std, max, min, and, or, sum, prod, wmean).
     use_threads : `bool`, optional
         Use multithreaded reading for parquet files.
+    hdf5_group: `str`, optional
+        Read from this hdf5 group (only used reading from an hdf5 file).
 
     Returns
     -------
@@ -47,17 +50,22 @@ def _read_map(healsparse_class, filename, nside_coverage=None, pixels=None, head
     """
     is_fits_file = False
     is_parquet_file = False
+    is_hdf5_file = False
 
     try:
         fits = HealSparseFits(filename)
         is_fits_file = True
         fits.close()
-    except OSError:
+    except (OSError, UnicodeDecodeError):
         pass
+    # UnicodeDecodeError occurs when trying to read an hdf5 file as if it's FITS
 
     if not is_fits_file:
         is_parquet_file = check_parquet_dataset(filename)
 
+    if not is_fits_file and not is_parquet_file:
+        is_hdf5_file = check_hdf5_file(filename)
+
     if is_fits_file:
         return _read_map_fits(healsparse_class, filename, nside_coverage=nside_coverage,
                               pixels=pixels, header=header, degrade_nside=degrade_nside,
@@ -67,13 +75,19 @@ def _read_map(healsparse_class, filename, nside_coverage=None, pixels=None, head
                                  pixels=pixels, header=header, degrade_nside=degrade_nside,
                                  weightfile=weightfile, reduction=reduction,
                                  use_threads=use_threads)
+    elif is_hdf5_file:
+        return _read_map_hdf5(healsparse_class, filename, hdf5_group=hdf5_group,
+                              pixels=pixels, header=header, degrade_nside=degrade_nside,
+                              weightfile=weightfile, reduction=reduction)
     elif not os.path.isfile(filename):
         raise IOError("Filename %s could not be found." % (filename))
     else:
-        raise NotImplementedError("HealSparse only supports fits and parquet files (with pyarrow).")
+        raise NotImplementedError(
+            "HealSparse only supports fits, hdf5 (with h5py), and parquet files (with pyarrow).")
 
 
-def _write_map(hsp_map, filename, clobber=False, nocompress=False, format='fits', nside_io=4):
+def _write_map(hsp_map, filename, clobber=False, nocompress=False, format='fits', nside_io=4,
+               hdf5_group='map'):
     """
     Internal method to write a HealSparseMap to a file, and check formats.
     Use the `metadata` property from
@@ -96,7 +110,7 @@ def _write_map(hsp_map, filename, clobber=False, nocompress=False, format='fits'
         This option only applies if format=``parquet``.
         Must be less than or equal to nside_coverage, and not greater than 16.
     format : `str`, optional
-        File format.  May be ``fits``, ``parquet``, or ``healpix``.
+        File format.  May be ``fits``, ``parquet``, ``hdf5`` or ``healpix``.
 
     Raises
     ------
@@ -109,8 +123,10 @@ def _write_map(hsp_map, filename, clobber=False, nocompress=False, format='fits'
         _write_map_parquet(hsp_map, filename, clobber=clobber, nside_io=nside_io)
     elif format == 'healpix':
         _write_map_healpix(hsp_map, filename, clobber=clobber)
+    elif format == 'hdf5':
+        _write_map_hdf5(hsp_map, filename, clobber=clobber, hdf5_group=hdf5_group)
     else:
-        raise NotImplementedError("Only 'fits', 'parquet' and 'healpix' file formats are supported.")
+        raise NotImplementedError("Only 'fits', 'parquet', 'hdf5' and 'healpix' file formats are supported.")
 
 
 def _write_moc(hsp_map, filename, clobber=False):