Merge pull request #282 from astronomy-commons/main

getting most recent version to test changes

src/.pylintrc

@@ -62,7 +62,7 @@ ignore-patterns=_version.py
 # (useful for modules/projects where namespaces are manipulated during runtime
 # and thus existing member attributes cannot be deduced by static analysis). It
 # supports qualified module names, as well as Unix pattern matching.
-ignored-modules=
+ignored-modules=astropy.units
 
 # Python code to execute, usually for sys.path manipulation such as
 # pygtk.require().

src/hipscat/catalog/catalog.py

@@ -13,13 +13,9 @@
 from hipscat.catalog.catalog_type import CatalogType
 from hipscat.catalog.healpix_dataset.healpix_dataset import HealpixDataset, PixelInputTypes
 from hipscat.pixel_math import HealpixPixel
-from hipscat.pixel_math.box_filter import filter_pixels_by_box, wrap_ra_angles
-from hipscat.pixel_math.cone_filter import filter_pixels_by_cone
-from hipscat.pixel_math.polygon_filter import (
-    CartesianCoordinates,
-    SphericalCoordinates,
-    filter_pixels_by_polygon,
-)
+from hipscat.pixel_math.box_filter import generate_box_moc, wrap_ra_angles
+from hipscat.pixel_math.cone_filter import generate_cone_moc
+from hipscat.pixel_math.polygon_filter import CartesianCoordinates, SphericalCoordinates, generate_polygon_moc
 from hipscat.pixel_math.validators import (
     validate_box_search,
     validate_declination_values,
@@ -85,7 +81,7 @@ def filter_by_cone(self, ra: float, dec: float, radius_arcsec: float) -> Catalog
         """
         validate_radius(radius_arcsec)
         validate_declination_values(dec)
-        return self.filter_from_pixel_list(filter_pixels_by_cone(self.pixel_tree, ra, dec, radius_arcsec))
+        return self.filter_by_moc(generate_cone_moc(ra, dec, radius_arcsec, self.get_max_coverage_order()))
 
     def filter_by_box(
         self, ra: Tuple[float, float] | None = None, dec: Tuple[float, float] | None = None
@@ -103,7 +99,7 @@ def filter_by_box(
         """
         ra = tuple(wrap_ra_angles(ra)) if ra else None
         validate_box_search(ra, dec)
-        return self.filter_from_pixel_list(filter_pixels_by_box(self.pixel_tree, ra, dec))
+        return self.filter_by_moc(generate_box_moc(ra, dec, self.get_max_coverage_order()))
 
     def filter_by_polygon(self, vertices: List[SphericalCoordinates] | List[CartesianCoordinates]) -> Catalog:
         """Filter the pixels in the catalog to only include the pixels that overlap
@@ -122,9 +118,11 @@ def filter_by_polygon(self, vertices: List[SphericalCoordinates] | List[Cartesia
             validate_declination_values(dec)
             # Get the coordinates vector on the unit sphere if we were provided
             # with polygon spherical coordinates of ra and dec
-            vertices = hp.ang2vec(ra, dec, lonlat=True)
-        validate_polygon(vertices)
-        return self.filter_from_pixel_list(filter_pixels_by_polygon(self.pixel_tree, vertices))
+            cart_vertices = hp.ang2vec(ra, dec, lonlat=True)
+        else:
+            cart_vertices = vertices
+        validate_polygon(cart_vertices)
+        return self.filter_by_moc(generate_polygon_moc(cart_vertices, self.get_max_coverage_order()))
 
     def generate_negative_tree_pixels(self) -> List[HealpixPixel]:
         """Get the leaf nodes at each healpix order that have zero catalog data.

src/hipscat/catalog/healpix_dataset/healpix_dataset.py

@@ -14,6 +14,7 @@
 from hipscat.io import FilePointer, file_io, paths
 from hipscat.pixel_math import HealpixPixel
 from hipscat.pixel_tree import PixelAlignment, PixelAlignmentType
+from hipscat.pixel_tree.moc_filter import filter_by_moc
 from hipscat.pixel_tree.pixel_alignment import align_with_mocs
 from hipscat.pixel_tree.pixel_tree import PixelTree
 
@@ -131,18 +132,45 @@ def _check_files_exist(cls, catalog_base_dir: FilePointer, storage_options: dict
                 f"_metadata or partition info file is required in catalog directory {catalog_base_dir}"
             )
 
+    def get_max_coverage_order(self) -> int:
+        """Gets the maximum HEALPix order for which the coverage of the catalog is known from the pixel
+        tree and moc if it exists"""
+        max_order = (
+            max(self.moc.max_order, self.pixel_tree.get_max_depth())
+            if self.moc is not None
+            else self.pixel_tree.get_max_depth()
+        )
+        return max_order
+
     def filter_from_pixel_list(self, pixels: List[HealpixPixel]) -> Self:
-        """Filter the pixels in the catalog to only include the requested pixels.
+        """Filter the pixels in the catalog to only include any that overlap with the requested pixels.
 
         Args:
             pixels (List[HealpixPixels]): the pixels to include
 
         Returns:
-            A new catalog with only those pixels. Note that we reset the total_rows
-            to None, instead of performing a scan over the new pixel sizes.
+            A new catalog with only the pixels that overlap with the given pixels. Note that we reset the
+            total_rows to None, as updating would require a scan over the new pixel sizes.
         """
+        orders = np.array([p.order for p in pixels])
+        pixel_inds = np.array([p.pixel for p in pixels])
+        max_order = np.max(orders) if len(orders) > 0 else 0
+        moc = MOC.from_healpix_cells(ipix=pixel_inds, depth=orders, max_depth=max_order)
+        return self.filter_by_moc(moc)
+
+    def filter_by_moc(self, moc: MOC) -> Self:
+        """Filter the pixels in the catalog to only include the pixels that overlap with the moc provided.
+
+        Args:
+            moc (mocpy.MOC): the moc to filter by
+
+        Returns:
+            A new catalog with only the pixels that overlap with the moc. Note that we reset the total_rows
+            to None, as updating would require a scan over the new pixel sizes."""
+        filtered_tree = filter_by_moc(self.pixel_tree, moc)
+        filtered_moc = self.moc.intersection(moc) if self.moc is not None else None
         filtered_catalog_info = dataclasses.replace(self.catalog_info, total_rows=None)
-        return self.__class__(filtered_catalog_info, pixels)
+        return self.__class__(filtered_catalog_info, filtered_tree, moc=filtered_moc)
 
     def align(
         self, other_cat: Self, alignment_type: PixelAlignmentType = PixelAlignmentType.INNER

src/hipscat/catalog/margin_cache/margin_catalog.py

@@ -1,11 +1,13 @@
 from __future__ import annotations
 
+import healpy as hp
 from mocpy import MOC
-from typing_extensions import TypeAlias
+from typing_extensions import Self, TypeAlias
 
 from hipscat.catalog.catalog_type import CatalogType
 from hipscat.catalog.healpix_dataset.healpix_dataset import HealpixDataset, PixelInputTypes
 from hipscat.catalog.margin_cache import MarginCacheCatalogInfo
+from hipscat.pixel_tree.moc_utils import copy_moc
 
 
 class MarginCatalog(HealpixDataset):
@@ -46,3 +48,27 @@ def __init__(
         super().__init__(
             catalog_info, pixels, catalog_path=catalog_path, moc=moc, storage_options=storage_options
         )
+
+    def filter_by_moc(self, moc: MOC) -> Self:
+        """Filter the pixels in the margin catalog to only include the margin pixels that overlap with the moc
+
+        For the case of margin pixels, this includes any pixels whose margin areas may overlap with the moc.
+        This is not always done with a high accuracy, but always includes any pixels that will overlap,
+        and may include extra partitions that do not.
+
+        Args:
+            moc (mocpy.MOC): the moc to filter by
+
+        Returns:
+            A new margin catalog with only the pixels that overlap or that have margin area that overlap with
+            the moc. Note that we reset the total_rows to None, as updating would require a scan over the new
+            pixel sizes."""
+        max_order = moc.max_order
+        max_order_size = hp.nside2resol(2**max_order, arcmin=True)
+        if self.catalog_info.margin_threshold > max_order_size * 60:
+            raise ValueError(
+                f"Cannot Filter Margin: Margin size {self.catalog_info.margin_threshold} is "
+                f"greater than the size of a pixel at the highest order {max_order}."
+            )
+        search_moc = copy_moc(moc).add_neighbours()
+        return super().filter_by_moc(search_moc)

src/hipscat/io/file_io/file_io.py

@@ -2,6 +2,7 @@
 
 import json
 import tempfile
+import warnings
 from typing import Any, Dict, Tuple, Union
 
 import healpy as hp
@@ -278,8 +279,17 @@ def read_fits_image(map_file_pointer: FilePointer, storage_options: Union[Dict[A
         with file_system.open(map_file_pointer, "rb") as _map_file:
             map_data = _map_file.read()
             _tmp_file.write(map_data)
-            map_fits_image = hp.read_map(_tmp_file.name)
-    return map_fits_image
+            map_fits_image = hp.read_map(_tmp_file.name, nest=True, h=True)
+            header_dict = dict(map_fits_image[1])
+            if header_dict["ORDERING"] != "NESTED":
+                warnings.warn(
+                    "point_map.fits file written in RING ordering, due to "
+                    "https://github.com/astronomy-commons/hipscat/issues/271. "
+                    "Converting to NESTED."
+                )
+                map_fits_image = hp.read_map(_tmp_file.name)
+                return map_fits_image
+            return map_fits_image[0]
 
 
 def write_fits_image(
@@ -296,7 +306,7 @@ def write_fits_image(
     file_system, map_file_pointer = get_fs(file_pointer=map_file_pointer, storage_options=storage_options)
     with tempfile.NamedTemporaryFile() as _tmp_file:
         with file_system.open(map_file_pointer, "wb") as _map_file:
-            hp.write_map(_tmp_file.name, histogram, overwrite=True, dtype=np.int64)
+            hp.write_map(_tmp_file.name, histogram, overwrite=True, dtype=np.int64, nest=True)
             _map_file.write(_tmp_file.read())
 
 

src/hipscat/io/write_metadata.py

@@ -12,6 +12,7 @@
 
 from hipscat.io import file_io, paths
 from hipscat.io.parquet_metadata import write_parquet_metadata as wpm
+from hipscat.pixel_math.healpix_pixel import HealpixPixel
 
 
 class HipscatEncoder(json.JSONEncoder):
@@ -26,6 +27,12 @@ def default(self, o):
             return str(o)
         if isinstance(o, (type, np.dtype, pd.core.dtypes.base.ExtensionDtype)):
             return str(o)
+        if isinstance(o, HealpixPixel):
+            return str(o)
+        if np.issubdtype(o.dtype, np.integer):
+            return int(o)
+        if isinstance(o, np.ndarray):
+            return o.tolist()
         return super().default(o)  # pragma: no cover
 
 

src/hipscat/pixel_math/box_filter.py

@@ -4,48 +4,34 @@
 
 import healpy as hp
 import numpy as np
+from mocpy import MOC
 
-from hipscat.pixel_math import HealpixPixel
-from hipscat.pixel_math.filter import get_filtered_pixel_list
 from hipscat.pixel_math.polygon_filter import SphericalCoordinates
-from hipscat.pixel_tree import PixelAlignmentType, align_trees
-from hipscat.pixel_tree.pixel_tree import PixelTree
 
 
-def filter_pixels_by_box(
-    pixel_tree: PixelTree, ra: Tuple[float, float] | None, dec: Tuple[float, float] | None
-) -> List[HealpixPixel]:
-    """Filter the leaf pixels in a pixel tree to return a partition_info dataframe
-    with the pixels that overlap with a right ascension and/or declination region.
+def generate_box_moc(ra: Tuple[float, float] | None, dec: Tuple[float, float] | None, order: int) -> MOC:
+    """Generates a MOC object that covers the specified box area
 
     Args:
-        pixel_tree (PixelTree): The catalog tree to filter pixels from
         ra (Tuple[float, float]): Right ascension range, in [0,360] degrees
         dec (Tuple[float, float]): Declination range, in [-90,90] degrees
+        order (int): Maximum order of the moc to generate the box at
 
     Returns:
-        List of HEALPix representing only the pixels that overlap with the right
-        ascension and/or declination region.
+        a MOC object that covers the specified box
     """
-    max_order = pixel_tree.get_max_depth()
-
-    filter_tree = None
-    ra_search_tree, dec_search_tree = None, None
+    filter_moc = None
+    ra_search_moc, dec_search_moc = None, None
 
     if ra is not None:
-        ra_search_tree = _generate_ra_strip_pixel_tree(ra, max_order)
-        filter_tree = ra_search_tree
+        ra_search_moc = _generate_ra_strip_moc(ra, order)
+        filter_moc = ra_search_moc
     if dec is not None:
-        dec_search_tree = _generate_dec_strip_pixel_tree(dec, max_order)
-        filter_tree = dec_search_tree
-    if ra_search_tree is not None and dec_search_tree is not None:
-        filter_tree = align_trees(
-            ra_search_tree, dec_search_tree, alignment_type=PixelAlignmentType.INNER
-        ).pixel_tree
-
-    result_pixels = get_filtered_pixel_list(pixel_tree, filter_tree)
-
-    return result_pixels
+        dec_search_moc = _generate_dec_strip_moc(dec, order)
+        filter_moc = dec_search_moc
+    if ra_search_moc is not None and dec_search_moc is not None:
+        filter_moc = ra_search_moc.intersection(dec_search_moc)
+    return filter_moc
 
 
 def wrap_ra_angles(ra: np.ndarray | Iterable | int | float) -> np.ndarray:
@@ -60,7 +46,7 @@ def wrap_ra_angles(ra: np.ndarray | Iterable | int | float) -> np.ndarray:
     return np.asarray(ra, dtype=float) % 360
 
 
-def _generate_ra_strip_pixel_tree(ra_range: Tuple[float, float], order: int) -> PixelTree:
+def _generate_ra_strip_moc(ra_range: Tuple[float, float], order: int) -> MOC:
     """Generates a pixel_tree filled with leaf nodes at a given order that overlap with the ra region."""
     # Subdivide polygons (if needed) in two smaller polygons of at most 180 degrees
     division_ra = _get_division_ra(ra_range)
@@ -90,11 +76,11 @@ def _generate_ra_strip_pixel_tree(ra_range: Tuple[float, float], order: int) ->
             order,
         )
 
-    pixel_list = [HealpixPixel(order, polygon_pixel) for polygon_pixel in pixels_in_range]
-    return PixelTree.from_healpix(pixel_list)
+    orders = np.full(pixels_in_range.shape, fill_value=order)
+    return MOC.from_healpix_cells(ipix=pixels_in_range, depth=orders, max_depth=order)
 
 
-def _generate_dec_strip_pixel_tree(dec_range: Tuple[float, float], order: int) -> PixelTree:
+def _generate_dec_strip_moc(dec_range: Tuple[float, float], order: int) -> MOC:
     """Generates a pixel_tree filled with leaf nodes at a given order that overlap with the dec region."""
     nside = hp.order2nside(order)
     # Convert declination values to colatitudes, in radians, and revert their order
@@ -103,8 +89,8 @@ def _generate_dec_strip_pixel_tree(dec_range: Tuple[float, float], order: int) -
     pixels_in_range = hp.ring2nest(
         nside, hp.query_strip(nside, theta1=colat_rad[0], theta2=colat_rad[1], inclusive=True)
     )
-    pixel_list = [HealpixPixel(order, polygon_pixel) for polygon_pixel in pixels_in_range]
-    return PixelTree.from_healpix(pixel_list)
+    orders = np.full(pixels_in_range.shape, fill_value=order)
+    return MOC.from_healpix_cells(ipix=pixels_in_range, depth=orders, max_depth=order)
 
 
 def _get_division_ra(ra_range: Tuple[float, float]) -> float | None:

src/hipscat/pixel_math/cone_filter.py

@@ -1,37 +1,7 @@
-from typing import List
+import astropy.units as u
+from mocpy import MOC
 
-import healpy as hp
-import numpy as np
 
-from hipscat.pixel_math import HealpixPixel
-from hipscat.pixel_math.filter import get_filtered_pixel_list
-from hipscat.pixel_tree.pixel_tree import PixelTree
-
-
-def filter_pixels_by_cone(
-    pixel_tree: PixelTree, ra: float, dec: float, radius_arcsec: float
-) -> List[HealpixPixel]:
-    """Filter the leaf pixels in a pixel tree to return a partition_info dataframe with the pixels
-    that overlap with a cone.
-
-    Args:
-        ra (float): Right Ascension of the center of the cone in degrees
-        dec (float): Declination of the center of the cone in degrees
-        radius_arcsec (float): Radius of the cone in arcseconds
-
-    Returns:
-        List of HealpixPixels representing only the pixels that overlap with the specified cone.
-    """
-    max_order = pixel_tree.get_max_depth()
-    cone_tree = _generate_cone_pixel_tree(ra, dec, radius_arcsec, max_order)
-    return get_filtered_pixel_list(pixel_tree, cone_tree)
-
-
-def _generate_cone_pixel_tree(ra: float, dec: float, radius_arcsec: float, order: int):
-    """Generates a pixel_tree filled with leaf nodes at a given order that overlap with a cone"""
-    n_side = hp.order2nside(order)
-    center_vec = hp.ang2vec(ra, dec, lonlat=True)
-    radius_radians = np.radians(radius_arcsec / 3600.0)
-    cone_pixels = hp.query_disc(n_side, center_vec, radius_radians, inclusive=True, nest=True)
-    pixel_list = [HealpixPixel(order, cone_pixel) for cone_pixel in cone_pixels]
-    return PixelTree.from_healpix(pixel_list)
+def generate_cone_moc(ra: float, dec: float, radius_arcsec: float, order: int) -> MOC:
+    """Generate a MOC object that covers the cone"""
+    return MOC.from_cone(lon=ra * u.deg, lat=dec * u.deg, radius=radius_arcsec * u.arcsec, max_depth=order)