From e072728c90d8ac9b6ce48dab22d5f71ac218ea8d Mon Sep 17 00:00:00 2001 From: James Chiang Date: Fri, 9 Feb 2024 11:03:15 -0800 Subject: [PATCH 1/2] Use Disk to prefilter objects for PolygonalRegion selection instead of simple cuts on RA, Dec since those do not handle coordinates on the sphere well because of RA +/- 360 deg ambiguities --- skycatalogs/skyCatalogs.py | 36 +++++++++++++++++++++++++----------- 1 file changed, 25 insertions(+), 11 deletions(-) diff --git a/skycatalogs/skyCatalogs.py b/skycatalogs/skyCatalogs.py index 4c35f560..ef1890a6 100644 --- a/skycatalogs/skyCatalogs.py +++ b/skycatalogs/skyCatalogs.py @@ -8,6 +8,7 @@ import numpy as np import numpy.ma as ma from astropy import units as u +import lsst.sphgeom from skycatalogs.objects.base_object import load_lsst_bandpasses, load_roman_bandpasses from skycatalogs.utils.catalog_utils import CatalogContext from skycatalogs.objects.base_object import ObjectList @@ -103,17 +104,30 @@ def _compress_via_mask(tbl, id_column, region, source_type={'galaxy'}, time_filter = ('start_mjd' in tbl) and ('end_mjd' in tbl) and mjd is not None if region is not None: - if isinstance(region, PolygonalRegion): # special case - # Find bounding box - vertices = region.get_vertices_radec() # list of (ra, dec) - ra_max = max([v[0] for v in vertices]) - ra_min = min([v[0] for v in vertices]) - dec_max = max([v[1] for v in vertices]) - dec_min = min([v[1] for v in vertices]) - bnd_box = Box(ra_min, ra_max, dec_min, dec_max) - # Compute mask for that box - mask = _compute_region_mask(bnd_box, tbl['ra'], tbl['dec']) - if all(mask): # even bounding box doesn't intersect table rows + 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 vetices])))*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_mask(disk_region, tbl['ra'], tbl['dec']) + if all(mask): if no_obj_type_return: return None, None, None, None else: From 3c1872f1d3f05eeedb8fe0e30552e60799926aef Mon Sep 17 00:00:00 2001 From: James Chiang Date: Fri, 9 Feb 2024 11:32:38 -0800 Subject: [PATCH 2/2] fix typo and use updated function name --- skycatalogs/skyCatalogs.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/skycatalogs/skyCatalogs.py b/skycatalogs/skyCatalogs.py index ef1890a6..10051c5e 100644 --- a/skycatalogs/skyCatalogs.py +++ b/skycatalogs/skyCatalogs.py @@ -118,7 +118,7 @@ def _compress_via_mask(tbl, id_column, region, source_type={'galaxy'}, # Find largest vertex offset from the mean direction in arcsec. max_offset = np.degrees(np.arccos(min([mean_dir.dot(_) - for _ in vetices])))*3600. + for _ in vertices])))*3600. # Construct a "Disk" enclosing the polygonal region. lon_lat = lsst.sphgeom.LonLat(mean_dir) @@ -126,7 +126,7 @@ def _compress_via_mask(tbl, id_column, region, source_type={'galaxy'}, dec = lon_lat.getLat().asDegrees() disk_region = Disk(ra, dec, max_offset) - mask = _compute_mask(disk_region, tbl['ra'], tbl['dec']) + mask = _compute_region_mask(disk_region, tbl['ra'], tbl['dec']) if all(mask): if no_obj_type_return: return None, None, None, None