Allow alternate skymaps (and unknown properties) (#511)

* Allow alternate skymaps (and unknown properties)

* Improve parsing of list of ints.

Author: delucchi-cmu

src/hats/catalog/dataset/table_properties.py

@@ -10,24 +10,6 @@
 from hats.catalog.catalog_type import CatalogType
 from hats.io import file_io
 
-## catalog_name, catalog_type, and total_rows are allowed for ALL types
-CATALOG_TYPE_ALLOWED_FIELDS = {
-    CatalogType.OBJECT: ["ra_column", "dec_column", "default_columns"],
-    CatalogType.SOURCE: ["primary_catalog", "ra_column", "dec_column", "default_columns"],
-    CatalogType.ASSOCIATION: [
-        "primary_catalog",
-        "primary_column",
-        "primary_column_association",
-        "join_catalog",
-        "join_column",
-        "join_column_association",
-        "contains_leaf_files",
-    ],
-    CatalogType.INDEX: ["primary_catalog", "indexing_column", "extra_columns"],
-    CatalogType.MARGIN: ["primary_catalog", "margin_threshold", "ra_column", "dec_column", "default_columns"],
-    CatalogType.MAP: ["default_columns"],
-}
-
 ## catalog_name, catalog_type, and total_rows are required for ALL types
 CATALOG_TYPE_REQUIRED_FIELDS = {
     CatalogType.OBJECT: ["ra_column", "dec_column"],
@@ -44,42 +26,6 @@
     CatalogType.MAP: [],
 }
 
-# All additional properties in the HATS recommendation.
-EXTRA_ALLOWED_FIELDS = [
-    "addendum_did",
-    "bib_reference",
-    "bib_reference_url",
-    "creator_did",
-    "data_ucd",
-    "hats_builder",
-    "hats_cols_sort",
-    "hats_cols_survey_id",
-    "hats_coordinate_epoch",
-    "hats_copyright",
-    "hats_creation_date",
-    "hats_creator",
-    "hats_estsize",
-    "hats_frame",
-    "hats_max_rows",
-    "hats_order",
-    "hats_progenitor_url",
-    "hats_release_date",
-    "hats_service_url",
-    "hats_status",
-    "hats_version",
-    "moc_sky_fraction",
-    "obs_ack",
-    "obs_copyright",
-    "obs_copyright_url",
-    "obs_description",
-    "obs_regime",
-    "obs_title",
-    "prov_progenitor",
-    "publisher_id",
-    "t_max",
-    "t_min",
-]
-
 
 class TableProperties(BaseModel):
     """Container class for catalog metadata"""
@@ -123,9 +69,6 @@ class TableProperties(BaseModel):
     extra_columns: Optional[list[str]] = Field(default=None, alias="hats_index_extra_column")
     """Any additional payload columns included in index."""
 
-    ## Allow any extra keyword args to be stored on the properties object.
-    model_config = ConfigDict(extra="allow", populate_by_name=True, use_enum_values=True)
-
     npix_suffix: str = Field(default=".parquet", alias="hats_npix_suffix")
     """Suffix of the Npix partitions.
     In the standard HATS directory structure, this is '.parquet' because there is a single file
@@ -136,6 +79,15 @@ class TableProperties(BaseModel):
     underneath (and then `npix_suffix` = '/').
     """
 
+    skymap_order: Optional[int] = Field(default=None, alias="hats_skymap_order")
+    """Nested Order of the healpix skymap stored in the default skymap.fits."""
+
+    skymap_alt_orders: Optional[list[int]] = Field(default=None, alias="hats_skymap_alt_orders")
+    """Nested Order (K) of the healpix skymaps stored in altnernative skymap.K.fits."""
+
+    ## Allow any extra keyword args to be stored on the properties object.
+    model_config = ConfigDict(extra="allow", populate_by_name=True, use_enum_values=True)
+
     @field_validator("default_columns", "extra_columns", mode="before")
     @classmethod
     def space_delimited_list(cls, str_value: str) -> list[str]:
@@ -146,28 +98,51 @@ def space_delimited_list(cls, str_value: str) -> list[str]:
         ## Convert empty strings and empty lists to None
         return str_value if str_value else None
 
-    @field_serializer("default_columns", "extra_columns")
-    def serialize_as_space_delimited_list(self, str_list: Iterable[str]) -> str:
+    @field_validator("skymap_alt_orders", mode="before")
+    @classmethod
+    def space_delimited_int_list(cls, str_value: str | list[int]) -> list[int]:
+        """Convert a space-delimited list string into a python list of integers.
+
+        Args:
+            str_value(str | list[int]): string representation of a list of integers, delimited by
+                space, comma, or semi-colon, or a list of integers.
+
+        Returns:
+            sorted list of unique integers, if all inputs are integers, None if the input is empty
+        Raises:
+            ValueError: if any non-digit characters are encountered, or
+        """
+        if not str_value:
+            return None
+        if isinstance(str_value, int):
+            return [str_value]
+        if isinstance(str_value, str):
+            # Split on a few kinds of delimiters (just to be safe)
+            int_list = [int(token) for token in list(filter(None, re.split(";| |,|\n", str_value)))]
+        elif isinstance(str_value, list) and all(isinstance(elem, int) for elem in str_value):
+            int_list = str_value
+        else:
+            raise ValueError(f"Unsupported type of skymap_alt_orders {type(str_value)}")
+        if len(int_list) == 0:
+            return None
+        int_list = list(set(int_list))
+        int_list.sort()
+        return int_list
+
+    @field_serializer("default_columns", "extra_columns", "skymap_alt_orders")
+    def serialize_as_space_delimited_list(self, str_list: Iterable) -> str:
         """Convert a python list of strings into a space-delimited string."""
         if str_list is None or len(str_list) == 0:
             return None
-        return " ".join(str_list)
+        return " ".join([str(element) for element in str_list])
 
     @model_validator(mode="after")
-    def check_allowed_and_required(self) -> Self:
+    def check_required(self) -> Self:
         """Check that type-specific fields are appropriate, and required fields are set."""
         explicit_keys = set(
             self.model_dump(by_alias=False, exclude_none=True).keys() - self.__pydantic_extra__.keys()
         )
 
-        allowed_keys = set(
-            CATALOG_TYPE_ALLOWED_FIELDS[self.catalog_type]
-            + ["catalog_name", "catalog_type", "total_rows", "npix_suffix"]
-        )
-        non_allowed = explicit_keys - allowed_keys
-        if len(non_allowed) > 0:
-            raise ValueError(f"Unexpected property for table type {self.catalog_type}: {non_allowed}")
-
         required_keys = set(
             CATALOG_TYPE_REQUIRED_FIELDS[self.catalog_type] + ["catalog_name", "catalog_type", "total_rows"]
         )
@@ -177,10 +152,6 @@ def check_allowed_and_required(self) -> Self:
                 f"Missing required property for table type {self.catalog_type}: {missing_required}"
             )
 
-        # Check against all known properties - catches typos.
-        non_allowed = set(self.__pydantic_extra__.keys()) - set(EXTRA_ALLOWED_FIELDS)
-        if len(non_allowed) > 0:
-            raise ValueError(f"Unexpected extra property: {non_allowed}")
         return self
 
     def copy_and_update(self, **kwargs):

src/hats/inspection/visualize_catalog.py

@@ -27,7 +27,7 @@
 from mocpy.moc.plot.utils import _set_wcs
 
 import hats.pixel_math.healpix_shim as hp
-from hats.io import file_io, paths
+from hats.io import skymap
 from hats.pixel_math import HealpixPixel
 from hats.pixel_tree.moc_filter import perform_filter_by_moc
 from hats.pixel_tree.pixel_tree import PixelTree
@@ -37,21 +37,9 @@
     from hats.catalog.healpix_dataset.healpix_dataset import HealpixDataset
 
 
-def _read_point_map(catalog_base_dir):
-    """Read the object spatial distribution information from a healpix FITS file.
-
-    Args:
-        catalog_base_dir: path to a catalog
-    Returns:
-        one-dimensional numpy array of long integers where the value at each index
-        corresponds to the number of objects found at the healpix pixel.
-    """
-    map_file_pointer = paths.get_point_map_file_pointer(catalog_base_dir)
-    return file_io.read_fits_image(map_file_pointer)
-
-
 def plot_density(catalog: Catalog, *, plot_title: str | None = None, order=None, unit=None, **kwargs):
     """Create a visual map of the density of input points of a catalog on-disk.
+
     Args:
         catalog (`hats.catalog.Catalog`) Catalog to display
         plot_title (str): Optional title for the plot
@@ -60,16 +48,9 @@ def plot_density(catalog: Catalog, *, plot_title: str | None = None, order=None,
     """
     if catalog is None or not catalog.on_disk:
         raise ValueError("on disk catalog required for point-wise visualization")
-    point_map = _read_point_map(catalog.catalog_base_dir)
-    map_order = hp.npix2order(len(point_map))
-
-    if order is not None:
-        if order > map_order:
-            raise ValueError(f"plotting order should be less than stored density map order ({map_order})")
-        ## Create larger pixel sums from the constituent pixels.
-        point_map = point_map.reshape(hp.order2npix(order), -1).sum(axis=1)
-    else:
-        order = map_order
+    point_map = skymap.read_skymap(catalog, order)
+    order = hp.npix2order(len(point_map))
+
     if unit is None:
         unit = u.deg * u.deg
 

src/hats/io/__init__.py

@@ -6,6 +6,7 @@
     get_parquet_metadata_pointer,
     get_partition_info_pointer,
     get_point_map_file_pointer,
+    get_skymap_file_pointer,
     pixel_catalog_file,
     pixel_directory,
 )

src/hats/io/paths.py

@@ -36,6 +36,7 @@
 PARQUET_COMMON_METADATA_FILENAME = "_common_metadata"
 DATA_THUMBNAIL_FILENAME = "data_thumbnail.parquet"
 POINT_MAP_FILENAME = "point_map.fits"
+SKYMAP_FILENAME = "skymap.fits"
 
 
 def pixel_directory(
@@ -223,6 +224,20 @@ def get_point_map_file_pointer(catalog_base_dir: str | Path | UPath) -> UPath:
     return get_upath(catalog_base_dir) / POINT_MAP_FILENAME
 
 
+def get_skymap_file_pointer(catalog_base_dir: str | Path | UPath, order: int | None = None) -> UPath:
+    """Get file pointer to `skymap.fits` or `skymap.K.fits` FITS image file.
+
+    Args:
+        catalog_base_dir: pointer to base catalog directory
+        order: healpix order of the desired down-sampled skymap
+    Returns:
+        File Pointer to the FITS image file.
+    """
+    if order is not None and order >= 0:
+        return get_upath(catalog_base_dir) / f"skymap.{order}.fits"
+    return get_upath(catalog_base_dir) / SKYMAP_FILENAME
+
+
 def get_partition_join_info_pointer(catalog_base_dir: str | Path | UPath) -> UPath:
     """Get file pointer to `partition_join_info.csv` association metadata file
 

src/hats/io/skymap.py

@@ -0,0 +1,92 @@
+from pathlib import Path
+from typing import Iterable
+
+import numpy as np
+from upath import UPath
+
+import hats.pixel_math.healpix_shim as hp
+from hats.io import file_io, paths
+
+
+def read_skymap(catalog, order):
+    """Read the object spatial distribution information from a healpix skymap FITS file.
+
+    Args:
+        catalog: Catalog object corresponding to an on-disk catalog.
+        order: healpix order to read the skymap at. If None, the order of the default
+            skymap will be used. We will try to load from alternative skymap orders,
+            where appropriate.
+    Returns:
+        one-dimensional numpy array of long integers where the value at each index
+        corresponds to the number of objects found at the healpix pixel.
+    """
+    if order is not None and catalog.catalog_info.skymap_alt_orders:
+        available_orders = catalog.catalog_info.skymap_alt_orders
+        best_order_idx = np.searchsorted(available_orders, order)
+        if best_order_idx < len(available_orders):
+            best_order = available_orders[best_order_idx]
+
+            ## We have a file with the same order - just use it
+            if best_order == order:
+                return file_io.read_fits_image(
+                    paths.get_skymap_file_pointer(catalog_base_dir=catalog.catalog_base_dir, order=order)
+                )
+
+            ## We have a file with a greater order - downsample it
+            skymap = file_io.read_fits_image(
+                paths.get_skymap_file_pointer(catalog_base_dir=catalog.catalog_base_dir, order=best_order)
+            )
+            return skymap.reshape(hp.order2npix(order), -1).sum(axis=1)
+
+    if catalog.catalog_info.skymap_order:
+        if order is None or order == catalog.catalog_info.skymap_order:
+            return file_io.read_fits_image(
+                paths.get_skymap_file_pointer(catalog_base_dir=catalog.catalog_base_dir)
+            )
+        if order > catalog.catalog_info.skymap_order:
+            raise ValueError(
+                f"order should be less than stored skymap order ({catalog.catalog_info.skymap_order})"
+            )
+        skymap = file_io.read_fits_image(
+            paths.get_skymap_file_pointer(catalog_base_dir=catalog.catalog_base_dir)
+        )
+        return skymap.reshape(hp.order2npix(order), -1).sum(axis=1)
+
+    ## Deprecated - prefer reading skymap.fits to reading point_map.fits
+    map_file_pointer = paths.get_point_map_file_pointer(catalog.catalog_base_dir)
+    point_map = file_io.read_fits_image(map_file_pointer)
+    point_map_order = hp.npix2order(len(point_map))
+    if order is None or order == point_map_order:
+        return point_map
+    if point_map_order < order:
+        raise ValueError(f"order should be less than stored skymap order ({point_map_order})")
+    return point_map.reshape(hp.order2npix(order), -1).sum(axis=1)
+
+
+def write_skymap(histogram: np.ndarray, catalog_dir: str | Path | UPath, orders: list | int | None = None):
+    """Write the object spatial distribution information to a healpix SKYMAP FITS file.
+
+    Args:
+        histogram (:obj:`np.ndarray`): one-dimensional numpy array of long integers where the
+            value at each index corresponds to the number of objects found at the healpix pixel.
+        catalog_dir (path-like): base directory of the catalog in which to write the skymap file(s)
+        order (list): list of orders to write additional skymap files. if provided and not empty,
+            we will write a `skymap.K.fits` for each integer K in the list. if empty or None,
+            we will not write additional files.
+    """
+    catalog_dir = file_io.get_upath(catalog_dir)
+    map_file_pointer = paths.get_skymap_file_pointer(catalog_dir)
+    file_io.write_fits_image(histogram=histogram, map_file_pointer=map_file_pointer)
+    if orders:
+        original_order = hp.npix2order(len(histogram))
+        if not isinstance(orders, Iterable):
+            ## allow input of a single order.
+            orders = [orders]
+        for order in orders:
+            if order > original_order:
+                raise ValueError(
+                    f"sub-sampling skymap order should be less than overal order ({original_order})"
+                )
+            sampled_histogram = histogram.reshape(hp.order2npix(order), -1).sum(axis=1)
+            map_file_pointer = paths.get_skymap_file_pointer(catalog_dir, order=order)
+            file_io.write_fits_image(histogram=sampled_histogram, map_file_pointer=map_file_pointer)

tests/conftest.py

@@ -235,6 +235,11 @@ def small_sky_source_dir(test_data_dir) -> str:
     return test_data_dir / "small_sky_source"
 
 
+@pytest.fixture
+def small_sky_source_skymap_dir(test_data_dir) -> str:
+    return test_data_dir / "small_sky_source_skymap"
+
+
 @pytest.fixture
 def small_sky_collection_dir(test_data_dir) -> str:
     return test_data_dir / "small_sky_o1_collection"

tests/data/small_sky_source_skymap/dataset/_common_metadata

@@ -0,0 +1,15 @@
+Norder,Npix
+0,4
+1,47
+2,176
+2,177
+2,178
+2,179
+2,180
+2,181
+2,182
+2,183
+2,184
+2,185
+2,186
+2,187

tests/data/small_sky_source_skymap/dataset/_metadata

@@ -0,0 +1,16 @@
+#HATS catalog
+obs_collection=small_sky_source
+dataproduct_type=source
+hats_nrows=17161
+hats_col_ra=source_ra
+hats_col_dec=source_dec
+hats_max_rows=3000
+hats_order=2
+moc_sky_fraction=0.16667
+hats_builder=hats-import v0.4.6.dev1+gf00cd7a
+hats_creation_date=2025-03-05T16\:16UTC
+hats_estsize=99275
+hats_release_date=2024-09-18
+hats_version=v0.1
+hats_skymap_order=6
+hats_skymap_alt_orders=2 4