Make point_map.fits plotting more friendly.

Author: delucchi-cmu

src/hats/inspection/__init__.py

@@ -1 +1 @@
-from .visualize_catalog import plot_pixel_list, plot_pixels, plot_points
+from .visualize_catalog import plot_density, plot_pixel_list, plot_pixels

src/hats/inspection/visualize_catalog.py

@@ -26,6 +26,7 @@
 from mocpy.moc.plot.culling_backfacing_cells import backface_culling
 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.pixel_math import HealpixPixel
 from hats.pixel_tree.moc_filter import perform_filter_by_moc
@@ -49,19 +50,42 @@ def _read_point_map(catalog_base_dir):
     return file_io.read_fits_image(map_file_pointer)
 
 
-def plot_points(catalog: Catalog, plot_title: str | None = None, **kwargs):
-    """Create a visual map of the input points of an in-memory catalog.
-
+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
+        order (int): Optionally reduce the display healpix order, and aggregate smaller tiles.
         kwargs: Additional args to pass to `plot_healpix_map`
     """
     if not catalog.on_disk:
         raise ValueError("on disk catalog required for point-wise visualization")
     point_map = _read_point_map(catalog.catalog_base_dir)
-    default_title = f"Catalog point density map - {catalog.catalog_name}"
-    return plot_healpix_map(point_map, title=default_title if plot_title is None else plot_title, **kwargs)
+    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), 2 ** (2 * (map_order - order)))).sum(axis=1)
+    else:
+        order = map_order
+    if unit is None:
+        unit = u.deg
+
+    pix_area = hp.order2pixarea(order, unit=unit)
+
+    point_map = point_map / pix_area
+    default_title = f"Angular density of catalog {catalog.catalog_name}"
+    fig, ax = plot_healpix_map(
+        point_map, title=default_title if plot_title is None else plot_title, cbar=False, **kwargs
+    )
+    col = ax.collections[0]
+    plt.colorbar(
+        col,
+        label=f"count / {unit} sq",
+    )
+    return fig, ax
 
 
 def plot_pixels(catalog: HealpixDataset, plot_title: str | None = None, **kwargs):

src/hats/pixel_math/healpix_shim.py

@@ -40,21 +40,45 @@ def order2npix(order: int) -> int:
     return 12 * (1 << (2 * order))
 
 
-def order2resol(order: int, arcmin: bool = False) -> float:
-    resol_rad = np.sqrt(order2pixarea(order))
-
+def order2resol(order: int, *, arcmin: bool = False, unit=None) -> float:
     if arcmin:
-        return np.rad2deg(resol_rad) * 60
+        unit = u.arcmin
+    else:
+        unit = _convert_unit(unit, default=u.rad)
+
+    return np.sqrt(order2pixarea(order, unit=unit))
 
-    return resol_rad
 
+def order2pixarea(order: int, *, degrees: bool = False, unit=None) -> float:
+    if degrees:
+        unit = u.deg
+    else:
+        unit = _convert_unit(unit, default=u.rad)
 
-def order2pixarea(order: int, degrees: bool = False) -> float:
     npix = order2npix(order)
     pix_area_rad = 4 * np.pi / npix
-    if degrees:
+
+    # Yes, we're using astropy units, but those are not for SOLID angles.
+    if unit == u.rad:
+        return pix_area_rad
+    if unit == u.deg:
         return pix_area_rad * (180 / np.pi) * (180 / np.pi)
-    return pix_area_rad
+    if unit == u.arcmin:
+        return pix_area_rad * (180 / np.pi) * (180 / np.pi) * 3600
+    if unit == u.arcsec:
+        return pix_area_rad * (180 / np.pi) * (180 / np.pi) * 12960000
+    raise ValueError("Unit must be angular measurement.")
+
+
+def _convert_unit(unit, default=u.rad):
+    if unit is not None:
+        unit = u.Unit(unit)
+    else:
+        unit = default
+
+    if unit not in (u.arcmin, u.arcsec, u.deg, u.rad):
+        raise ValueError("Unit must be angular measurement.")
+    return unit
 
 
 def radec2pix(order: int, ra: float, dec: float) -> np.ndarray[np.int64]:

tests/hats/inspection/test_visualize_catalog.py

@@ -15,7 +15,8 @@
 from mocpy.moc.plot.culling_backfacing_cells import from_moc
 from mocpy.moc.plot.utils import build_plotting_moc
 
-from hats.inspection import plot_pixels
+from hats import read_hats
+from hats.inspection import plot_density, plot_pixels
 from hats.inspection.visualize_catalog import (
     compute_healpix_vertices,
     cull_from_pixel_map,
@@ -831,6 +832,22 @@ def test_plot_existing_wrong_axes():
     np.testing.assert_array_equal(col.get_array(), pix_map)
 
 
+def test_catalog_plot_density(small_sky_source_dir):
+    """Test plotting pixel-density for on-disk catalog.
+    Confirm plotting at lower order doesn't have a warning, and creates fewer plot paths."""
+    small_sky_source_catalog = read_hats(small_sky_source_dir)
+    with pytest.warns(match="smaller"):
+        _, ax = plot_density(small_sky_source_catalog)
+    order10_paths = ax.collections[0].get_paths()
+    assert "Angular density of catalog small_sky_source" == ax.get_title()
+
+    _, ax = plot_density(small_sky_source_catalog, order=3)
+    order3_paths = ax.collections[-1].get_paths()
+    assert "Angular density of catalog small_sky_source" == ax.get_title()
+
+    assert len(order3_paths) < len(order10_paths)
+
+
 def test_catalog_plot(small_sky_order1_catalog):
     fig, ax = plot_pixels(small_sky_order1_catalog)
     pixels = sorted(small_sky_order1_catalog.get_healpix_pixels())

tests/hats/pixel_math/test_healpix_shim.py

@@ -1,3 +1,4 @@
+import astropy.units as u
 import cdshealpix
 import numpy as np
 import pytest
@@ -112,6 +113,17 @@ def test_order2pixarea_degrees():
     assert pix_area_test == pix_area_expected
 
 
+def test_order2pixarea_arcmin():
+    orders = [20, 29]
+    npix = [12 * (4**order) for order in orders]
+    pix_area_expected = [np.rad2deg(np.rad2deg(4 * np.pi / x)) * 3600 for x in npix]
+    pix_area_test = [hps.order2pixarea(order, unit=u.arcmin) for order in orders]
+    assert_allclose(pix_area_test, pix_area_expected)
+
+    pix_area_test = [hps.order2pixarea(order, unit=u.arcminute) for order in orders]
+    assert_allclose(pix_area_test, pix_area_expected)
+
+
 def test_order2resol():
     orders = [0, 1, 5, 10, 20, 29]
     resol_expected = [np.sqrt(hps.order2pixarea(order)) for order in orders]
@@ -123,7 +135,20 @@ def test_order2resol_arcmin():
     orders = [0, 1, 5, 10, 20, 29]
     resol_expected = [np.rad2deg(np.sqrt(hps.order2pixarea(order))) * 60 for order in orders]
     resol_test = [hps.order2resol(order, arcmin=True) for order in orders]
-    assert resol_test == resol_expected
+    assert_allclose(resol_test, resol_expected)
+
+
+def test_order2resol_degree():
+    orders = [0, 1, 5, 10, 20, 29]
+    resol_expected = [np.rad2deg(np.sqrt(hps.order2pixarea(order))) for order in orders]
+    resol_test = [hps.order2resol(order, unit=u.deg) for order in orders]
+    assert_allclose(resol_test, resol_expected)
+
+    resol_test = [hps.order2resol(order, unit=u.degree) for order in orders]
+    assert_allclose(resol_test, resol_expected)
+
+    resol_test = [hps.order2resol(order, unit="deg") for order in orders]
+    assert_allclose(resol_test, resol_expected)
 
 
 def test_radec2pix_lonlat():