From 081961abb2ece87e1ffaefff2bfabdfccab6a581 Mon Sep 17 00:00:00 2001 From: Melissa DeLucchi <113376043+delucchi-cmu@users.noreply.github.com> Date: Wed, 4 Dec 2024 09:22:50 -0500 Subject: [PATCH] Make point_map.fits plotting more friendly. (#439) * Make point_map.fits plotting more friendly. * Unit test coverage. * Use astropy for unit conversion --- src/hats/inspection/__init__.py | 2 +- src/hats/inspection/visualize_catalog.py | 36 +++++++++++++++---- src/hats/pixel_math/healpix_shim.py | 21 +++++------ .../hats/inspection/test_visualize_catalog.py | 31 +++++++++++++++- tests/hats/pixel_math/test_healpix_shim.py | 32 +++++++++++++++-- 5 files changed, 102 insertions(+), 20 deletions(-) diff --git a/src/hats/inspection/__init__.py b/src/hats/inspection/__init__.py index e1786986..dac60516 100644 --- a/src/hats/inspection/__init__.py +++ b/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 diff --git a/src/hats/inspection/visualize_catalog.py b/src/hats/inspection/visualize_catalog.py index cce75e75..f6291974 100644 --- a/src/hats/inspection/visualize_catalog.py +++ b/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: + 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) - 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), -1).sum(axis=1) + else: + order = map_order + if unit is None: + unit = u.deg * 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): diff --git a/src/hats/pixel_math/healpix_shim.py b/src/hats/pixel_math/healpix_shim.py index b5705a77..f7457f78 100644 --- a/src/hats/pixel_math/healpix_shim.py +++ b/src/hats/pixel_math/healpix_shim.py @@ -40,21 +40,22 @@ 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=u.rad) -> float: if arcmin: - return np.rad2deg(resol_rad) * 60 + unit = u.arcmin + unit = u.Unit(unit) - return resol_rad + return np.sqrt(order2pixarea(order, unit=unit * unit)) -def order2pixarea(order: int, degrees: bool = False) -> float: - npix = order2npix(order) - pix_area_rad = 4 * np.pi / npix +def order2pixarea(order: int, *, degrees: bool = False, unit=u.sr) -> float: if degrees: - return pix_area_rad * (180 / np.pi) * (180 / np.pi) - return pix_area_rad + unit = "deg**2" + unit = u.Unit(unit) + + npix = order2npix(order) + pix_area_rad = 4 * np.pi / npix * u.steradian + return pix_area_rad.to_value(unit) def radec2pix(order: int, ra: float, dec: float) -> np.ndarray[np.int64]: diff --git a/tests/hats/inspection/test_visualize_catalog.py b/tests/hats/inspection/test_visualize_catalog.py index 8cfdf1a1..cdaeb87a 100644 --- a/tests/hats/inspection/test_visualize_catalog.py +++ b/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,34 @@ 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_density_errors(small_sky_source_dir): + small_sky_source_catalog = read_hats(small_sky_source_dir) + with pytest.raises(ValueError, match="plotting order"): + plot_density(small_sky_source_catalog, order=13) + + with pytest.raises(ValueError, match="not convertible"): + plot_density(small_sky_source_catalog, unit="arcmin") + + with pytest.raises(ValueError, match="catalog required"): + plot_density(None) + + 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()) diff --git a/tests/hats/pixel_math/test_healpix_shim.py b/tests/hats/pixel_math/test_healpix_shim.py index ba4b39ac..c33115a3 100644 --- a/tests/hats/pixel_math/test_healpix_shim.py +++ b/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 @@ -104,13 +105,27 @@ def test_order2pixarea(): assert pix_area_test == pix_area_expected -def test_order2pixarea_degrees(): +def test_order2pixarea_units(): orders = [0, 1, 5, 10, 20, 29] npix = [12 * (4**order) for order in orders] pix_area_expected = [np.rad2deg(np.rad2deg(4 * np.pi / x)) for x in npix] pix_area_test = [hps.order2pixarea(order, degrees=True) for order in orders] assert pix_area_test == pix_area_expected + pix_area_expected = [np.rad2deg(np.rad2deg(4 * np.pi / x)) * 3600 for x in npix] + pix_area_test = [hps.order2pixarea(order, unit="arcmin2") for order in orders] + assert_allclose(pix_area_test, pix_area_expected) + + pix_area_test = [hps.order2pixarea(order, unit=u.arcmin * u.arcmin) for order in orders] + assert_allclose(pix_area_test, pix_area_expected) + + pix_area_expected = [np.rad2deg(np.rad2deg(4 * np.pi / x)) * 12960000 for x in npix] + pix_area_test = [hps.order2pixarea(order, unit="arcsec2") for order in orders] + assert_allclose(pix_area_test, pix_area_expected) + + with pytest.raises(ValueError, match="not convertible"): + hps.order2pixarea(10, unit="arcmin") + def test_order2resol(): orders = [0, 1, 5, 10, 20, 29] @@ -123,7 +138,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():