From 7add8f0c0440e9b45790779b5cae8a87252ac4ce Mon Sep 17 00:00:00 2001 From: Sean McGuire Date: Fri, 1 Nov 2024 14:38:54 -0400 Subject: [PATCH 1/4] plot lower order pixels with correct borders --- src/hats/inspection/visualize_catalog.py | 85 +++++++++++---- .../hats/inspection/test_visualize_catalog.py | 100 ++++++++++++++---- 2 files changed, 140 insertions(+), 45 deletions(-) diff --git a/src/hats/inspection/visualize_catalog.py b/src/hats/inspection/visualize_catalog.py index 227a9425..866c93af 100644 --- a/src/hats/inspection/visualize_catalog.py +++ b/src/hats/inspection/visualize_catalog.py @@ -24,7 +24,6 @@ from matplotlib.path import Path from mocpy import MOC, WCS from mocpy.moc.plot.culling_backfacing_cells import backface_culling -from mocpy.moc.plot.fill import compute_healpix_vertices from mocpy.moc.plot.utils import _set_wcs from hats.io import file_io, paths @@ -311,30 +310,26 @@ def cull_from_pixel_map(depth_ipix_d: Dict[int, Tuple[np.ndarray, np.ndarray]], for depth in range(min_depth, max_split_depth + 1): # for each depth, check if pixels are too large, or wrap around projection, and split into pixels at # higher order - if depth < 3: - too_large_ipix = ipixels - too_large_vals = vals - else: - ipix_lon, ipix_lat = cdshealpix.vertices(ipixels, depth) + ipix_lon, ipix_lat = cdshealpix.vertices(ipixels, depth) - ipix_lon = ipix_lon[:, [2, 3, 0, 1]] - ipix_lat = ipix_lat[:, [2, 3, 0, 1]] - ipix_vertices = SkyCoord(ipix_lon, ipix_lat, frame=ICRS()) + ipix_lon = ipix_lon[:, [2, 3, 0, 1]] + ipix_lat = ipix_lat[:, [2, 3, 0, 1]] + ipix_vertices = SkyCoord(ipix_lon, ipix_lat, frame=ICRS()) - # Projection on the given WCS - xp, yp = skycoord_to_pixel(coords=ipix_vertices, wcs=wcs) - _, _, frontface_id = backface_culling(xp, yp) + # Projection on the given WCS + xp, yp = skycoord_to_pixel(coords=ipix_vertices, wcs=wcs) + _, _, frontface_id = backface_culling(xp, yp) - # Get the pixels which are backfacing the projection - backfacing_ipix = ipixels[~frontface_id] # backfacing - backfacing_vals = vals[~frontface_id] - frontface_ipix = ipixels[frontface_id] - frontface_vals = vals[frontface_id] + # Get the pixels which are backfacing the projection + backfacing_ipix = ipixels[~frontface_id] # backfacing + backfacing_vals = vals[~frontface_id] + frontface_ipix = ipixels[frontface_id] + frontface_vals = vals[frontface_id] - ipix_d.update({depth: (frontface_ipix, frontface_vals)}) + ipix_d.update({depth: (frontface_ipix, frontface_vals)}) - too_large_ipix = backfacing_ipix - too_large_vals = backfacing_vals + too_large_ipix = backfacing_ipix + too_large_vals = backfacing_vals next_depth = depth + 1 @@ -358,6 +353,46 @@ def cull_from_pixel_map(depth_ipix_d: Dict[int, Tuple[np.ndarray, np.ndarray]], return ipix_d +def compute_healpix_vertices(depth, ipix, wcs, step=1): + """Compute HEALPix vertices. + + Parameters + ---------- + depth : int + The depth of the HEALPix cells. + ipix : `numpy.ndarray` + The HEALPix cell index given as a `np.uint64` numpy array. + wcs : `astropy.wcs.WCS` + A WCS projection + + Returns + ------- + path_vertices, codes : (`np.array`, `np.array`) + """ + + depth = int(depth) + + ipix_lon, ipix_lat = cdshealpix.vertices(ipix, depth, step=step) + indices = np.concatenate([np.arange(2 * step, 4 * step), np.arange(2 * step)]) + + ipix_lon = ipix_lon[:, indices] + ipix_lat = ipix_lat[:, indices] + ipix_boundaries = SkyCoord(ipix_lon, ipix_lat, frame=ICRS()) + # Projection on the given WCS + xp, yp = skycoord_to_pixel(ipix_boundaries, wcs=wcs) + + raw_cells = [np.vstack((xp[:, i], yp[:, i])).T for i in range(4 * step)] + + cells = np.hstack(raw_cells + [np.zeros((raw_cells[0].shape[0], 2))]) + + path_vertices = cells.reshape(((4 * step + 1) * raw_cells[0].shape[0], 2)) + single_code = np.array([Path.MOVETO] + [Path.LINETO] * (step * 4 - 1) + [Path.CLOSEPOLY]) + + codes = np.tile(single_code, raw_cells[0].shape[0]) + + return path_vertices, codes + + def plot_healpix_map( healpix_map: np.ndarray, *, @@ -527,9 +562,15 @@ def _plot_healpix_value_map(ipix, depth, values, ax, wcs, cmap="viridis", norm=N plt_paths = [] cum_vals = [] for d, (ip, vals) in culled_d.items(): - vertices, codes = compute_healpix_vertices(depth=d, ipix=ip, wcs=wcs) + step = 1 if d >= 3 else 2 ** (3 - d) + vertices, codes = compute_healpix_vertices(depth=d, ipix=ip, wcs=wcs, step=step) for i in range(len(ip)): - plt_paths.append(Path(vertices[5 * i : 5 * (i + 1)], codes[5 * i : 5 * (i + 1)])) + plt_paths.append( + Path( + vertices[(4 * step + 1) * i : (4 * step + 1) * (i + 1)], + codes[(4 * step + 1) * i : (4 * step + 1) * (i + 1)], + ) + ) cum_vals.append(vals) col = PathCollection(plt_paths, cmap=cmap, norm=norm, **kwargs) col.set_array(np.concatenate(cum_vals)) diff --git a/tests/hats/inspection/test_visualize_catalog.py b/tests/hats/inspection/test_visualize_catalog.py index 762dd4f8..3216fc71 100644 --- a/tests/hats/inspection/test_visualize_catalog.py +++ b/tests/hats/inspection/test_visualize_catalog.py @@ -1,20 +1,29 @@ from unittest.mock import MagicMock import astropy.units as u +import cdshealpix import matplotlib.pyplot as plt import numpy as np import pytest from astropy.coordinates import Angle, SkyCoord from astropy.visualization.wcsaxes.frame import EllipticalFrame, RectangularFrame +from cdshealpix.ring import vertices from matplotlib.colors import LogNorm, Normalize +from matplotlib.path import Path from matplotlib.pyplot import get_cmap from mocpy import MOC, WCS +from mocpy.moc.plot import fill from mocpy.moc.plot.culling_backfacing_cells import from_moc -from mocpy.moc.plot.fill import compute_healpix_vertices from mocpy.moc.plot.utils import build_plotting_moc from hats.inspection import plot_pixels -from hats.inspection.visualize_catalog import cull_from_pixel_map, cull_to_fov, plot_healpix_map, plot_moc +from hats.inspection.visualize_catalog import ( + cull_from_pixel_map, + cull_to_fov, + plot_healpix_map, + compute_healpix_vertices, + plot_moc, +) # pylint: disable=no-member @@ -27,6 +36,48 @@ DEFAULT_PROJECTION = "MOL" +def test_healpix_vertices(): + depth = 3 + ipix = np.array([10, 11]) + fig = plt.figure() + wcs = WCS( + fig, + fov=DEFAULT_FOV, + center=DEFAULT_CENTER, + coordsys=DEFAULT_COORDSYS, + rotation=DEFAULT_ROTATION, + projection=DEFAULT_PROJECTION, + ).w + paths, codes = compute_healpix_vertices(depth, ipix, wcs) + mocpy_paths, mocpy_codes = fill.compute_healpix_vertices(depth, ipix, wcs) + np.testing.assert_array_equal(paths, mocpy_paths) + np.testing.assert_array_equal(codes, mocpy_codes) + + +def test_healpix_vertices_step(): + depth = 1 + ipix = np.array([10, 11]) + fig = plt.figure() + step = 4 + wcs = WCS( + fig, + fov=DEFAULT_FOV, + center=DEFAULT_CENTER, + coordsys=DEFAULT_COORDSYS, + rotation=DEFAULT_ROTATION, + projection=DEFAULT_PROJECTION, + ).w + paths, codes = compute_healpix_vertices(depth, ipix, wcs, step=step) + np.testing.assert_array_equal( + codes, np.tile(np.array([Path.MOVETO] + [Path.LINETO] * (step * 4 - 1) + [Path.CLOSEPOLY]), len(ipix)) + ) + mocpy_paths, _ = fill.compute_healpix_vertices(depth, ipix, wcs) + first_path_vertex_indices = np.array([0, step, 2 * step, 3 * step, 4 * step]) + start_path_index = np.array(([0] * 5) + ([first_path_vertex_indices[-1] + 1] * 5)) + vertex_indices = start_path_index + np.tile(first_path_vertex_indices, len(ipix)) + np.testing.assert_array_equal(paths[vertex_indices], mocpy_paths) + + def test_plot_healpix_pixels(): order = 3 length = 10 @@ -89,7 +140,7 @@ def test_plot_healpix_pixels_different_order(): np.testing.assert_array_equal(col.get_array(), pix_map) -def test_order_0_pixels_split_to_order_3(): +def test_order_0_pixel_plots_with_step(): map_value = 0.5 order_0_pix = 4 ipix = np.array([order_0_pix]) @@ -99,8 +150,7 @@ def test_order_0_pixels_split_to_order_3(): assert len(ax.collections) == 1 col = ax.collections[0] paths = col.get_paths() - length = 4**3 - order3_ipix = np.arange(length * order_0_pix, length * (order_0_pix + 1)) + length = 1 assert len(paths) == length wcs = WCS( fig, @@ -110,11 +160,10 @@ def test_order_0_pixels_split_to_order_3(): rotation=DEFAULT_ROTATION, projection=DEFAULT_PROJECTION, ).w - all_verts, all_codes = compute_healpix_vertices(3, order3_ipix, wcs) - for i, (path, ipix) in enumerate(zip(paths, order3_ipix)): - verts, codes = all_verts[i * 5 : (i + 1) * 5], all_codes[i * 5 : (i + 1) * 5] - np.testing.assert_array_equal(path.vertices, verts) - np.testing.assert_array_equal(path.codes, codes) + all_verts, all_codes = compute_healpix_vertices(0, ipix, wcs, step=2**3) + assert len(all_verts) == len(ipix) * 4 * (2**3) + 1 + np.testing.assert_array_equal(paths[0].vertices, all_verts) + np.testing.assert_array_equal(paths[0].codes, all_codes) np.testing.assert_array_equal(col.get_array(), np.full(length, fill_value=map_value)) @@ -130,9 +179,9 @@ def test_edge_pixels_split_to_order_7(): for iter_ord in range(1, 8): edge_pixels[iter_ord] = [p * 4 + i for p in edge_pixels[iter_ord - 1] for i in (2, 3)] non_edge_pixels = {} - pixels_ord3 = np.arange(4**3 * order_0_pix, 4**3 * (order_0_pix + 1)) - non_edge_pixels[3] = pixels_ord3[~np.isin(pixels_ord3, edge_pixels[3])] - for iter_ord in range(4, 8): + pixels_ord1 = np.arange(4 * order_0_pix, 4 * (order_0_pix + 1)) + non_edge_pixels[1] = pixels_ord1[~np.isin(pixels_ord1, edge_pixels[1])] + for iter_ord in range(2, 8): pixels_ord = np.concatenate([np.arange(4 * pix, 4 * (pix + 1)) for pix in edge_pixels[iter_ord - 1]]) non_edge_pixels[iter_ord] = pixels_ord[~np.isin(pixels_ord, edge_pixels[iter_ord])] col = ax.collections[0] @@ -150,14 +199,15 @@ def test_edge_pixels_split_to_order_7(): ords = np.concatenate([np.full(len(x), fill_value=o) for o, x in non_edge_pixels.items()]) pixels = np.concatenate([np.array(x) for _, x in non_edge_pixels.items()]) for path, iter_ord, pix in zip(paths, ords, pixels): - verts, codes = compute_healpix_vertices(iter_ord, np.array([pix]), wcs) + step = 1 if iter_ord >= 3 else 2 ** (3 - iter_ord) + verts, codes = compute_healpix_vertices(iter_ord, np.array([pix]), wcs, step=step) np.testing.assert_array_equal(path.vertices, verts) np.testing.assert_array_equal(path.codes, codes) np.testing.assert_array_equal(col.get_array(), np.full(length, fill_value=map_value)) def test_cull_from_pixel_map(): - order = 1 + order = 3 ipix = np.arange(12 * 4**order) pix_map = np.arange(12 * 4**order) map_dict = {order: (ipix, pix_map)} @@ -309,9 +359,14 @@ def test_plot_healpix_map(): all_vals = [] start_i = 0 for iter_ord, (pixels, pix_map) in culled_dict.items(): - all_verts, all_codes = compute_healpix_vertices(iter_ord, pixels, wcs) + step = 1 if iter_ord >= 3 else 2 ** (3 - iter_ord) + vert_len = step * 4 + 1 + all_verts, all_codes = compute_healpix_vertices(iter_ord, pixels, wcs, step=step) for i, _ in enumerate(pixels): - verts, codes = all_verts[i * 5 : (i + 1) * 5], all_codes[i * 5 : (i + 1) * 5] + verts, codes = ( + all_verts[i * vert_len : (i + 1) * vert_len], + all_codes[i * vert_len : (i + 1) * vert_len], + ) path = paths[start_i + i] np.testing.assert_array_equal(path.vertices, verts) np.testing.assert_array_equal(path.codes, codes) @@ -652,11 +707,9 @@ def test_plot_kwargs(): 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()) - order_3_pixels = [p for pix in pixels for p in pix.convert_to_higher_order(3 - pix.order)] - order_3_orders = [pix.order for pix in pixels for _ in pix.convert_to_higher_order(3 - pix.order)] col = ax.collections[0] paths = col.get_paths() - assert len(paths) == len(order_3_pixels) + assert len(paths) == len(pixels) wcs = WCS( fig, fov=DEFAULT_FOV, @@ -665,11 +718,12 @@ def test_catalog_plot(small_sky_order1_catalog): rotation=DEFAULT_ROTATION, projection=DEFAULT_PROJECTION, ).w - for p, path in zip(order_3_pixels, paths): - verts, codes = compute_healpix_vertices(p.order, np.array([p.pixel]), wcs) + for p, path in zip(pixels, paths): + step = 2 ** (3 - p.order) + verts, codes = compute_healpix_vertices(p.order, np.array([p.pixel]), wcs, step=step) np.testing.assert_array_equal(path.vertices, verts) np.testing.assert_array_equal(path.codes, codes) - np.testing.assert_array_equal(col.get_array(), np.array(order_3_orders)) + np.testing.assert_array_equal(col.get_array(), np.array([p.order for p in pixels])) assert ax.get_title() == f"Catalog pixel density map - {small_sky_order1_catalog.catalog_name}" From bbaa15f652479d734ab9e068e53be0f61325dfb8 Mon Sep 17 00:00:00 2001 From: Sean McGuire Date: Fri, 1 Nov 2024 14:59:23 -0400 Subject: [PATCH 2/4] isort --- tests/hats/inspection/test_visualize_catalog.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/tests/hats/inspection/test_visualize_catalog.py b/tests/hats/inspection/test_visualize_catalog.py index 3216fc71..5effe5b2 100644 --- a/tests/hats/inspection/test_visualize_catalog.py +++ b/tests/hats/inspection/test_visualize_catalog.py @@ -1,13 +1,11 @@ from unittest.mock import MagicMock import astropy.units as u -import cdshealpix import matplotlib.pyplot as plt import numpy as np import pytest from astropy.coordinates import Angle, SkyCoord from astropy.visualization.wcsaxes.frame import EllipticalFrame, RectangularFrame -from cdshealpix.ring import vertices from matplotlib.colors import LogNorm, Normalize from matplotlib.path import Path from matplotlib.pyplot import get_cmap @@ -18,10 +16,10 @@ from hats.inspection import plot_pixels from hats.inspection.visualize_catalog import ( + compute_healpix_vertices, cull_from_pixel_map, cull_to_fov, plot_healpix_map, - compute_healpix_vertices, plot_moc, ) From b601f237295702cc4d46d3622a72c62db8d7e721 Mon Sep 17 00:00:00 2001 From: Sean McGuire Date: Fri, 1 Nov 2024 16:26:43 -0400 Subject: [PATCH 3/4] add creditation --- src/hats/inspection/visualize_catalog.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/hats/inspection/visualize_catalog.py b/src/hats/inspection/visualize_catalog.py index 866c93af..06344c99 100644 --- a/src/hats/inspection/visualize_catalog.py +++ b/src/hats/inspection/visualize_catalog.py @@ -356,6 +356,8 @@ def cull_from_pixel_map(depth_ipix_d: Dict[int, Tuple[np.ndarray, np.ndarray]], def compute_healpix_vertices(depth, ipix, wcs, step=1): """Compute HEALPix vertices. + Modified from mocpy.moc.plot.fill.compute_healpix_vertices + Parameters ---------- depth : int From 1ff39ad28818160c5cc852ad73255a5dcfb004af Mon Sep 17 00:00:00 2001 From: Sean McGuire Date: Fri, 1 Nov 2024 17:04:31 -0400 Subject: [PATCH 4/4] add comments --- tests/hats/inspection/test_visualize_catalog.py | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/tests/hats/inspection/test_visualize_catalog.py b/tests/hats/inspection/test_visualize_catalog.py index 5effe5b2..397c8162 100644 --- a/tests/hats/inspection/test_visualize_catalog.py +++ b/tests/hats/inspection/test_visualize_catalog.py @@ -66,10 +66,13 @@ def test_healpix_vertices_step(): projection=DEFAULT_PROJECTION, ).w paths, codes = compute_healpix_vertices(depth, ipix, wcs, step=step) + # checks the codes match the expected matplotlib path codes np.testing.assert_array_equal( codes, np.tile(np.array([Path.MOVETO] + [Path.LINETO] * (step * 4 - 1) + [Path.CLOSEPOLY]), len(ipix)) ) mocpy_paths, _ = fill.compute_healpix_vertices(depth, ipix, wcs) + # mocpy only generates path points at the healpix corner vertices. So we subsample our generated vertices + # to check the corners match the expected mocpy generated ones first_path_vertex_indices = np.array([0, step, 2 * step, 3 * step, 4 * step]) start_path_index = np.array(([0] * 5) + ([first_path_vertex_indices[-1] + 1] * 5)) vertex_indices = start_path_index + np.tile(first_path_vertex_indices, len(ipix)) @@ -159,6 +162,8 @@ def test_order_0_pixel_plots_with_step(): projection=DEFAULT_PROJECTION, ).w all_verts, all_codes = compute_healpix_vertices(0, ipix, wcs, step=2**3) + # assert the number of vertices == number of pixels * 4 sides per pixel * steps per side + 1 for the + # close polygon assert len(all_verts) == len(ipix) * 4 * (2**3) + 1 np.testing.assert_array_equal(paths[0].vertices, all_verts) np.testing.assert_array_equal(paths[0].codes, all_codes) @@ -173,9 +178,15 @@ def test_edge_pixels_split_to_order_7(): depth = np.array([0]) fig, ax = plot_healpix_map(pix_map, ipix=ipix, depth=depth) assert len(ax.collections) == 1 + + # Generate a dictionary of pixel indices that have sides that align with the meridian at ra = 180deg, the + # right edge of the plot edge_pixels = {0: [order_0_pix]} for iter_ord in range(1, 8): edge_pixels[iter_ord] = [p * 4 + i for p in edge_pixels[iter_ord - 1] for i in (2, 3)] + + # Generate a dictionary of subpixels of the order 0 pixel that are not on the edge of the plot, i.e. the + # pixels that should be in the culled plot non_edge_pixels = {} pixels_ord1 = np.arange(4 * order_0_pix, 4 * (order_0_pix + 1)) non_edge_pixels[1] = pixels_ord1[~np.isin(pixels_ord1, edge_pixels[1])] @@ -184,6 +195,8 @@ def test_edge_pixels_split_to_order_7(): non_edge_pixels[iter_ord] = pixels_ord[~np.isin(pixels_ord, edge_pixels[iter_ord])] col = ax.collections[0] paths = col.get_paths() + + # Check that the plotted paths match the non_edge_pixels generated length = sum(len(x) for x in non_edge_pixels.values()) assert len(paths) == length wcs = WCS(