refactor region handling

skycatalogs/objects/gaia_object.py

@@ -19,7 +19,6 @@
 import lsst.afw.table as afwtable
 from lsst.sphgeom import HtmPixelization, Circle, LonLat
 
-from skycatalogs.utils.shapes import Disk, PolygonalRegion, compute_region_mask
 from .base_object import BaseObject, ObjectCollection
 from .base_config_fragment import BaseConfigFragment
 
@@ -178,38 +177,10 @@ def _read_fits(htm_id, gaia_config, sky_root, out_dict, logger, region=None):
     ra_full_deg = np.degrees(ra_full)
     dec_full_deg = np.degrees(dec_full)
 
-    if isinstance(region, Disk):
-        mask = compute_region_mask(region, ra_full_deg, dec_full_deg)
-        if all(mask):
-            return
-
-    elif isinstance(region, PolygonalRegion):
-        # First mask off points outside a bounding circle because it's
-        # relatively fast
-        circle = region._convex_polygon.getBoundingCircle()
-        v = circle.getCenter()
-        ra_c = LonLat.longitudeOf(v).asDegrees()
-        dec_c = LonLat.latitudeOf(v).asDegrees()
-        rad_as = circle.getOpeningAngle().asDegrees() * 3600
-        disk = Disk(ra_c, dec_c, rad_as)
-        mask = compute_region_mask(disk, ra_full_deg, dec_full_deg)
-        if all(mask):
-            return
-        if any(mask):
-            ra_compress = ma.array(ra_full, mask=mask).compressed()
-            dec_compress = ma.array(dec_full, mask=mask).compressed()
-        else:
-            ra_compress = ra_full
-            dec_compress = dec_full
-        poly_mask = compute_region_mask(region, np.degrees(ra_compress),
-                                        np.degrees(dec_compress))
-        if all(poly_mask):
-            return
-
-        # Get indices of unmasked
-        ixes = np.where(~mask)
+    mask = region.get_mask(ra_full_deg, dec_full_deg)
 
-        mask[ixes] |= poly_mask
+    if all(mask):
+        return
 
     if any(mask):
         ra_compress = ma.array(ra_full, mask=mask).compressed()
@@ -260,16 +231,7 @@ def load_collection(region, skycatalog, mjd=None, exposure=None):
         '''
         if not skycatalog:
             raise ValueError('GaiaCollection.load_collection: skycatalog cannot be None')
-        if isinstance(region, Disk):
-            ra = lsst.geom.Angle(region.ra, lsst.geom.degrees)
-            dec = lsst.geom.Angle(region.dec, lsst.geom.degrees)
-            center = lsst.geom.SpherePoint(ra, dec)
-            radius = lsst.geom.Angle(region.radius_as, lsst.geom.arcseconds)
-            refcat_region = Circle(center.getVector(), radius)
-        elif isinstance(region, PolygonalRegion):
-            refcat_region = region._convex_polygon
-        else:
-            raise TypeError(f'GaiaCollection.load_collection: {region} not supported')
+        refcat_region = region.sphgeom_region()
         if GaiaCollection.get_config() is None:
             gaia_section = skycatalog.raw_config['object_types']['gaia_star']
             GaiaCollection.set_config(gaia_section)

skycatalogs/scripts/rotate.py

@@ -11,7 +11,7 @@
 from skycatalogs.skyCatalogs import open_catalog
 from skycatalogs.utils.parquet_schema_utils import make_galaxy_schema
 from skycatalogs.utils.creator_utils import make_MW_extinction_av
-from skycatalogs.utils.shapes import Disk
+from skycatalogs.utils import Disk
 
 DC2_RA = [68, 64, 58,
           68.5, 62, 55,

skycatalogs/skyCatalogs.py

@@ -18,8 +18,6 @@
 from skycatalogs.utils.sed_tools import SsoSedFactory
 from skycatalogs.utils.sed_tools import MilkyWayExtinction
 from skycatalogs.utils.config_utils import Config
-from skycatalogs.utils.shapes import Box, Disk, PolygonalRegion
-from skycatalogs.utils.shapes import compute_region_mask
 from skycatalogs.objects.star_object import StarObject
 from skycatalogs.objects.galaxy_object import GalaxyObject
 from skycatalogs.objects.diffsky_object import DiffskyObject
@@ -28,47 +26,6 @@
 __all__ = ['SkyCatalog', 'open_catalog']
 
 
-# This function should maybe be moved to utils
-def _get_intersecting_hps(hp_ordering, nside, region):
-    '''
-    Given healpixel structure defined by hp_ordering and nside, find
-    all healpixels which intersect region, where region may be a box,
-    a disk or a polygonal region.
-    Return as some kind of iterable
-    Note it's possible extra hps which don't actually intersect the region
-    will be returned
-    '''
-    # First convert region description to an array (4,3) with (x,y,z) coords
-    # for each vertex
-    if isinstance(region, Box):
-        vec = healpy.pixelfunc.ang2vec([region.ra_min, region.ra_max,
-                                        region.ra_max, region.ra_min],
-                                       [region.dec_min, region.dec_min,
-                                        region.dec_max, region.dec_max],
-                                       lonlat=True)
-
-        pixels = healpy.query_polygon(nside, vec, inclusive=True, nest=False)
-    elif isinstance(region, Disk):
-        # Convert inputs to the types query_disk expects
-        center = healpy.pixelfunc.ang2vec(region.ra, region.dec,
-                                          lonlat=True)
-        radius_rad = (region.radius_as * u.arcsec).to_value('radian')
-
-        pixels = healpy.query_disc(nside, center, radius_rad, inclusive=True,
-                                   nest=False)
-
-    elif isinstance(region, PolygonalRegion):
-        pixels = healpy.query_polygon(nside, region.get_vertices(),
-                                      inclusive=True, nest=False)
-    else:
-        raise ValueError('Unsupported region type')
-
-    # ensure pixels are always presented in the same order
-    pixels = list(pixels)
-    pixels.sort()
-    return pixels
-
-
 def _compress_via_mask(tbl, id_column, region, source_type='galaxy',
                        mjd=None, exposure=EXPOSURE_DEFAULT):
     '''
@@ -106,47 +63,7 @@ def _compress_via_mask(tbl, id_column, region, source_type='galaxy',
     variable_filter = ('mjd' in tbl)
 
     if region is not None:
-        if isinstance(region, PolygonalRegion):
-            # Using native PolygonalRegion selection is slow, so
-            # pre-mask the RA, Dec values using an acceptance
-            # cone that encloses all of the vertices.
-
-            # Compute the mean direction from the region vertices.
-            vertices = region.get_vertices()
-            mean_dir = vertices[0]
-            for vertex in vertices[1:]:
-                mean_dir += vertex
-            mean_dir = lsst.sphgeom.UnitVector3d(mean_dir)
-
-            # Find largest vertex offset from the mean direction in arcsec.
-            max_offset = np.degrees(np.arccos(min([mean_dir.dot(_)
-                                                   for _ in vertices])))*3600.
-
-            # Construct a "Disk" enclosing the polygonal region.
-            lon_lat = lsst.sphgeom.LonLat(mean_dir)
-            ra = lon_lat.getLon().asDegrees()
-            dec = lon_lat.getLat().asDegrees()
-            disk_region = Disk(ra, dec, max_offset)
-
-            mask = compute_region_mask(disk_region, tbl['ra'], tbl['dec'])
-            if all(mask):
-                if no_obj_type_return and no_mjd_return:
-                    return None, None, None, None
-                else:    # currently if object type is returned, mjd is not
-                    return None, None, None, None, None
-
-            # Get compressed ra, dec
-            ra_compress = ma.array(tbl['ra'], mask=mask).compressed()
-            dec_compress = ma.array(tbl['dec'], mask=mask).compressed()
-            poly_mask = compute_region_mask(region, ra_compress, dec_compress)
-
-            # Get indices of unmasked
-            ixes = np.where(~mask)
-
-            # Update bounding box mask with polygonal mask
-            mask[ixes] |= poly_mask
-        else:
-            mask = compute_region_mask(region, tbl['ra'], tbl['dec'])
+        mask = region.compute_mask(tbl['ra'], tbl['dec'])
 
         if transient_filter:
             time_mask = _compute_transient_mask(mjd, tbl['start_mjd'],
@@ -512,8 +429,8 @@ def get_hps_by_region(self, region, object_type='galaxy'):
             partition = self._config['object_types'][object_type]['area_partition']
         else:
             partition = self._global_partition
-        return _get_intersecting_hps(partition['ordering'], partition['nside'],
-                                     region)
+        return region.get_intersecting_hps(partition['nside'],
+                                           partition['ordering'])
 
     def get_object_type_names(self):
         '''