Repository: astronomy-commons/hats

Issue:
Implement coverage computation function for skymap
**Feature request**

Currently we have `hats.catalog.partition_info.PartitionInfo.calculate_fractional_coverage`, which gives fraction of sky covered by the partition tiles. It would also very useful to have a function which would give an estimate based on `skymap` files instead of partitions, e.g.
```python
def skymap_coverage(catalog, order=None):
    if order is None:
        order = order_of_skymap
    return np.mean(hats.io.skymap.read_skymap(catalog, order) > 0)
```

For most of the surveys (but maybe very sparse and small ones) this estimate will be closer to the actual coverage.

**Before submitting**
Please check the following:

- [x] I have described the purpose of the suggested change, specifying what I need the enhancement to accomplish, i.e. what problem it solves.
- [x] I have included any relevant links, screenshots, environment information, and data relevant to implementing the requested feature, as well as pseudocode for how I want to access the new functionality.
- [x] If I have ideas for how the new feature could be implemented, I have provided explanations and/or pseudocode and/or task lists for the steps.

Discussion:

Fix (patch):
diff --git a/src/hats/io/__init__.py b/src/hats/io/__init__.py
index fc24ca97..0d717d71 100644
--- a/src/hats/io/__init__.py
+++ b/src/hats/io/__init__.py
@@ -10,3 +10,4 @@
     pixel_catalog_file,
     pixel_directory,
 )
+from .skymap import skymap_coverage
diff --git a/src/hats/io/skymap.py b/src/hats/io/skymap.py
index acaa5532..431ff142 100644
--- a/src/hats/io/skymap.py
+++ b/src/hats/io/skymap.py
@@ -69,6 +69,28 @@ def read_skymap(catalog, order):
     return point_map.reshape(hp.order2npix(order), -1).sum(axis=1)
 
 
+def skymap_coverage(catalog, order=None):
+    """Compute the fractional sky coverage of a catalog from its healpix skymap file.
+
+    Parameters
+    ----------
+    catalog : Catalog
+        Catalog object corresponding to an on-disk catalog with skymap files.
+    order : int, optional
+        healpix order to read the skymap at. If None, the order of the default
+        skymap will be used.
+
+    Returns
+    -------
+    float
+        fractional sky coverage between 0.0 and 1.0.
+    """
+    if catalog.catalog_info.skymap_order is None and catalog.catalog_info.skymap_alt_orders is None:
+        raise ValueError("Catalog does not have skymap information. Cannot compute skymap-based coverage.")
+    skymap = read_skymap(catalog, order)
+    return float(np.mean(skymap > 0))
+
+
 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.
 


Test changes:
diff --git a/tests/hats/io/test_skymap.py b/tests/hats/io/test_skymap.py
index f402fedb..b8b48257 100644
--- a/tests/hats/io/test_skymap.py
+++ b/tests/hats/io/test_skymap.py
@@ -6,7 +6,7 @@
 import hats.pixel_math.healpix_shim as hp
 from hats import read_hats
 from hats.io.file_io import read_fits_image
-from hats.io.skymap import read_skymap, write_skymap
+from hats.io.skymap import read_skymap, skymap_coverage, write_skymap
 
 
 def test_write_skymap_roundtrip(tmp_path):
@@ -162,3 +162,62 @@ def test_read_noalt_skymap(small_sky_source_dir, mocker):
     ## Requesting waaay too-high order should error.
     with pytest.raises(ValueError, match="order should be less"):
         skymap = read_skymap(catalog, 13)
+
+
+def test_skymap_coverage_no_skymap():
+    """Test that skymap_coverage raises an error for catalogs without skymap files."""
+    catalog = type(
+        "MockCatalog",
+        (),
+        {
+            "catalog_base_dir": "/nonexistent",
+            "catalog_info": type(
+                "MockInfo",
+                (),
+                {
+                    "skymap_order": None,
+                    "skymap_alt_orders": None,
+                },
+            )(),
+        },
+    )()
+    with pytest.raises(ValueError, match="does not have skymap information"):
+        skymap_coverage(catalog)
+
+
+def test_skymap_coverage_synthetic(tmp_path):
+    """Test skymap_coverage with a known synthetic skymap."""
+    num_pixels = hp.order2npix(2)
+    histogram = np.zeros(num_pixels, dtype=np.int64)
+    histogram[: num_pixels // 2] = 10
+
+    write_skymap(histogram, tmp_path)
+
+    catalog = type(
+        "MockCatalog",
+        (),
+        {
+            "catalog_base_dir": tmp_path,
+            "catalog_info": type(
+                "MockInfo",
+                (),
+                {
+                    "skymap_order": 2,
+                    "skymap_alt_orders": None,
+                },
+            )(),
+        },
+    )()
+
+    coverage = skymap_coverage(catalog)
+    assert coverage == pytest.approx(0.5)
+
+
+def test_skymap_coverage_changes_with_order(small_sky_source_dir):
+    """Test that coverage estimate changes when different orders are requested."""
+    catalog = read_hats(small_sky_source_dir)
+
+    coverage_default = skymap_coverage(catalog)
+    coverage_order4 = skymap_coverage(catalog, order=4)
+
+    assert coverage_default != coverage_order4